决策树模型实现黔东南州水稻提取
依据作物在不同物候期内卫星影像的光谱存在差异的特征和地形因子，可建立水稻提取算法，进行水稻提取。
初始化环境
import aie
aie.Authenticate()
aie.Initialize()指定需要检索的区域
feature_collection = aie.FeatureCollection('China_City') \
.filter(aie.Filter.eq('city', '黔东南苗族侗族自治州'))
region = feature_collection.geometry()DEM处理
# 指定检索数据集，可设置检索的空间范围
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
这里说一下，我的谷歌浏览器上传影像数据，不知道为什么失败，我是用edge浏览器上传的，我把我的谷歌浏览器版本上传给了官方。
Landsat 8 数据处理
# 插秧期影像
# 指定检索数据集，可设置检索的空间和时间范围，以及属性过滤条件（如云量过滤等）
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)
提取规则
# 水稻提取规则集
## 水稻一般生长在海拔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()) # 绿色区域为水稻
task = aie.Export.image.toAsset(rice,'rice_extract',30)
task.start()精度评价
这一部分在ArcGIS和Excel里面完成，查找统计年鉴可知黔东南州水稻种植面积S~1~=233 万亩,提取出来的面积S~2~= 252万亩，提取下来结果如表所示，总的来说还是比较粗糙，希望大家有更好的算法。
$$总体误差= \frac{\left | S2-S1 \right |}{S1} = \frac{\left | 252-233 \right |}{233}=8.1\%$$
