决策树模型实现冬小麦提取
依据作物在不同物候期内卫星影像的光谱存在差异的特征，可建立冬小麦提取算法，进行像元尺度冬小麦提取。
初始化环境
import aie
aie.Authenticate()
aie.Initialize()
指定需要检索的区域
feature_collection = aie.FeatureCollection('China_Province') \
.filter(aie.Filter.eq('province', '贵州省'))
region = feature_collection.geometry()
去云
因为SENTINEL_MSIL2A数据集现在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)
# 拔节期影像
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：拔节期抽穗期 NDVI大
# 条件3：成熟期NDVI小于拔节期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)
数据可视化
之前我运行过了，等一下我们之间看本地效果。
# 结果可视化
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()) # 绿色区域为小麦
map
重分类
wheat = wheat.where(wheat.eq(aie.Image.constant(0)),aie.Image(0))\
.where(wheat.eq(aie.Image.constant(1)),aie.Image(1))
导出数据
task = aie.Export.image.toAsset(wheat,'wheat',100)
task.start()
精度评价
由于分类后处理的很多函数，aie都还没有，所以，可以去ArcGIS来看看结果。
贵州省冬小麦提取结果
这里可以计算一下面积，然后和贵州省冬小麦的播种的实际面积进行比较。
面积
贵州统计年鉴
pie提取结果
以S2为真值，则
总结一下吧，我也不是比较，完全没有黑的意思，也没有踩和贬的意思。就是一样的数据，差不多的代码，跑下来结果相差有一点点大。原因的话，AIE里面没有去云，没有分类后处理，还有就是AIE里面我运行下来，空值较多，当然这也是和我的搜索条件有关。
本案例主要引用了AIE和PIE里面的案例，然后自己修改的。
网友评论