PIE-Engine 教程:水稻面积提取4—阈值筛选(宿迁市)

这里水稻面积提取用阈值法进行，这里我们需要用到的是将随机森林和svm支持向量机的监督分类分析的结果用mask去除和提取，这里我们所需要用到几个函数

where(condition,value)

执行值的有条件替换。

方法参数：

- image(Image)

Image实例。

- condition(Image|Double)

测试图像。

- value(Image|Double)

测试不为零时使用的输出值。

返回值：Image

水体指数

LSWI=（NIR-SWIR）/（NIR+SWIR式）中，NIR和SWIR分别表示MODIS近红外和短波红外波段处的地表反射率；
2）、提取水稻像元：利用水指数LSWI与植被指数EVI的关系，提取水稻像元，其过程为：在灌水移栽期LSWI >0.12，EVI <0.26，（LSWI+0.05）>EVI，且其后第6到11个8天合成的图像的平均EVI>0.35；

代码：

/**
 * @Name    :   基于 PIE Engine Studio 的水稻自动提取
 * @Author  :   武汉大学VHR队
 * @Desc    :   淮安市水稻提取 - 阈值筛选
 * @Source  :   航天宏图第四届 “航天宏图杯”PIE软件二次开发大赛云开发组三等奖获奖作品
 */

//4.阈值筛选
//加载淮安市矢量
var roi = pie.FeatureCollection("NGCC/CHINA_CITY_BOUNDARY").filter(pie.Filter.eq('name', '淮安市')).first().geometry()
Map.centerObject(roi, 8)

//导入监督分类结果
var rf = pie.Image('user/pieadmin/saas/06/rf').select('B1');
var svm = pie.Image('user/pieadmin/saas/06/svm').select('B1');

//可视化参数
var visParam = {
    opacity: 1,
    min: 0,
    max: 1,
    palette: '00FFFF, 0000FF'
};
//设置mask，提取监督分类结果中水稻像素，并将两种分类结果取并集
var mask = rf.lt(2)
var rf_rice = mask.where(mask.eq(0), 0)
var rf_ricenew = mask.updateMask(mask)
var mask2 = svm.lt(2)
var svm_rice = mask2.where(mask2.eq(0), 0)
var svm_ricenew = mask2.updateMask(mask2)
//并集
var symm = rf_rice.add(svm_rice)
var symmetricDifference = symm.gt(0)
symmetricDifference = symmetricDifference.updateMask(symmetricDifference)

//计算阈值
function VI(image) {
    var lswi = image.select('B1');
    var evi = image.select('B2');
    var vi = evi.subtract(lswi);
    return image.addBands(vi.rename('VI'));
}

function negative(image) {
    var lswi = image.select('B1');
    var evi = image.select('B2');
    var negvi = lswi.subtract(evi);
    return image.addBands(negvi.rename('negVI'));
}
//doy:14 45 74 105 135 166 196 227 258 288 319 349
//选择位于水稻插秧期的月度影像
var col1 = pie.ImageCollection.fromImages([
    pie.Image("user/pieadmin/saas/06/s2_huaian_apr"),
    pie.Image("user/pieadmin/saas/06/s2_huaian_may"),
    pie.Image("user/pieadmin/saas/06/s2_huaian_june"),
    pie.Image("user/pieadmin/saas/06/s2_huaian_july")
])
col1 = col1.select(['B1', 'B2']).map(negative)
var apr = col1.getAt(0).addBands(pie.Image().constant(105).rename('doy'))
var may = col1.getAt(1).addBands(pie.Image().constant(135).rename('doy'))
var june = col1.getAt(2).addBands(pie.Image().constant(166).rename('doy'))
var july = col1.getAt(3).addBands(pie.Image().constant(196).rename('doy'))
var transplant = pie.ImageCollection.fromImages([apr, may, june, july]) //B1,B2,negVI,doy

//选择位于水稻生长期的月度影像                                    
var col2 = pie.ImageCollection.fromImages([
    pie.Image("user/pieadmin/saas/06/s2_huaian_mar"),
    pie.Image("user/pieadmin/saas/06/s2_huaian_apr"),
    pie.Image("user/pieadmin/saas/06/s2_huaian_may"),
    pie.Image("user/pieadmin/saas/06/s2_huaian_june"),
    pie.Image("user/pieadmin/saas/06/s2_huaian_july"),
    pie.Image("user/pieadmin/saas/06/s2_huaian_aug"),
    pie.Image("user/pieadmin/saas/06/s2_huaian_sept"),
    pie.Image("user/pieadmin/saas/06/s2_huaian_oct")
]).select(['B1', 'B2'])
var mar = col2.getAt(0).addBands(pie.Image().constant(74).rename('doy'))
var apr2 = col2.getAt(1).addBands(pie.Image().constant(105).rename('doy'))
var may2 = col2.getAt(2).addBands(pie.Image().constant(135).rename('doy'))
var june2 = col2.getAt(3).addBands(pie.Image().constant(166).rename('doy'))
var july2 = col2.getAt(4).addBands(pie.Image().constant(196).rename('doy'))
var aug = col2.getAt(5).addBands(pie.Image().constant(227).rename('doy'))
var sept = col2.getAt(6).addBands(pie.Image().constant(258).rename('doy'))
var oct = col2.getAt(7).addBands(pie.Image().constant(288).rename('doy'))
var growth = pie.ImageCollection.fromImages([mar, apr2, may2, june2, july2, aug, sept, oct]) //B1,B2,doy   

阈值筛选/
//插秧期合成negVI最大值，并记录doy1
var flood = transplant.qualityMosaic('negVI').select(['B2', 'negVI', 'doy']).clip(roi)
//生长期合成EVI最大值，并记录doy2
var peak = growth.qualityMosaic('B2').select(['B2', 'doy']).clip(roi) //B2 doy
//计算水稻插秧期快速生长后EVI值的变化
var datediffer = (peak.select('doy').subtract(flood.select('doy'))).abs()
var ratio = (peak.select('B2').subtract(flood.select('B2'))).divide(datediffer.divide(30)).rename('ratio') //.clamp(-1,1)//(0,1)

//确定阈值并筛选
var floodmask = flood.select('negVI').gte(-0.26) //flood:LSWI-EVI>-0.26排除其他植被、城镇像素
var peakmask = peak.select('B2').gte(0.47) //peak:EVI>0.47排除水体、城镇像素
var ratiomask = ratio.select('ratio').gte(0.13) //(peak_EVI - flood_EVI)/(date/30)>0.13 排除其他植被中EVI变化速率小但均值较高的像素
var thersold = symmetricDifference.updateMask(floodmask).updateMask(peakmask).updateMask(ratiomask)
//加载显示
Map.addLayer(thersold, visParam, '阈值筛选后结果')
Map.addLayer(symmetricDifference, { min: 0, max: 1, palette: '00FFFF, FFE100' }, '并集')

//导出阈值筛选后图像
Export.image({
    image: thersold,
    description: "result",
    assetId: "result_huaian",
    region: roi,
    scale: 30
});