決策樹模型實現(xiàn)冬小麥提取

依據作物在不同物候期內衛(wèi)星影像的光譜存在差異的特征，可建立冬小麥提取算法，進行像元尺度冬小麥提取。

初始化環(huán)境

import aie

aie.Authenticate()
aie.Initialize()

指定需要檢索的區(qū)域

feature_collection = aie.FeatureCollection('China_Province') \
                        .filter(aie.Filter.eq('province', '貴州省'))
region = feature_collection.geometry()

去云

因為SENTINEL_MSIL2A數(shù)據集現(xiàn)在AIE還沒有去云波段，所以這一步目前還不能做。

影像檢索

貴州省冬小麥典型物候期。播種期 10中上旬-11月上旬，旺長期2月上旬-3月，成熟期4月下旬-5月上旬

# 播種期影像
img1 = aie.ImageCollection('SENTINEL_MSIL2A') \
             .filterBounds(region) \
             .filterDate('2021-10-01', '2021-11-11') \
             .filter(aie.Filter.lte('eo:cloud_cover',60.0)) \
             .select(["B11","B8","B4","B3","B2"])\
             .median()
             # .map(removeLandsatCloud)          

# 拔節(jié)期影像
img2 = aie.ImageCollection('SENTINEL_MSIL2A')\
    .filterDate("2021-02-01", "2021-03-01")\
    .filterBounds(region)\
    .filter(aie.Filter.lt('eo:cloud_cover', 60))\
    .select(["B11","B8","B4","B3"])\
    .median()

# 成熟收獲期影像
img3 = aie.ImageCollection('SENTINEL_MSIL2A')\
    .filterDate("2021-04-20", "2021-05-10")\
    .filterBounds(region)\
    .filter(aie.Filter.lt('eo:cloud_cover', 60))\
    .select(["B11","B8","B4","B3","B2",])\
    .median()

波段提取

red1 = img1.select("B4")
nir1 = img1.select("B8")
swir1 = img1.select("B11")
red2 = img2.select("B4")
nir2 = img2.select("B8")
red3 = img3.select("B4")
nir3 = img3.select("B8")

ndvi1 = (nir1.subtract(red1)).divide(nir1.add(red1)).rename(["NDVI"]).select("NDVI")
nbr1 = (nir1.subtract(swir1)).divide(nir1.add(swir1)).rename(["NBR"]).select("NBR")
ndvi2 = (nir2.subtract(red2)).divide(nir2.add(red2)).rename(["NDVI"]).select("NDVI")
ndvi3 = (nir3.subtract(red3)).divide(nir3.add(red3)).rename(["NDVI"]).select("NDVI")

# 小麥在10月份的近紅外波段更大，短波紅外波段更小
# 條件1：播種期NDVI小，NBR小
# 條件2：拔節(jié)期抽穗期 NDVI大
# 條件3：成熟期NDVI小于拔節(jié)期NDVI

wheat = ndvi1.lt(aie.Image.constant(0.3))\
            .And(nbr1.lt(aie.Image.constant(0.07)))\
            .And(ndvi2.gt(aie.Image.constant(0.32)))\
            .And(ndvi3.lt(ndvi2))\
            .clip(region)

數(shù)據可視化

之前我運行過了，等一下我們之間看本地效果。

# 結果可視化

map = aie.Map(
    center=region.getCenter(),
    height=800,
    zoom=7
)

vis_params = {
    'color': '#00FF00'
}

map.addLayer(
    region,
    vis_params,
    'region',
    bounds=region.getBounds()
)

mask_vis  = {
    'min': 0,
    'max': 1,
    'palette': ['#ffffff', '#008000']    # 0:白色, 1:綠色
}

ndvi_vis  = {
    'min': -0.2,
    'max': 0.6,
    'palette': ['#d7191c', '#fdae61', '#ffffc0', '#a6d96a', '#1a9641']
}

map.addLayer(ndvi1,ndvi_vis, 'ndvi1', bounds=region.getBounds())
map.addLayer(nir1,ndvi_vis, 'nir1', bounds=region.getBounds())
map.addLayer(ndvi2,ndvi_vis, 'ndvi2', bounds=region.getBounds())
map.addLayer(ndvi3,ndvi_vis, 'ndvi3', bounds=region.getBounds())
map.addLayer(wheat,mask_vis, 'wheat', bounds=region.getBounds())    # 綠色區(qū)域為小麥

map

重分類

wheat = wheat.where(wheat.eq(aie.Image.constant(0)),aie.Image(0))\
             .where(wheat.eq(aie.Image.constant(1)),aie.Image(1))

導出數(shù)據

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

精度評價

由于分類后處理的很多函數(shù)，aie都還沒有，所以，可以去ArcGIS來看看結果。

貴州省冬小麥提取結果

這里可以計算一下面積，然后和貴州省冬小麥的播種的實際面積進行比較。

面積

貴州統(tǒng)計年鑒

pie提取結果

總結一下吧，我也不是比較，完全沒有黑的意思，也沒有踩和貶的意思。就是一樣的數(shù)據，差不多的代碼，跑下來結果相差有一點點大。原因的話，AIE里面沒有去云，沒有分類后處理，還有就是AIE里面我運行下來，空值較多，當然這也是和我的搜索條件有關。
本案例主要引用了AIE和PIE里面的案例，然后自己修改的。