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

import aie

aie.Authenticate()
aie.Initialize()

# 比值植被指数
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

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 = 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

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_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()

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