植被覆盖度(Fractional Vegetation Cover,FVC),是指植被(包括叶、茎、枝)在地面的垂直投影面积占统计区总面积的百分比,范围在 [0,100%] 之间。FVC 是刻画地表植被覆盖的重要参数,能够直观的反映一个地区绿的程度,是反应植被生长状态的重要指标,在植被变化、生态环境研究、水土保持、城市宜居等方面问题研究中起到重要作用。本案例以 Landsat-8 数据为例,计算浙江省区域的 FVC 指数。
import aie
aie.Authenticate()
aie.Initialize()
指定区域、时间、云量检索 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 = 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)
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.05)
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