基于AIE的贵州省FVC提取

植被覆盖度获取

植被覆盖度(Fractional Vegetation Cover,FVC)，是指植被（包括叶、茎、枝）在地面的垂直投影面积占统计区总面积的百分比，范围在 [0,1] 之间。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()

裁剪影像

image = image.clip(geometry)

计算 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 实现输入影像的直方图统计。通过 numpy 调用数组运算，计算生长季的 NDVI 像元百分比统计中 5% 位置 NDVI 值作为土壤部分 NDVIsoil 、95% 位置的 NDVI 值作为植被部分 NDVIveg ，并通过 FVC = （NDVI - NDVIsoil）/ （NDVIveg - NDVIsoil ） 计算 FCV ，得出 FVC 。

import numpy as np
import pandas as pd

def calculateFVC(image, scale):
    histogram = image.reduceRegion(aie.Reducer.histogram(2000), None, scale)
    histogram_info = histogram.getInfo()
    # print(histogram_info)


    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

导出数据

task = aie.Export.image.toAsset(FVC, 'FVC_export_result', 100)
task.start()

后记

AIE进行遥感云计算的时候还是很方便，可能刚刚出来，很多地方还是需要完善，这个案例里面，我导出数据以后要到ArcGIS里面再出来一下下。接下来，我利用自然间断法分成了五类，然后再统计这五类的面积，这ArcGIS操作都很简单了，这里就不多说，还有就是阿里云的小哥哥特别有耐心，特别负责任，计算也很强。本案例主要引用AIE的官方例子。