决策树模型实现冬小麦提取
决策树模型实现冬小麦提取
依据作物在不同物候期内卫星影像的光谱存在差异的特征，可建立冬小麦提取算法，进行像元尺度冬小麦提取。
初始化环境
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来看看结果。
这里可以计算一下面积，然后和贵州省冬小麦的播种的实际面积进行比较。
A
I
E
提取面积
S
1
=
257356
∗
100
∗
100
/
10000
=
257356
h
a
=
257
k
h
a
AIE提取面积S1= 257356 * 100 *100 / 10000 =257356 ha = 257 kha
AIE提取面积S1=257356∗100∗100/10000=257356ha=257kha
P
I
E
提取面积
S
2
=
114
k
h
a
PIE提取面积S2= 114kha
PIE提取面积S2=114kha
统计年鉴播种面积
S
3
=
138.05
k
h
a
统计年鉴播种面积S3= 138.05kha
统计年鉴播种面积S3=138.05kha
以S2为真值，则
∣
S
1
−
S
3
∣
S
3
=
∣
257
−
138
∣
138
=
86.2
%
frac{left | S1-S3 right |}{S3} = frac{left | 257-138 right |}{138}=86.2%
S3∣S1−S3∣=138∣257−138∣=86.2%
∣
S
2
−
S
3
∣
S
3
=
∣
114
−
138
∣
138
=
17.3
%
frac{left | S2-S3 right |}{S3} = frac{left | 114-138 right |}{138}=17.3%
S3∣S2−S3∣=138∣114−138∣=17.3%
总结一下吧，我也不是比较，完全没有黑的意思，也没有踩和贬的意思。就是一样的数据，差不多的代码，跑下来结果相差有一点点大。原因的话，AIE里面没有去云，没有分类后处理，还有就是AIE里面我运行下来，空值较多，当然这也是和我的搜索条件有关。
本案例主要引用了AIE和PIE里面的案例，然后自己修改的。