植被覆盖度获取¶
植被覆盖度(Fractional Vegetation Cover,FVC)，是指植被（包括叶、茎、枝）在地面的垂直投影面积占统计区总面积的百分比，范围在 [0,100%] 之间。FVC 是刻画地表植被覆盖的重要参数，能够直观的反映一个地区绿的程度，是反应植被生长状态的重要指标，在植被变化、生态环境研究、水土保持、城市宜居等方面问题研究中起到重要作用。本案例以 Landsat-8 数据为例，计算浙江省区域的 FVC 指数。
初始化环境¶
import aie aie.Authenticate() aie.Initialize()
Landsat-8 数据检索¶
指定区域、时间、云量检索 Landsat-8 ，并对数据进行去云处理。
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)
feature_collection = aie.FeatureCollection('China_Province') \ .filter(aie.Filter.eq('province', '浙江省')) geometry = feature_collection.geometry() dataset = aie.ImageCollection('LANDSAT_LC08_C02_T1_L2') \ .filterBounds(geometry) \ .filterDate('2021-05-01', '2021-10-31') \ .filter(aie.Filter.lte('eo:cloud_cover', 30.0)) print(dataset.size().getInfo())
dataset = dataset.map(removeLandsatCloud) image = dataset.median()
计算 NDVI 指数¶
ndvi = image.normalizedDifference(['SR_B5', 'SR_B4']).rename(['NDVI']) ndvi_vis = { 'min': -0.2, 'max': 0.6, 'palette': ['#d7191c', '#fdae61', '#ffffc0', '#a6d96a', '#1a9641'] } map = aie.Map( center=ndvi.getCenter(), height=800, zoom=6 ) map.addLayer( ndvi, ndvi_vis, 'NDVI', bounds=ndvi.getBounds() ) map
定义植被覆盖度算法¶
使用像元二分模型法进行 FVC 估算。 利用 aie.Reducer.histogram 实现输入影像的直方图统计。通过 import numpy 调用数组运算，计算生长季的 NDVI 像元百分比统计中 5% 位置 NDVI 值作为土壤部分 NDVIsoil 、95% 位置的 NDVI 值作为植被部分 NDVIveg ，并通过 FVC = （NDVI - NDVIsoil）/ （NDVIveg - NDVIsoil ） 计算 FCV ，得出 FVC 。
def calculateFVC(image, scale): histogram = image.reduceRegion(aie.Reducer.histogram(2000), None, scale) histogram_info = histogram.getInfo() # print(histogram_info) #值得注意的是这里的数学运算调用了numpy的部分，所以这里不同于直接在JavaScript中集成的部分 import numpy as np bucketKey = histogram_info['NDVI_range'] bucketValue = histogram_info['NDVI_counts'] key = np.array(bucketValue) accSum = np.cumsum(key) # print(accSum[20]) # print(accSum[-1]) accPercent = accSum / accSum[-1] p5 = np.searchsorted(accPercent, 0.5) min_ndvi = bucketKey[p5 + 1] # print(min_ndvi) p95 = np.searchsorted(accPercent, 0.95) max_ndvi = bucketKey[p95] # print(max_ndvi) higher_ndvi_mask = image.gt(aie.Image(max_ndvi)) lower_ndvi_mask = image.lt(aie.Image(min_ndvi)) middle_ndvi_mask = aie.Image(1).subtract(higher_ndvi_mask).subtract(lower_ndvi_mask) tmp = image.subtract(aie.Image(min_ndvi)).divide(aie.Image(max_ndvi).subtract(aie.Image(min_ndvi))) FVC = aie.Image(1).multiply(higher_ndvi_mask).add(aie.Image(0).multiply(lower_ndvi_mask)).add(tmp.multiply(middle_ndvi_mask)) return FVC
数据可视化¶
FVC = calculateFVC(ndvi, 1000) vis_params = { 'min': 0, 'max': 1, 'palette': [ '#a1a1a1', '#008000' ] } map.addLayer( FVC, vis_params, 'fvc', bounds=ndvi.getBounds() ) map
