植被覆蓋度獲取

植被覆蓋度(Fractional Vegetation Cover,FVC)，是指植被（包括葉、莖、枝）在地面的垂直投影面積占統(tǒng)計(jì)區(qū)總面積的百分比，范圍在 [0,1] 之間。FVC 是刻畫地表植被覆蓋的重要參數(shù)，能夠直觀的反映一個(gè)地區(qū)綠的程度，是反應(yīng)植被生長(zhǎng)狀態(tài)的重要指標(biāo)，在植被變化、生態(tài)環(huán)境研究、水土保持、城市宜居等方面問題研究中起到重要作用。本案例以 Landsat-8 數(shù)據(jù)為例，計(jì)算貴州省區(qū)域的 FVC 指數(shù)。

初始化環(huán)境

import aie

aie.Authenticate()
aie.Initialize()

Landsat-8 數(shù)據(jù)檢索

指定區(qū)域、時(shí)間、云量檢索 Landsat-8 ，并對(duì)數(shù)據(jù)進(jìn)行去云處理。

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)

計(jì)算 NDVI 指數(shù)

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

定義植被覆蓋度算法

使用像元二分模型法進(jìn)行 FVC 估算。 利用 aie.Reducer.histogram 實(shí)現(xiàn)輸入影像的直方圖統(tǒng)計(jì)。通過 numpy 調(diào)用數(shù)組運(yùn)算，計(jì)算生長(zhǎng)季的 NDVI 像元百分比統(tǒng)計(jì)中 5% 位置 NDVI 值作為土壤部分 NDVIsoil 、95% 位置的 NDVI 值作為植被部分 NDVIveg ，并通過 FVC = （NDVI - NDVIsoil）/ （NDVIveg - NDVIsoil ） 計(jì)算 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

數(shù)據(jù)可視化

FVC = calculateFVC(ndvi, 1000)

vis_params = {
    'min': 0,
    'max': 1,
    'palette': [
        '#a1a1a1', '#008000'
    ]
}

map.addLayer(
    FVC,
    vis_params,
    'fvc',
    bounds=ndvi.getBounds()
)
map

導(dǎo)出數(shù)據(jù)

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

貴州FVC原始數(shù)據(jù)

后記

AIE進(jìn)行遙感云計(jì)算的時(shí)候還是很方便，可能剛剛出來，很多地方還是需要完善，這個(gè)案例里面，我導(dǎo)出數(shù)據(jù)以后要到ArcGIS里面再出來一下下。接下來，我利用自然間斷法分成了五類，然后再統(tǒng)計(jì)這五類的面積，這ArcGIS操作都很簡(jiǎn)單了，這里就不多說，還有就是阿里云的小哥哥特別有耐心，特別負(fù)責(zé)任，計(jì)算也很強(qiáng)。