-
Notifications
You must be signed in to change notification settings - Fork 6
/
Figure6a.js
401 lines (329 loc) · 17.8 KB
/
Figure6a.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
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
/*
MIT License
Copyright (c) 2021 Mong-Han Huang and Alexander L. Handwerger
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/
// #############################################################################
// ### Event info ###
// #############################################################################
/*
/*
* SAR-based landslide (and other ground surface) detection
* Study Site: Hokkaido, Japan
* Event: Mw 6.7 earthquake on 6 September 2018
* Code written by Mong-Han Huang, Department of Geology, [email protected]
* and Alexander L. Handwerger, UCLA JIFRESSE/NASA JPL, [email protected] and [email protected] with help from many users on stackexchange!
*/
// #############################################################################
// ### Set AOI ###
// #############################################################################
var AOI = AOI; // define Area of Interest
Map.setOptions("TERRAIN");
Map.centerObject(AOI, 12); //zooms to center of AOI after clicking "run". The number determines the zoom level.
// #############################################################################
// ### Define pre-event and post-event time periods ###
// #############################################################################
//define pre-event stack time period
var PreEventTime_1 = '2015-08-01T23:59'; // format: yyyy-mm-dd-HH:MM
var PreEventTime_2 = '2018-09-05T23:59'; // format: yyyy-mm-dd-HH:MM
//define post-event stack time period for event inventory detection
//for this approach use all of the available post-event imagery
var PostEventTime_1 = '2018-09-06T23:59'; // format: yyyy-mm-dd-HH:MM
var PostEventTime_2 = '2018-09-21T23:59'; // format: yyyy-mm-dd-HH:MM
var PostEventTime_S2longer = '2019-09-25T23:59'; // format: yyyy-mm-dd-HH:MM
// #############################################################################
// ### slope and curvature thresholds ###
// #############################################################################
// DEM thresholds
var slope_threshold = 10; // Exclude areas with hillslope angle < slope_threshold, unit: degree
var curv_threshold = -0.005; // Exclude areas with hillslope curvature > curv_threshold, unit: m/m^2
// #############################################################################
// ### NASA DEM ###
// #############################################################################
// load NASA Digital Elevation Model (DEM)
var NASADEM_dataset = ee.Image('NASA/NASADEM_HGT/001')
var elevation = NASADEM_dataset.select('elevation');
elevation=elevation.clip(AOI); //clip to AOI
var waterZones = NASADEM_dataset.select('swb'); //load in water body data
var waterMask = waterZones.eq(0) // Create a binary water mask.
var exaggeration = 1; //changes the appearance of the hillshade
var hillshade = ee.Terrain.hillshade(elevation.multiply(exaggeration));
hillshade=hillshade.clip(AOI);
var slope = ee.Terrain.slope(elevation); // slope angle in degrees
slope=slope.clip(AOI);
var mask_slope = slope.gte(slope_threshold); // slope mask with values 0 or 1, removes values less than or equal to threshold
var slope_masked = slope.updateMask(mask_slope); // slope angle with values < threshold excluded
// Calculate hillslope curvature
// Define a Gaussian kernel for smoothing. This step helps reduce noise in the curvature maps
var smooth_curv = ee.Kernel.gaussian({
radius: 60,
sigma: 30,
units: 'meters',
normalize: true,
});
// Smoothing the DEM with the gaussian kernel.
var elevation_smooth= elevation.convolve(smooth_curv).resample("bilinear");
var xyDemGrad = elevation_smooth.gradient().resample("bilinear");
var xGradient = xyDemGrad.select('x').gradient().resample("bilinear");
var yGradient = xyDemGrad.select('y').gradient().resample("bilinear");
var curvature = xGradient.select('x').add(yGradient.select('y'));
var mask_curvature = curvature.gte(curv_threshold);
var curvature_masked = curvature.updateMask(mask_curvature);
//var DEMMask= mask_slope.multiply(mask_curvature).multiply(waterMask); //slope, curvature, and water mask
// #############################################################################
// ### Sentinel-2 Optical Images ###
// #############################################################################
// cloud filter and mask for Sentinel-2 optical data
function maskS2clouds(image) {
var qa = image.select('QA60');
// Bits 10 and 11 are clouds and cirrus, respectively.
var cloudBitMask = 1 << 10;
var cirrusBitMask = 1 << 11;
// Both flags should be set to zero, indicating clear conditions.
var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
.and(qa.bitwiseAnd(cirrusBitMask).eq(0));
return image.updateMask(mask).divide(10000);
}
// Load Sentinel-2 reflectance data.
var S2_PreEvent = ee.ImageCollection('COPERNICUS/S2')
.filterDate(PreEventTime_1,PreEventTime_2)
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 10)) //only include images with less than 10% clouds
.filterBounds(AOI)
.map(maskS2clouds);
var S2_PreEvent_median=S2_PreEvent.median();
S2_PreEvent_median=S2_PreEvent_median.clip(AOI);
var S2_PostEvent = ee.ImageCollection('COPERNICUS/S2')
.filterDate(PostEventTime_1, PostEventTime_S2longer)
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 10)) //only include images with less than 10%
.filterBounds(AOI)
.map(maskS2clouds);
var S2_PostEvent_median = S2_PostEvent.median();
S2_PostEvent_median=S2_PostEvent_median.clip(AOI);
// #############################################################################
// ### Sentinel-1 SAR Images ###
// #############################################################################
// LOAD Sentinel-1 (S1) Backscatter Coefficient data VH polarization
var imgVH = ee.ImageCollection('COPERNICUS/S1_GRD')
.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
.filter(ee.Filter.eq('instrumentMode', 'IW'))
.select('VH') // VH polarization only
.filterBounds(AOI)
.map(function(image) {
var edge = image.lt(-30.0); //remove low edge values as suggested by GEE
var maskedImage = image.mask().and(edge.not());
return image.updateMask(maskedImage);
});
var desc = imgVH.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING')); //descending acquisition geometry data
var asc = imgVH.filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING')); //ascending acquisition geometry data
var PreEventPeriod = ee.Filter.date(PreEventTime_1,PreEventTime_2);
var PostEventPeriod = ee.Filter.date(PostEventTime_1,PostEventTime_2);
// calculate the median of the Pre-event S1 SAR Backscatter Intensity
var PreEventPeriod_asc = ee.Image.cat(asc.filter(PreEventPeriod).median());
var PreEventPeriod_desc = ee.Image.cat(desc.filter(PreEventPeriod).median());
// calculate the median of the Post-event S1 SAR Backscatter Intensity
var PostEventPeriod_asc = ee.Image.cat(asc.filter(PostEventPeriod).median());
var PostEventPeriod_desc = ee.Image.cat(desc.filter(PostEventPeriod).median());
// print out image information (number of images, image file name)
//Note this counts two images acquired on the same day as two separate images, to avoid double counting create very small AOI called "AOI_4_totalimagecount"
var num_asc_pre = asc.filter(PreEventPeriod).filterBounds(AOI_4_totalimagecount);
var num_desc_pre = desc.filter(PreEventPeriod).filterBounds(AOI_4_totalimagecount);
var num_asc_post = asc.filter(PostEventPeriod).filterBounds(AOI_4_totalimagecount);
var num_desc_post = desc.filter(PostEventPeriod).filterBounds(AOI_4_totalimagecount);
var count_asc_pre = num_asc_pre.sort('system:time_start').toList(5000,0).length(); // 5000 controls size of the list
var count_desc_pre = num_desc_pre.sort('system:time_start').toList(5000,0).length(); //
var count_asc_post = num_asc_post.sort('system:time_start').toList(5000,0).length(); //
var count_desc_post = num_desc_post.sort('system:time_start').toList(5000,0).length(); //
print("Number of pre-event ascending images: ", count_asc_pre );
print("Number of pre-event descending images: ", count_desc_pre );
print("Number of post-event ascending images: ", count_asc_post );
print("Number of post-event descending images: ", count_desc_post );
print(num_asc_post,'Post-event ascending'); //print out to check date of image collection
print(num_desc_post,'Post-event descending'); //print out to check date of image collection
// #############################################################################
// ### I_ratio ###
// #############################################################################
// calculate the log ratio (using subtraction since data are in log scale) for Pre- and Post-event S1 SAR Backscatter
var I_ratio_desc = PreEventPeriod_desc.subtract(PostEventPeriod_desc);
I_ratio_desc=I_ratio_desc.clip(AOI);
var I_ratio_asc = PreEventPeriod_asc.subtract(PostEventPeriod_asc);
I_ratio_asc=I_ratio_asc.clip(AOI);
var I_ratio_avg_desc_asc = (I_ratio_asc.add(I_ratio_desc)).divide(2); // calculate the mean backscatter change for ascending and descending scenes combined
// #############################################################################
// ### Define the I_ratio to be used for landslide detection. This can be "I_ratio_desc", "I_ratio_asc", or "I_ratio_avg_desc_asc".
// ### We recommend using "I_ratio_avg_desc_asc" if both ascending and descending data exist. Otherwise use either "I_ratio_desc" or "I_ratio_asc" ###
// #############################################################################
//slope and water mask
var I_ratio_avg_masked = I_ratio_avg_desc_asc.updateMask(mask_slope).updateMask(waterMask);
//slope, curvature, and water mask
//var I_ratio_avg_masked = I_ratio_avg_desc_asc.updateMask(mask_slope).updateMask(mask_curvature).updateMask(waterMask);
// #############################################################################
// ### Calculate I_ratio statistics ###
// #############################################################################
// Calculate percentiles of I_ratio
var I_ratio_Percentiles = I_ratio_avg_masked.reduceRegion({
reducer:ee.Reducer.percentile([90, 95, 99]),
geometry: AOI,
scale: 10,
bestEffort: true
});
print('I_ratio Percentiles',I_ratio_Percentiles)
var I_ratio_90thPtile = I_ratio_avg_masked.gte(ee.Number(I_ratio_Percentiles.get("VH_p90")))
var I_ratio_95thPtile = I_ratio_avg_masked.gte(ee.Number(I_ratio_Percentiles.get("VH_p95")))
var I_ratio_99thPtile = I_ratio_avg_masked.gte(ee.Number(I_ratio_Percentiles.get("VH_p99")))
// #############################################################################
// ### Heatmap and threshold-based landslide detection ###
// #############################################################################
var backscatter_chng_threshold = ee.Number(I_ratio_Percentiles.get("VH_p99")); // backscatter change threshold. This determines the values that correspond to potential landslides. Used for heatmap, polygons, points, and raster of possible landslide locations
print('I_ratio threshold',I_ratio_Percentiles.get("VH_p99"))
// identify possible landslide zones
var pos_slide_zones = I_ratio_avg_masked.gt(backscatter_chng_threshold);
pos_slide_zones = pos_slide_zones.updateMask(pos_slide_zones.neq(0));
// Convert the pos_slide_zones of the thresholded to vectors.
var pos_slide_polygons = pos_slide_zones.addBands(I_ratio_avg_masked).reduceToVectors({
geometry: AOI,
crs: I_ratio_avg_masked.projection(),
scale: 10,
bestEffort: true,
geometryType: 'polygon',
eightConnected: false,
labelProperty: 'zone',
reducer: ee.Reducer.mean()
});
// Convert from polygons to points
var S1_points = imgVH;
var first = ee.Image(S1_points.first()).clip(pos_slide_polygons)
// get image projection
var proj = first.select([0]).projection()
// get coordinates image
var latlon = ee.Image.pixelLonLat().reproject(proj)
// put each lon lat in a list
var coords = latlon.select(['longitude', 'latitude'])
.reduceRegion({
reducer: ee.Reducer.toList(),
geometry: pos_slide_polygons,
scale: 10,
bestEffort: true,
})
// get lat & lon
var lat = ee.List(coords.get('latitude'))
var lon = ee.List(coords.get('longitude'))
// zip them.
var point_list = lon.zip(lat)
// Create points
var mp = ee.Geometry.MultiPoint(point_list)
//******************** Identify points for heatmap to be exported as a KML ***********************//
var mp2 = ee.FeatureCollection(point_list.map(function(p){
var point = ee.Feature(ee.Geometry.Point(p), {})
return point
}))
var points_heatmap = mp2.map(function(feature){
return feature.set('dummy',1);
});
// #############################################################################
// ### Color Palettes for Maps ###
// #############################################################################
// define color palette ('ffffff' is white, 'ff0000' is red, '0000ff' is blue, '000000' is black)
var ColorScale = {min: -10, max: 10, palette: ['0013ff','8178ff','ffffff','ff7e7e','ff0000']}; // for Backscatter Intensity change (I_ratio)
var rgbVis = {min: 0.0, max: 0.18, bands: ['B4', 'B3', 'B2']}; // for sentinel-2
var rgbVis_wClouds = {min: 0.0, max: 2000, bands: ['B4', 'B3', 'B2']}; // for sentinel-2 with clouds
var ColorCurv = {min: -0.02, max:0.02, palette: ['ff0000','ffffff','0000ff']}; // for curvature
var MaskColor = {min: 0, max:1, palette: ['ffffff','0050d6']}; // for curvature
var percentileColor = {min: 0, max:1, palette: ['ffffff','c000d6']}; // for showing percentile
// #############################################################################
// ### Add layers to maps ###
// #############################################################################
//DEM layers
Map.addLayer(hillshade, null, 'NASADEM Hillshade');
//topo lines
var lines = ee.List.sequence(0, 4000, 100)
//var lines = ee.List.sequence(0, 4000, 200)
var contourlines = lines.map(function(line) {
var mycontour = elevation
.convolve(ee.Kernel.gaussian(5, 3))
.subtract(ee.Image.constant(line)).zeroCrossing()
.multiply(ee.Image.constant(line)).toFloat();
return mycontour.mask(mycontour);
})
contourlines = ee.ImageCollection(contourlines).mosaic()
Map.addLayer(contourlines, {min: 0, max: 1000, palette:['000000', '000000']}, 'topo lines', false, 0.5)
//Sentinel-2 optical images
Map.addLayer(S2_PreEvent_median, rgbVis , 'S2 pre-event',false);
Map.addLayer(S2_PostEvent_median, rgbVis, 'S2 post-event',false);
//Possible landslide areas defined based on backscatter change threshold
Map.addLayer(pos_slide_zones.clip(AOI), {min: 0, max: 1, palette: ['ffffff','FF0000']}, 'possible landslide zones',false);
Map.addLayer(I_ratio_avg_masked, ColorScale, 'I_ratio masked', false);
//Map.addLayer(I_ratio_90thPtile, percentileColor, 'I_ratio >= 90th percentile', false);
Map.addLayer(I_ratio_95thPtile, percentileColor, 'I_ratio >= 95th percentile', false);
Map.addLayer(I_ratio_99thPtile, percentileColor, 'I_ratio >= 99th percentile', false, 0.75);
Map.addLayer(Zhangetal2018_Inventory,{},'Zhang et al. (2018) landslide inventory',true,0.9)
// #############################################################################
// ### Generate heatmap visualization in GEE ###
// #############################################################################
//NOTE This often causes GEE to timeout. Until we can optimize the code, we leave this commented out and perform the heatmap calculation in QGIS ***********************//
var radius_heatmap2 = 20; // heatmap kernal radius
function heatmap(mp2,radius_heatmap){
var ptImg = points_heatmap.reduceToImage(['dummy'],ee.Reducer.first()).unmask(0);
var kernel = ee.Kernel.circle(radius_heatmap);
var result = ptImg.convolve(kernel);
return result.clip(AOI); //clip to AOI
//return result.updateMask(result.neq(0)); //uncomment this if want to remove green area
}
var heatmapImg = heatmap(points_heatmap,radius_heatmap2);
var gradientcolormap = ['lightgreen','yellow','red']; //heatmap gradient
var heatmapColor = {
min: 0,
max: .1,
palette: gradientcolormap,
};
//heatmap
Map.addLayer(heatmapImg,heatmapColor,'Heatmap',true, 0.6);
// #############################################################################
// ### save data for export ###
// #############################################################################
Export.image.toDrive({
image: S2_PostEvent_median,
description: 'Hokkaido_S2_PostEvent',
scale: 10,
fileFormat: 'GeoTIFF',
region: AOI
});
Export.image.toDrive({
image: S2_PreEvent_median,
description: 'Hokkaido_S2_PreEvent',
scale: 10,
fileFormat: 'GeoTIFF',
region: AOI
});
Export.image.toDrive({
image: I_ratio_99thPtile,
description: 'Hokkaido_I_ratio>=99thPtile',
scale: 10,
fileFormat: 'GeoTIFF',
region: AOI
});
Export.image.toDrive({
image: I_ratio_avg_masked,
description: 'Hokkaido_I_ratio_masked',
scale: 10,
fileFormat: 'GeoTIFF',
region: AOI
});
Export.table.toDrive({
collection: points_heatmap,
description:'Hokkaido_Points4Heatmap',
fileFormat: 'KML'
});