-
Notifications
You must be signed in to change notification settings - Fork 0
/
Flavours_of_Unmixing_with_rmse_C2.js
298 lines (223 loc) · 11.5 KB
/
Flavours_of_Unmixing_with_rmse_C2.js
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
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
/*
#+
# :AUTHOR: Katarzyna Ewa Lewinska [[email protected]]
# :URL: https://github.com/erfea/Google_Earth_Engine/blob/master/Flavours_of_Unmixing_with_rmse.js
# :DATE: 5 February 2019
#
# :Description: Script allows on comparison between spectral unmixing executed with unmix and matrixSolve function.
# Due to limitations the latter approach shade endmember is introduced as [1,1,1,1,1,1] (an endmember cannot be [0,0,0,0,0,0]).
# The code comes with the solveMartix solution being commented out.
# In the second part of a code, a comparison among results of the unmix approach run for 4 combinations of
# sumToOne and nonNegative constrains is shown. In addition to endmembers' abundances, rmse for each solution is calculated
# through a 'reversed' approach. (A reversed image is constructed based on endmembers and their obtained abundances.
# Band-wise differences between the original image and the reversed image are used to calculate rmse for each model).
# The script allows on visualization on all results in a quadruple display, where each considered model
# is presented in a different map window.
# In the current example the code runs for a single Landsat 8 OLI data, converting it to the ETM+ space. The code can be
# adjusted to work for TM and ETM+ scenes, as well as long time series of Landsat data. Implemented cloud masking is optional.
# A scene to be used is identified as a first image in the defined search window overlapping with a used-defined geometry point.
# Endmembers definitions are based on [10.1016/j.rse.2005.07.013].
# Scaling factor for endmembers 0.0001. Endmembers are defined in C1 range
#
# :Input: point - geometry point to identify a corresponding Landsat tile
# start - a string YYYY-MM-DD defining the first day of the scene search window
# end - a string YYYY-MM-DD defining the last day of the scene search window
# soil, gv, npv, shade - endmembers introduced as separate lists
#
# :Output: Dynamic quadruple linked display. Values can be inspected in the most-left map view.
#
# :Updates: 2019-02-05: quadruple linked display; RMSE
# 2019-04-29: map0.style().set('cursor', 'crosshair')
# 2021-10-24: update to Collection 2
#
# :2Do:
#
# :Disclaimer: The author of this code accepts no responsibility for errors or omissions in this work
# and shall not be liable for any damage caused by these.
#-
*/
// ### INPUTS ### \\
var point = ee.Geometry.Point([49.0214, 40.4356]);
var start = '2019-06-10';
var end = '2019-06-20';
var endmembers = ee.Dictionary({
'soil': [1670, 2600, 3400, 3900, 4430, 3750],
'gv': [500, 900, 400, 6100, 3000, 1000],
'npv': [1080, 1326, 1770, 2480, 4150, 3346],
'shade': [1, 1, 1, 1, 1, 1],
});
var vis = {"opacity":1,"bands":["RMSE"],"max":1500.0,"gamma":1};
// ### ENVIRONMENT ### \\
// ETM+ cloud mask
var cloudMaskL8 = function(in_image) {
var qa = in_image.select('QA_PIXEL');
var mask = qa.bitwiseAnd(1 << 0)
.or(qa.bitwiseAnd(1 << 1)) // dialated cloud
.or(qa.bitwiseAnd(1 << 2)) // cirrus
.or(qa.bitwiseAnd(1 << 14).and(qa.bitwiseAnd(1 << 15))) // cirrus confidence
.or(qa.bitwiseAnd(1 << 3)) //cloud
.or(qa.bitwiseAnd(1 << 8).and(qa.bitwiseAnd(1 << 9))) // cloud confidence
.or(qa.bitwiseAnd(1 << 4)) //shadow
.or(qa.bitwiseAnd(1 << 10).and(qa.bitwiseAnd(1 << 11))) //cloud shadow confidence
.or(qa.bitwiseAnd(1 << 5))// snow
.or(qa.bitwiseAnd(1 << 12).and(qa.bitwiseAnd(1 << 13))) // Snow confidence
return in_image.updateMask(mask.not())//.updateMask(mask2);
};
// OLI to ETM function
// after doi:10.1016/j.rse.2015.12.024
var OLI2ETMf = function(image) {
var blue = image.expression('0.0183*10000 + 0.8850 * (SR_B2*0.0000275-0.2)*10000', {'SR_B2':image.select('SR_B2')});
var green = image.expression('0.0123*10000 + 0.9317 * (SR_B3*0.0000275-0.2)*10000', {'SR_B3':image.select('SR_B3')});
var red = image.expression('0.0123*10000 + 0.9372 * (SR_B4*0.0000275-0.2)*10000', {'SR_B4':image.select('SR_B4')});
var nir = image.expression('0.0448*10000 + 0.8339 * (SR_B5*0.0000275-0.2)*10000', {'SR_B5':image.select('SR_B5')});
var swir1 = image.expression('0.0306*10000 + 0.8639 * (SR_B6*0.0000275-0.2)*10000', {'SR_B6':image.select('SR_B6')});
var swir2 = image.expression('0.0116*10000 + 0.9165 * (SR_B7*0.0000275-0.2)*10000', {'SR_B7':image.select('SR_B7')});
// *0.0000275 -0.2
var temp = ee.Image(image).addBands(blue).addBands(green).addBands(red).addBands(nir).addBands(swir1).addBands(swir2);
return temp.select(['constant', 'constant_1', 'constant_2', 'constant_3', 'constant_4', 'constant_5' ],
['B1','B2','B3','B4','B5','B7']).int16();
};
// ### CODE ### \\
var l8 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2');
var ts_l8 = ee.ImageCollection(l8.filterDate(start, end).filterBounds(point));
// ts_l8 = ts_l8.map(cloudMaskL8); // cloud mask
ts_l8 = ts_l8.map(OLI2ETMf); // OLI -> ETM+ normalization
var image = ts_l8.first();
print(image);
var bands = image.bandNames();
// // unmixing through matrixSolve()
// var end_arr = ee.Array.cat([soil, gv, npv, shade],1);
// print('end_arr', end_arr, end_arr.length()); //[6x5]
// // Make an Array Image with a 2-D Array per pixel, 6x1.
// var arrayImage = image.toArray().toArray(1); //[6x1]
// print('arrayImage', arrayImage);
// // [4x1] [6x4] [6x1] -> [6x4]x[4x1]=[6x1]
// var unmixed = ee.Image(end_arr).matrixSolve(arrayImage);
// // var unmixedImage = unmixed.arrayProject([0]).arrayFlatten([['soil', 'gv', 'npv', 'shade']]); // 0 image axis
// print('unmixedImage', unmixedImage);
// Map.addLayer(unmixedImage, {}, 'matrixSolve');
var u_soil = endmembers.get('soil');
var u_gv = endmembers.get('gv');
var u_npv =endmembers.get('npv');
var u_shade = endmembers.get('shade');
// sumToOne = F and nonNegative=F
var unmixFF = image.unmix([u_soil, u_gv, u_npv, u_shade], false, false).rename(['soil','gv','npv','shade']);
// sumToOne = T and nonNegative=T
var unmixTT = image.unmix([u_soil, u_gv, u_npv, u_shade], true, true).rename(['soil','gv','npv','shade']);
// sumToOne = T and nonNegative=false
var unmixTF = image.unmix([u_soil, u_gv, u_npv, u_shade], true, false).rename(['soil','gv','npv','shade']);
// sumToOne = F and nonNegative=T
var unmixFT = image.unmix([u_soil, u_gv, u_npv, u_shade], false, true).rename(['soil','gv','npv','shade']);
// ### ERROR TERMS ### \\
// var end_arr = ee.Array.cat([soil, gv, npv, shade],1);
var end_arr = ee.Array.cat([u_soil, u_gv, u_npv, u_shade],1);
// reverse the image from spectra and endmembers
var unmixArFF = unmixFF.toArray().toArray(1); //[5x1]
var revImFF = ee.Image(end_arr).matrixMultiply(unmixArFF); // [6x5] x [5x1] -> [6x1]
var revImFF = revImFF.arrayProject([0]).arrayFlatten([bands]);
var unmixArTT = unmixTT.toArray().toArray(1); //[5x1]
var revImTT = ee.Image(end_arr).matrixMultiply(unmixArTT); // [6x5] x [5x1] -> [6x1]
var revImTT = revImTT.arrayProject([0]).arrayFlatten([bands]);
var unmixArTF = unmixTF.toArray().toArray(1); //[5x1]
var revImTF = ee.Image(end_arr).matrixMultiply(unmixArTF); // [6x5] x [5x1] -> [6x1]
var revImTF = revImTF.arrayProject([0]).arrayFlatten([bands]);
var unmixArFT = unmixFT.toArray().toArray(1); //[5x1]
var revImFT = ee.Image(end_arr).matrixMultiply(unmixArFT); // [6x5] x [5x1] -> [6x1]
var revImFT = revImFT.arrayProject([0]).arrayFlatten([bands]);
var difFF = image.subtract(revImFF);
var difTT = image.subtract(revImTT);
var difTF = image.subtract(revImTF);
var difFT = image.subtract(revImFT);
var difFFsq = ((difFF.select('B1').pow(ee.Image(2)))
.add(difFF.select('B2').pow(ee.Image(2)))
.add(difFF.select('B3').pow(ee.Image(2)))
.add(difFF.select('B4').pow(ee.Image(2)))
.add(difFF.select('B5').pow(ee.Image(2)))
.add(difFF.select('B7').pow(ee.Image(2))).divide(ee.Image(6))).sqrt().rename('RMSE');
var difTTsq = ((difTT.select('B1').pow(ee.Image(2)))
.add(difTT.select('B2').pow(ee.Image(2)))
.add(difTT.select('B3').pow(ee.Image(2)))
.add(difTT.select('B4').pow(ee.Image(2)))
.add(difTT.select('B5').pow(ee.Image(2)))
.add(difTT.select('B7').pow(ee.Image(2))).divide(ee.Image(6))).sqrt().rename('RMSE');
var difTFsq = ((difTF.select('B1').pow(ee.Image(2)))
.add(difTF.select('B2').pow(ee.Image(2)))
.add(difTF.select('B3').pow(ee.Image(2)))
.add(difTF.select('B4').pow(ee.Image(2)))
.add(difTF.select('B5').pow(ee.Image(2)))
.add(difTF.select('B7').pow(ee.Image(2))).divide(ee.Image(6))).sqrt().rename('RMSE');
var difFTsq = ((difFT.select('B1').pow(ee.Image(2)))
.add(difFT.select('B2').pow(ee.Image(2)))
.add(difFT.select('B3').pow(ee.Image(2)))
.add(difFT.select('B4').pow(ee.Image(2)))
.add(difFT.select('B5').pow(ee.Image(2)))
.add(difFT.select('B7').pow(ee.Image(2))).divide(ee.Image(6))).sqrt().rename('RMSE');
/// ### MAPS ### \\\
var long1= point.getInfo()['coordinates'][0];
var lat1=point.getInfo()['coordinates'][1];
function initMap(map) {
map.setCenter(long1,lat1, 9);
}
// Initialize Map
initMap(Map);
function createMap(title) {
var map = ui.Map();
ui.Label(title, {position:'bottom-center'});
map.add(title);
return map;
}
function getMapSize() {
var scale = Map.getScale();
var bounds = ee.Geometry(Map.getBounds(true)).coordinates().get(0).getInfo();
var ll = bounds[0];
var ur = bounds[2];
var width = (ur[0] - ll[0]) / scale;
var height = (ur[1] - ll[1]) / scale;
return { w: Math.floor(width), h: Math.floor(height) };
}
var height = getMapSize().h;
var map0 = ui.Map();
var map1 = ui.Map();
var map2 = ui.Map();
var map3 = ui.Map();
map0.addLayer(difFTsq, vis, 'FT Error', false);
map0.addLayer(difTFsq, vis, 'TF Error', false);
map0.addLayer(difTTsq, vis, 'TT Error', false);
map0.addLayer(difFFsq, vis, 'FF Error', false);
map0.addLayer(unmixFT, {}, 'FT', false);
map0.addLayer(unmixTF, {}, 'TF', false);
map0.addLayer(unmixTT, {}, 'TT', false);
map0.addLayer(unmixFF, {}, 'FF');
// map0.layers().set(10, dot);
map1.addLayer(difTTsq, vis, 'TT Error');
map1.addLayer(unmixTT, {}, 'TT');
map2.addLayer(difTFsq, vis, 'TF Error');
map2.addLayer(unmixTF, {}, 'TF');
map3.addLayer(difFTsq, vis, 'FT Error');
map3.addLayer(unmixFT, {}, 'FT');
map0.add(ui.Label('SumToOne= F nonNegative= F', {position:'bottom-center'}));
map1.add(ui.Label('SumToOne= T nonNegative= T', {position:'bottom-center'}));
map2.add(ui.Label('SumToOne= T nonNegative= F', {position:'bottom-center'}));
map3.add(ui.Label('SumToOne= F nonNegative= T', {position:'bottom-center'}));
var maps = [map0, map1, map2, map3];
// Create a panel with vertical flow layout.
var panel = ui.Panel({
layout: ui.Panel.Layout.flow('horizontal'),
style: {width: '100vw', height: height + '300px'}
});
var linker = ui.Map.Linker(maps);
maps.map(function(map) {
initMap(map);
// map.setOptions('HYBRID');
panel.add(map);
});
ui.root.clear();
ui.root.add(maps[0]);
ui.root.add(maps[1]);
ui.root.add(maps[2]);
ui.root.add(maps[3]);
var linker = ui.Map.Linker(maps);
map0.style().set('cursor', 'crosshair');
/*
*/
// this is the end of the code \\