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

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

初始化环境

In [1]:
import aie

aie.Authenticate()
aie.Initialize()

典型光谱指数算法

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

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

In [2]:
# 比值植被指数
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 ,并对数据进行去云处理。

In [3]:
region = aie.FeatureCollection('China_Province') \
            .filter(aie.Filter.eq('province', '浙江省'))
In [4]:
def l8Collection(startdate, enddate):
    images = aie.ImageCollection('LANDSAT_LC08_C02_T1_L2') \
            .filterBounds(region) \
            .filterDate(startdate, enddate)
    return images
In [5]:
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)
In [6]:
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 计算为例输出指数计算成果,并地图可视化展示。

In [7]:
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.reduceRegionsaie.Reducer.mean 实现对影像进行指定区域范围均值统计。 当在较大范围内执行 ReduceRegion 或者 ReduceRegions 函数时,可能存在较为耗时的情况。开发者根据实际需求调整 scale单位:米),scale 越大,耗时越少。

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

In [8]:
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()

NDVI时间序列分析

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

In [9]:
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_Province') \
          .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()

影像输出

In [10]:
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). https://doi.org/10.1038/s43017-022-00298-5