決策樹模型實現(xiàn)黔東南州水稻提取
依據(jù)作物在不同物候期內(nèi)衛(wèi)星影像的光譜存在差異的特征和地形因子，可建立水稻提取算法，進行水稻提取。
初始化環(huán)境
import aie
aie.Authenticate()
aie.Initialize()
指定需要檢索的區(qū)域
feature_collection = aie.FeatureCollection('China_City') \
.filter(aie.Filter.eq('city', '黔東南苗族侗族自治州'))
region = feature_collection.geometry()
DEM處理
# 指定檢索數(shù)據(jù)集，可設(shè)置檢索的空間范圍
elevation = aie.ImageCollection('JAXA_ALOS_AW3D30_V3_2') \
.filterBounds(region)\
.select(['DSM'])\
.mosaic()\
.clip(region)
map = aie.Map(
center=feature_collection.getCenter(),
height=800,
zoom=7
)
vis_params = {
'bands': 'DSM',
'min': 100,
'max': 2200,
'palette': [
'#0000ff', '#00ffff', '#ffff00', '#ff0000', '#ffffff'
]
}
map.addLayer(
elevation,
vis_params,
'Elevation',
bounds=elevation.getBounds()
)
map
黔東南州高程分布圖
task = aie.Export.image.toAsset(elevation,'dem_qdn',30)
task.start()
## 坡度，下載aie還不能計算，我這里使用ArcGIS運算
slope = aie.Image('user/c7ec068793e54fccb9ba8692ed9d0b91').clip(region)
vis_params = {
'min': 0,
'max': 80,
'palette': [
'#0000ff', '#00ffff', '#ffff00', '#ff0000', '#ffffff'
]
}
map.addLayer(
slope,
vis_params,
'slope',
bounds=slope.getBounds()
)
map
黔東南州坡度分布圖
這里說一下，我的谷歌瀏覽器上傳影像數(shù)據(jù)，不知道為什么失敗，我是用edge瀏覽器上傳的，我把我的谷歌瀏覽器版本上傳給了官方。
Landsat 8 數(shù)據(jù)處理
# 插秧期影像
# 指定檢索數(shù)據(jù)集，可設(shè)置檢索的空間和時間范圍，以及屬性過濾條件（如云量過濾等）
img1 = aie.ImageCollection('LANDSAT_LC08_C02_T1_L2') \
.filterBounds(region) \
.filterDate('2021-4-01', '2021-6-10') \
.filter(aie.Filter.lte('eo:cloud_cover', 20.0))\
.median()\
.clip(region)
# print(img1.size().getInfo())
vis_params = {
'bands': ['SR_B4', 'SR_B3', 'SR_B2'],
'min': 8000,
'max': 13000,
}
map.addLayer(
img1,
vis_params,
'img1',
bounds=img1.getBounds()
)
map
# 生長期影像
img2 = aie.ImageCollection('LANDSAT_LC08_C02_T1_L2') \
.filterBounds(region) \
.filterDate('2021-6-20', '2021-8-30') \
.filter(aie.Filter.lte('eo:cloud_cover', 45.0))\
print(img2.size().getInfo())
img2 = img2.median().clip(region)
vis_params = {
'bands': ['SR_B4', 'SR_B3', 'SR_B2'],
'min': 8000,
'max': 13000,
}
map.addLayer(
img2,
vis_params,
'img2',
bounds=img2.getBounds()
)
map
NDVI 計算
# NDVI擴大10，好比較
NDVI1 = img1.normalizedDifference(['SR_B5', 'SR_B4'])\
.multiply(aie.Image.constant(10)).rename(['NDVI'])
NDVI2 = img2.normalizedDifference(['SR_B5', 'SR_B4'])\
.multiply(aie.Image.constant(10)).rename(['NDVI'])
NDVI_diff = NDVI2.subtract(NDVI1).rename(['Diff'])
import numpy as np
scale = 1000
histogram = NDVI1.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]
p2 = np.searchsorted(accPercent, 0.2)
min_ndvi = bucketKey[p2 + 1]
print('min_ndvi1:%f' % min_ndvi)
p98 = np.searchsorted(accPercent, 0.98)
max_ndvi = bucketKey[p98]
print('max_ndvi1:%f' % max_ndvi)
查看NDVI
提取規(guī)則
# 水稻提取規(guī)則集
## 水稻一般生長在海拔900m以下，坡度在20度以下
mask1 = elevation.lt(aie.Image.constant(900)).clip(region)
mask2 = slope.lt(aie.Image.constant(20)).clip(region)
## 水稻播種期NDVI一般在0.32至0.38，每個地方可能有差異
mask3 = NDVI1.gt(aie.Image.constant(3.2)).And(NDVI1.lt(aie.Image.constant(3.8)))
## 水稻生長期NDVI和播種期NDVI一般在-0.9至0.6，每個地方可能有差異
mask4 = NDVI_diff.gt(aie.Image.constant(-0.9)).And(NDVI_diff.lt(aie.Image.constant(0.6)))
rice = mask1.And(mask2).And(mask3).And(mask4)
mask_vis = {
'min': 0,
'max': 1,
'palette': ['#ffffff', '#008000'] # 0:白色, 1:綠色
}
map.addLayer(rice,mask_vis, 'wheat', bounds=region.getBounds()) # 綠色區(qū)域為水稻
水稻提取結(jié)果
task = aie.Export.image.toAsset(rice,'rice_extract',30)
task.start()
精度評價
這一部分在ArcGIS和Excel里面完成，查找統(tǒng)計年鑒可知黔東南州水稻種植面積S1=233 萬畝,提取出來的面積S2= 252萬畝，提取下來結(jié)果如表所示，總的來說還是比較粗糙，希望大家有更好的算法。
模型構(gòu)建器
提取結(jié)果
