对检索的影像(以 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 ,并对数据进行去云处理。
region = aie.FeatureCollection('China_Province') \
.filter(aie.Filter.eq('province', '浙江省'))
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 = 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 进行均值统计。使用 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()
在指定空间范围内实现时间序列统计分析,并绘制折线图。
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()
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