AI Earth ——开发者模式案例7：植被覆盖度提取

此星光明

发布于 2024-02-02 10:04:09

1850

发布于 2024-02-02 10:04:09

植被覆盖度获取¶

植被覆盖度(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