#冲刺创作新星#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
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.
32.
33.
34.
35.
36.
37.
38.
39.
40.
41.
42.
43.
44.
45.
46.
47.
48.
49.
50.
51.
52.
53.
54.
55.
56.
57.
58.
59.
60.
61.
62.
63.
64.
65.
66.
67.
68.
69.
70.
71.
72.
73.
74.
75.
76.
77.
78.
79.
80.
81.
82.
83.
84.
85.
86.
87.
88.
89.
90.
91.
92.
93.
94.
95.
96.
97.
98.
99.
100.
101.
102.
103.
104.
105.
106.
107.
108.
109.
110.
111.
112.
113.
114.
115.
116.
117.
118.
119.
120.
121.
122.
123.
124.
125.
126.
127.
128.
129.
130.
131.
132.
133.
134.
135.
136.
137.
138.
139.
140.
141.
142.
143.
144.
145.
146.
147.
148.
149.
150.
151.
152.
153.
154.
155.
156.
157.
158.
159.
160.
161.
162.
163.
164.
165.
166.
167.
168.
169.
170.
171.
172.
173.
174.
175.
176.
177.
178.
179.
180.
181.
182.
183.
184.
185.
186.
187.
188.
189.
190.
191.
192.
193.
194.
195.
196.
197.
198.
199.
200.
201.
202.
203.
204.
205.
206.
207.
208.
209.
210.
211.
212.
213.
214.
215.
216.
217.
218.
219.
220.
221.
222.
223.
224.
225.
226.
227.
228.
229.
230.
231.
232.
233.
234.
235.
236.
237.
238.