AI Earth ——开发者模式案例3：典型植被指数计算及区域统计

典型植被指数计算及区域统计¶

对检索的影像（以 Landsat-8 为例），通过波段运算计算常见的指数。并以归一化植被指数（ NDVI ）为例，进行区域均值统计以及时序折线图制作。

初始化环境¶

import aie
 
aie.Authenticate()
aie.Initialize()

典型光谱指数算法¶

定义典型指数计算方法。使用 aie.Image.add 、 aie.Image.subtract 、 aie.Image.multiply 和 aie.Image.divide 实现影像波段运算。另外可使用 aie.Image.normalizedDifference 实现两个波段的归一化差值运算 (Band1-Band2)/(Band1+Band2) ，使用 aie.Image.expression 可实现构建表达式对影像进行波段运算。

如切换卫星数据源，需要调整对应的波段名称。

# 比值植被指数
def getRVI(image):
    nir = image.select(['SR_B5'])
    red = image.select(['SR_B4'])
    rvi = nir.divide(red)
    return rvi.rename(['RVI'])
 
# 增强型植被指数
def getEVI(image):
    evi = image.expression(
        '(2.5 * (nir - red)) /(nir + 6 * red - 7.5 * blue + 1)', 
        {
            'nir': image.select(['SR_B5']),
            'red': image.select(['SR_B4']),
            'blue': image.select(['SR_B2'])
    }).rename('EVI')
    return evi
 
# 归一化植被指数
def getNDVI(image):
    ndvi = image.normalizedDifference(['SR_B5', 'SR_B4']).rename(['NDVI'])
    return ndvi
 
# 近红外植被指数
def getNIRv(image):
    nir = image.select(['SR_B5'])
    nirv = nir.multiply(image.normalizedDifference(['SR_B5', 'SR_B4'])).rename('NIRv')
    return nirv                                                 
  
# 土壤调整植被指数
def getSAVI(image):
    nir = image.select(['SR_B5'])
    red = image.select(['SR_B4'])                                                
    savi = ((nir.subtract(red)).multiply(aie.Image.constant(1.5))).divide((nir.add(red)).add(aie.Image.constant(0.5))).rename('SAVI')
    return savi
 
# 归一化水体指数
def getNDWI(image):
    ndwi = image.normalizedDifference(['SR_B3', 'SR_B5']).rename('NDWI')
    return ndwi

Landsat-8 数据检索¶

指定区域、时间、云量检索 Landsat-8 ，并对数据进行去云处理。

region = aie.FeatureCollection('China_Province') \
            .filter(aie.Filter.eq('province', '浙江省')) \
            .geometry()

def l8Collection(startdate, enddate):
    images = aie.ImageCollection('LANDSAT_LC08_C02_T1_L2') \
            .filterBounds(region) \
            .filterDate(startdate, enddate)
    return images

def removeLandsatCloud(image):
    cloudShadowBitMask = (1 << 4)
    cloudsBitMask = (1 << 3)
    qa = image.select('QA_PIXEL')
    mask = qa.bitwiseAnd(aie.Image(cloudShadowBitMask)).eq(aie.Image(0)).And(qa.bitwiseAnd(aie.Image(cloudsBitMask)).eq(aie.Image(0)))
    return image.updateMask(mask)

lc8_collection = l8Collection('2021-08-01', '2021-08-31')
lc8_collection.map(removeLandsatCloud)
print(lc8_collection.size().getInfo())
lc8_img = lc8_collection.max()

NDVI 计算及可视化¶

以 NDVI 计算为例输出指数计算成果，并地图可视化展示。

ndvi = getNDVI(lc8_img)
 
ndvi_vis = {
    'min': -0.2,
    'max': 0.6,
    'palette': [
        '#2B83BA', '#ABDDA4', '#FFFFBF', '#FDAE61', '#D7191C'
    ]
}
 
map = aie.Map(
    center=ndvi.getCenter(),
    height=800,
    zoom=5
)
 
map.addLayer(
    ndvi,
    ndvi_vis,
    'NDVI',
    bounds=ndvi.getBounds()
)
 
map

NDVI 区域统计¶

使用中国市级行政区划数据，按照市域范围对 NDVI 进行均值统计。使用 aie.Image.reduceRegions 和 aie.Reducer.mean 实现对影像进行指定区域范围均值统计。 当在较大范围内执行 ReduceRegion 或者 ReduceRegions 函数时，可能存在较为耗时的情况。开发者根据实际需求调整 scale（单位：米），scale 越大，耗时越少。

通过引用 Python 的 pyplot 绘制浙江各地市区域 NDVI 均值统计图。

zone = aie.FeatureCollection('China_City') \
         .filter(aie.Filter.eq('province', '浙江省'))
 
zone_mean = ndvi.reduceRegions(zone, aie.Reducer.mean(), 1000)
zone_info = zone_mean.getInfo()
 
x_list = []
y_list = []
for feature in zone_info['features']:
    x_list.append(feature['properties']['city'])
    y_list.append(feature['properties']['NDVI_mean'])
    
# print(x_list)
# print(y_list)
 
from bqplot import pyplot as plt
plt.figure(1, title='2021年浙江省各市NDVI均值统计')
plt.bar(x_list, y_list)   #colors=['MediumSeaGreen']
plt.show()

杭州市宁波市温州市嘉兴市湖州市绍兴市金华市衢州市舟山市台州市丽水市00.020.040.060.080.10.120.140.160.180.20.220.240.260.282021年浙江省各市NDVI均值统计

NDVI时间序列分析¶

在指定空间范围内实现时间序列统计分析，并绘制折线图。

def doSeries(start_time, end_time, zone):
    lc8_col = l8Collection(start_time, end_time)
    lc8_col.map(removeLandsatCloud)
    lc8_img = lc8_col.mosaic()
    ndvi = getNDVI(lc8_img)
    return ndvi.reduceRegion(aie.Reducer.mean(), zone, 1000)
 
zone = aie.FeatureCollection('China_City') \
          .filter(aie.Filter.eq('province', '浙江省')) \
          .geometry()
 
x_ndvi_series = []
y_ndvi_series = []
 
year = '2021'
mon = ['01','02','03','04','05','06','07','08','09','10','11','12']
lday = ['31','28','31','30','31','30','31','31','30','31','30','31']
 
for i in range(0,12):
    startdate = year + '-' + mon[i] + '-01' 
    enddate = year + '-' + mon[i] + '-' + lday[i]
    
    lc8_ndvi_mon = doSeries(startdate, enddate , zone)
    x_ndvi_series.append(mon[i] + '月')
    y_ndvi_series.append(lc8_ndvi_mon.getInfo()['NDVI_mean'])
 
 
# print(x_ndvi_series)
# print(y_ndvi_series)
 
 
from bqplot import pyplot as plt
plt.figure(2, title='2021年浙江省逐月NDVI均值统计')
plt.plot(x_ndvi_series, y_ndvi_series)
plt.show()

01月02月03月04月05月06月07月08月09月10月11月12月0.060.080.10.120.140.160.180.20.220.242021年浙江省逐月NDVI均值统计

影像输出¶

task = aie.Export.image.toAsset(ndvi, 'NDVI', 30)
task.start()

参考文献：¶

Zeng, Y., Hao, D., Huete, A. et al. Optical vegetation indices for monitoring terrestrial ecosystems globally. Nat Rev Earth Environ 3, 477–493 (2022). Optical vegetation indices for monitoring terrestrial ecosystems globally | Nature Reviews Earth & Environment