GEE python:geemap进行ndvi和ndwi计算并设定指定的阈值提取土地类型

安装地球引擎API和geemap
安装地球引擎的Python API和geemap。geemap Python包是建立在ipyleaflet和folium包之上的，它实现了几个与地球引擎数据层交互的方法，比如Map.addLayer()、Map.setCenter()和Map.centerObject()。下面的脚本检查geemap包是否已经安装。如果没有，它将安装geemap，它会自动安装其依赖项，包括earthengine-api、folium和ipyleaflet。

# Installs geemap package
import subprocess

try:
    import geemap
except ImportError:
    print('Installing geemap ...')
    subprocess.check_call(["python", '-m', 'pip', 'install', 'geemap'])

按照惯例依旧进行地图的加载来确保可以进行正常的使用，这里主要的就是进行GEE的授权。

Map = geemap.Map(center=[40,-100], zoom=4)
Map

​我们都知道NDVI可以进行波段运算获得，最简单的方式就是通过iGEE中的函数来获取，页可以通过四则运算来获取，我们进行对比分析看下，下面这个是归一化植被指数计算公式
**landsat2008.normalizedDifference(*args, kwargs)

Computes the normalized difference between two bands. If the bands to use
are not specified, uses the first two bands. The normalized difference is
computed as (first − second) / (first + second). Note that the returned
image band name is ‘nd’, the input image properties are not retained in the
output image, and a negative pixel value in either input band will cause
the output pixel to be masked. To avoid masking negative input values, use
ee.Image.expression() to compute normalized difference.

Args:
input: The input image.
bandNames: A list of names specifying the bands to use. If
not specified, the first and second bands are used.

https://www.cbedai.net/xg

代码：

# Add Earth Engine dataset
# Load two 5-year Landsat 7 composites.
landsat1999 = ee.Image('LANDSAT/LE7_TOA_5YEAR/1999_2003')
landsat2008 = ee.Image('LANDSAT/LE7_TOA_5YEAR/2008_2012')

# Compute NDVI the hard way.
ndvi1999 = landsat1999.select('B4').subtract(landsat1999.select('B3')) \
               .divide(landsat1999.select('B4').add(landsat1999.select('B3')))

# Compute NDVI the easy way.
ndvi2008 = landsat2008.normalizedDifference(['B4', 'B3'])

# Compute the multi-band difference image.
diff = landsat2008.subtract(landsat1999)
Map.addLayer(diff,
             {'bands': ['B4', 'B3', 'B2'], 'min': -32, 'max': 32},
             'difference')

# Compute the squared difference in each band.
squaredDifference = diff.pow(2)
Map.addLayer(squaredDifference,
             {'bands': ['B4', 'B3', 'B2'], 'max': 1000},
             'squared diff.')

最后进行结果展示：

Map.addLayerControl() # This line is not needed for ipyleaflet-based Map.
Map

这里我们进行指定阈值范围的加载，设定ndvi小于0.2，ndwi小于0的部分定义为裸地，代码如下：

# Add Earth Engine dataset
# Load a Landsat 8 image.
image = ee.Image('LANDSAT/LC08/C01/T1_TOA/LC08_044034_20140318')

# Create NDVI and NDWI spectral indices.
ndvi = image.normalizedDifference(['B5', 'B4'])
ndwi = image.normalizedDifference(['B3', 'B5'])

# Create a binary layer using logical operations.
bare = ndvi.lt(0.2).And(ndwi.lt(0))

# Mask and display the binary layer.
Map.setCenter(-122.3578, 37.7726, 12)
Map.setOptions('satellite')
Map.addLayer(bare.updateMask(bare), {}, 'bare')

最终加载裸地的结果：