基於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官方案例。