#冲刺创作新星#PIE-Engine：广东省生态遥感指数研究 原创

/**
 * @Name    :   广东省遥感生态指数研究
 * @Author  :   广州大学张三的组
 * @Source  :   航天宏图第四届 “航天宏图杯”PIE软件二次开发大赛云开发组三等奖获奖作品
 * @Description

var featureCollection0 = pie.FeatureCollection("NGCC/CHINA_PROVINCE_BOUNDARY");
var gd = featureCollection0.filter(pie.Filter.eq("name", "广东省")).first().geometry();

//在地图上显示广东省区域矢量图
Map.addLayer(gd, { color: "ffff00ff", fillColor: "00000000" }, "广东省行政区划")
Map.centerObject(gd, 5)

//导入所需数据集 
var MOD09A1 = pie.ImageCollection("USGS/MOD09A1/006");
var MOD11A2 = pie.ImageCollection("USGS/MOD11A2/006");
var MOD13A1 = pie.ImageCollection("USGS/MOD13A1/006");

//构建云掩膜函数
function cloudfree_mod09a1(image) {
    var qa = image.select('sur_refl_state_500m')
    var cloudState0 = 1 << 0;
    var cloudState1 = 1 << 1;
    var cloudShadowState = 1 << 2;
    var cirrusState8 = 1 << 8;
    var cirrusState9 = 1 << 9;
    var mask = qa.bitwiseAnd(cloudState0).eq(0)
        .and(qa.bitwiseAnd(cloudState1).eq(0))
        .and(qa.bitwiseAnd(cloudShadowState).eq(0)) // No cloud shadow
        .and(qa.bitwiseAnd(cirrusState8).eq(0))
        .and(qa.bitwiseAnd(cirrusState9).eq(0)) // No cirrus
    return image.updateMask(mask)
}

//筛选09A1影像准备真色彩波段验证云掩膜
var recentimage = MOD09A1
    .filterDate("2020-11-01", "2020-12-30")
    .filterBounds(gd)
    .select(["sur_refl_b01", "sur_refl_b03", "sur_refl_b04", "sur_refl_b06", "sur_refl_state_500m"])
    .map(cloudfree_mod09a1)
    .mosaic()
    .clip(gd);
var visParam = {
    min: -100,
    max: 3000,
    bands: ["sur_refl_b01", "sur_refl_b04", "sur_refl_b03"]
}

//构建水体指数函数
function MNDWI(img) {
    var blue = img.select("sur_refl_b03")
    var green = img.select("sur_refl_b04");
    var swir1 = img.select("sur_refl_b06");
    var swir2 = img.select("sur_refl_b07");
    var nir = img.select("sur_refl_b02");
    var mndwi = green.subtract(swir1).divide(green.add(swir1)).rename(mndwi)
    var aweinsh = blue.add(green.multiply(2.5))
        .subtract(nir.add(swir1).multiply(1.5))
        .subtract(swir2.multiply(0.25))
        .rename(aweinsh);
    var awei = green.subtract(swir1)
        .multiply(4)
        .subtract(nir.multiply(0.25).add(swir2.multiply(2.75)))
        .rename(awei);
    var ndwi = green.subtract(nir)
        .divide(green.add(nir))
        .rename(ndwi);
    return mndwi
}

//以阈值为0创建水体掩膜
var WaterMask = MNDWI(recentimage).lt(0);
//对影像进行水体掩膜运算
var wa = recentimage.updateMask(WaterMask);
//影像可视化渲染
Map.addLayer(wa, visParam, "ss")

//筛选时间地区，水体掩膜处理
var LST = MOD11A2.filterDate("2020-12-01", "2020-12-30")
    .filterBounds(gd)
    .select('LST_Day_1km')
    .mean()
    .clip(gd)
    .multiply(0.02)
    .subtract(273.15)
    .rename('lst')
    .updateMask(WaterMask)
//计算热度指数最大值
var maxLst = pie.Number(LST.reduceRegion(pie.Reducer.max(), gd, 1000).get('lst')) 
//计算热度指数平均值
var aveLst = pie.Number(LST.reduceRegion(pie.Reducer.mean(), gd, 1000).get('lst'))
//计算热度指数最小值 
var minLst = pie.Number(LST.reduceRegion(pie.Reducer.min(), gd, 1000).get('lst'))
//通过最大最小值计算归一化热度指数平均值 
var noraveLst = aveLst.subtract(minLst).divide(maxLst.subtract(minLst)) 

var LSTVis = {
    min: 4.6,
    max: 27.5,
    palette: [
        '040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6',
        '0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef',
        '3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f',
        'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d',
        'ff0000', 'de0101', 'c21301', 'a71001', '911003'
    ]
};
//热度指数可视化渲染
Map.addLayer(LST, LSTVis, 'LST'); 

var NDVI = MOD13A1.filterDate("2020-12-01", "2020-12-30")
    .filterBounds(gd)
    .select('NDVI')
    .mean()
    .clip(gd)
    .multiply(0.0001)
    .rename('ndvi');
//计算绿度指数最大值    
var maxNDVI = pie.Number(NDVI.reduceRegion(pie.Reducer.max(), gd, 500).get('ndvi')) 
//计算绿度指数平均值
var aveNDVI = pie.Number(NDVI.reduceRegion(pie.Reducer.mean(), gd, 500).get('ndvi'))
//计算绿度指数最小值 
var minNDVI = pie.Number(NDVI.reduceRegion(pie.Reducer.min(), gd, 500).get('ndvi'))
//通过最大最小值计算归一化植被指数平均值 
var noraveNDVI = aveNDVI.subtract(minNDVI).divide(maxNDVI.subtract(minNDVI)) 

var ndviVis = {
    min: -0.19,
    max: 0.9,
    palette: [
        'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
        '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
        '012E01', '011D01', '011301'
    ],
};
//绿度指数可视化渲染
Map.addLayer(NDVI, ndviVis, "ndvi") 

//筛选影像，去云处理，水体掩膜处理

var srIMG0 = MOD09A1
    .filterDate("2020-11-01", "2020-12-30")
    .filterBounds(gd)
    .select(["sur_refl_b01", "sur_refl_b02", "sur_refl_b03", "sur_refl_b04",
        "sur_refl_b05", "sur_refl_b06", "sur_refl_b07", "sur_refl_state_500m"
    ])
    .map(cloudfree_mod09a1)
    .mosaic()
    .clip(gd)
    .updateMask(WaterMask);

var WET = srIMG0.select("sur_refl_b01").multiply(0.1147) //红波段
    .add(srIMG0.select("sur_refl_b02").multiply(0.2489)) //近红外波段1
    .add(srIMG0.select("sur_refl_b03").multiply(0.2408)) //蓝波段
    .add(srIMG0.select("sur_refl_b04").multiply(0.3132)) //绿波段
    .add(srIMG0.select("sur_refl_b05").multiply(-0.3122)) //近红外波段2
    .add(srIMG0.select("sur_refl_b06").multiply(-0.6416)) //短波红外1
    .add(srIMG0.select("sur_refl_b07").multiply(-0.5087)) //短波红外2
    .multiply(0.0001)
    .rename("wet")
//计算湿度指数
var maxWET = pie.Number(WET.reduceRegion(pie.Reducer.max(), gd, 500).get('wet')) //计算湿度指数最大值
var aveWET = pie.Number(WET.reduceRegion(pie.Reducer.mean(), gd, 500).get('wet')) //计算湿度指数平均值
var minWET = pie.Number(WET.reduceRegion(pie.Reducer.min(), gd, 500).get('wet')) //计算湿度指数最小值
var noraveWET = aveWET.subtract(minWET).divide(maxWET.subtract(minWET)) //通过最大最小值计算归一化湿度指数平均值

var wetVis = {
    min: -0.3,
    max: 0.01,
    palette: [
        '040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6',
        '0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef',
        '3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f',
        'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d',
        'ff0000', 'de0101', 'c21301', 'a71001', '911003'
    ]
}
Map.addLayer(WET, wetVis, "wet") //湿度指数可视化渲染

var sswir1 = srIMG0.select("sur_refl_b06");
var rred = srIMG0.select("sur_refl_b01");
var nnir1 = srIMG0.select("sur_refl_b02");
var bblue = srIMG0.select("sur_refl_b03");
var ggreen = srIMG0.select("sur_refl_b04");
var si = sswir1.add(rred)
        .subtract(nnir1.add(bblue))
        .divide(sswir1.add(rred).add(nnir1.add(bblue)));
//计算干度指数分量裸土指数
var ibi = sswir1.multiply(2).divide(sswir1.add(nnir1))
    .subtract(nnir1.divide(nnir1.add(rred)).add(ggreen.divide(ggreen.add(sswir1))))
    .divide(sswir1.multiply(2).divide(sswir1.add(nnir1))
        .add(nnir1.divide(nnir1.add(rred)).add(ggreen.divide(ggreen.add(sswir1)))));
//计算干度指数分量城市建筑指数
var NDBSI = si.add(ibi).divide(2).rename('ndbsi');
var maxNDBSI = pie.Number(NDBSI.reduceRegion(pie.Reducer.max(), gd, 500).get('ndbsi')) //计算干度指数最大值
var aveNDBSI = pie.Number(NDBSI.reduceRegion(pie.Reducer.mean(), gd, 500).get('ndbsi')) //计算干度指数平均值
var minNDBSI = pie.Number(NDBSI.reduceRegion(pie.Reducer.min(), gd, 500).get('ndbsi')) //计算干度指数最小值
var noraveNDBSI = aveNDBSI.subtract(minNDBSI).divide(maxNDBSI.subtract(minNDBSI)) //通过最大最小值计算归一化干度指数平均值

var ndbsiVis = {
    min: -0.3,
    max: 0.01,
    palette: [
        '040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6',
        '0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef',
        '3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f',
        'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d',
        'ff0000', 'de0101', 'c21301', 'a71001', '911003'
    ]
}

Map.addLayer(NDBSI, ndbsiVis, "ndbsi") //干度指数可视化渲染

var index = LST.addBands(NDVI).addBands(WET).addBands(NDBSI)

var bar_options = {
    title: "广东省RSEI分量归一化平均值",
    legend: ["2020年12月"],
    yAxis: ["LST", "NDVI", "WET", "NDBSI"],
    xAxisName: "指数值",
    yAxisName: "指数",
    series: [
        [0.577, 0.729, 0.728, 0.321]
    ],
    chartType: "bar",
};

//调用绘制方法，输出图表显示在结果面板中
var chart = ui.Chart.array(bar_options);
chart.setStyle({ width: "380px" });
print(chart);

Export.image({
    image: index,
    description: "index",
    assetId: "result",
    region: gd,
    scale: 500