PIE-Engine 教程:水稻面积提取1(宿迁市)

我们这次用国家的边界通过name进行筛选到宿迁市进行分析，本次计算的水稻面积用到的指数是NDVI，LSWI，EVI三种指数进行分析，用到的影像是Landsat 8数据，DEM数据（计算坡度、坡向）还有指定5月和8月份SAR影像（"VH1","VH2"）最后分别加入到Landsat波段影像，可导出我们所合成的影像。先看所用到的函数：

aspect(input)

计算地形DEM数据的坡向信息（以度为单位），局部梯度是用每个像素的4个连通邻域来计算的，因此图像的边缘会出现缺失值。

方法参数：

- terrain(Terrain)

DEM数据处理方法。

- input(Image)

地形DEM数据，以米为单位。

返回值：Image

slope(image)

计算地形DEM数据的坡度信息（以度为单位），局部梯度是用每个像素的4个连通邻域来计算的，因此图像的边缘会出现缺失值。

方法参数：

- terrain(Terrain)

DEM数据处理方法。

- image(Image)

地形DEM数据，以米为单位。

返回值：Image

addBands(srcImg,names,overwrite)

添加影像中的波段到当前影像。

方法参数：

- dstImg(Image)

Image实例。

- srcImg(Image)

要添加波段的影像对象。

- names(String|Array<String>)

波段名称列表，默认为null，表示全部添加。

- overwrite(Boolean)

相同名称是否覆盖，默认是false。

返回值：Image

map(algorithm,dropNulls)

针对ImageCollection中的每个Image进行循环计算。

方法参数：

- imageCollection(ImageCollection)

ImageCollection实例。

- algorithm(Function)

参数是Image的方法，需要返回Image

- dropNulls(Boolean)

是否允许返回空值

返回值：ImageCollection

代码：

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

//1.合成监督分类影像
var area = pie.Geometry.Polygon([
    [
        [
            118.19042604840365,
            32.717477656827796
        ],
        [
            119.67193792508522,
            32.717477656827796
        ],
        [
            119.67193792508522,
            34.12687250105918
        ],
        [
            118.19042604840365,
            34.12687250105918
        ],
        [
            118.19042604840365,
            32.717477656827796
        ]
    ]
], null);

//加载淮安市矢量
var roi = pie.FeatureCollection("NGCC/CHINA_CITY_BOUNDARY")
    .filter(pie.Filter.eq('name', '淮安市'))
    .first()
    .geometry()
Map.centerObject(roi, 8)
Map.addLayer(roi, { color: "red", fillColor: "00000000", width: 2 }, "Huaian")

//NDVI
function NDVI(image) {
    var red = image.select("B4")
    var nir = image.select("B5")
    var ndvi = (nir.subtract(red)).divide(nir.add(red)).multiply(10000)
    return image.addBands(ndvi.rename("NDVI"));
}

//LSWI
function LSWI(image) {
    var swir = image.select("B6")
    var nir = image.select("B5")
    var lswi = (nir.subtract(swir)).divide(nir.add(swir)).multiply(10000)
    return image.addBands(lswi.rename("LSWI"));
}

//EVI
function EVI(image) {
    var red = image.select("B4").divide(10000)
    var nir = image.select("B5").divide(10000)
    var blue = image.select("B2").divide(10000)
    var evi = ((nir.subtract(red)).multiply(2.5))
        .divide(nir.add(red.multiply(6)).subtract(blue.multiply(7.5)).add(1))
        .multiply(10000)
    return image.addBands(evi.rename("EVI"));
}

//去云处理
function cloudMask(image) {
    var qa = image.select('BQA');
    var cloudMask = qa.bitwiseAnd(1 << 4).eq(0);
    return image.updateMask(cloudMask);
};

//选择用于分类的影像并裁剪
var l8col = pie.ImageCollection("LC08/01/T1")
    .filterBounds(roi)
    .filterDate("2020-3-1", "2020-12-1")
    .select(["B2", "B3", "B4", "B5", "B6", "B7", "BQA"])
    .map(cloudMask)
    .map(NDVI)
    .map(LSWI)
    .map(EVI)
    .median()

//选取dem
var dem = pie.ImageCollection('DEM/SRTM_30')
    .filterBounds(roi)
    .select('elevation')
    .mosaic()
    .clip(roi);
var aspect = pie.Terrain.aspect(dem);
var slope = pie.Terrain.slope(dem).multiply(180 / 3.1415926);

//选取SAR 5月8月影像
var sarmay = pie.Image('user/pieadmin/saas/06/sar_may').select('B1').multiply(10000).clip(roi);
var saraug = pie.Image('user/pieadmin/saas/06/sar_aug').select('B1').multiply(10000).clip(roi);

//选择用于分类的波段         
var bands = [
    "B2", "B3", "B4", "B5", "B6", "B7", "NDVI", "LSWI", "EVI"
]; 

//"VH1","VH2","ASPECT","SLOPE"
//合成用于监督分类的影像
var l8Image = l8col
    .addBands(aspect.rename("ASPECT"))
    .addBands(slope.rename("SLOPE"))
    .addBands(sarmay.rename("VH1"))
    .addBands(saraug.rename("VH2"))
    .select(bands)
    .clip(roi)

//导出影像
Export.image({
    image: l8Image,
    description: "l8Image",
    assetId: "l8Image",
    region: area,
    scale: 30
});