#!/usr/bin/env python # coding: utf-8 # #### post flood standing water # In[2]: import ee import geemap import geopandas as gpd ee.Authenticate() ee.Initialize() # In[3]: aoi = gpd.read_file('C:/Users/jtrum/Wash scan/data/aoiLuanda.geojson') luanda_geometry = ee.Geometry.Point([13.2344, -8.8115]) luanda_extent = ee.Geometry.Rectangle([12.8, -9.1, 13.6, -8.5]) pt = luanda_geometry AOI = luanda_extent id = "L1a" ''' # return a list of the coordinates of the bounding box of the aoi coords = aoi.bounds.values[0] # Reorder the coordinates to match Earth Engine's convention (west, south, east, north) coords_ee = [coords[0], coords[1], coords[2], coords[1], coords[2], coords[3], coords[0], coords[3], coords[0], coords[1]] # Create a polygon geometry using the reordered bounding box coordinates polygon = ee.Geometry.Polygon(coords_ee) # Convert the geometry to a feature aoi_feature = ee.Feature(polygon) # Set AOI as the geometry of the feature AOI = aoi_feature.geometry() id = 'L1a' ''' srtm = ee.Image("USGS/SRTMGL1_003") slope = ee.Terrain.slope(srtm) grade10 = slope.gt(15) gtGrade10 = grade10.updateMask(grade10.neq(0)) slopeMask = gtGrade10.clip(AOI).unmask(0).subtract(1).multiply(-1) BaseImageStartDate = '2016-11-01' BaseImageEndDate = '2017-02-28' EventStartDate = '2017-03-13' EventEndDate = '2017-03-27' PostEventSearchStartDate = '2017-04-10' PostEventSearchEndDate = '2017-04-24' print('Event started ' + EventStartDate + ' and ended ' + EventEndDate) waterMask = ee.Image("JRC/GSW1_0/GlobalSurfaceWater") # Select transition band from water mask waterMask = waterMask.select('transition') #not transition # Create a blank image blank = ee.Image(0) # Create non-water mask nonWater = blank.addBands(waterMask).unmask().select('transition').eq(0).rename('non_water') # 0 is no change # Define the rescale function def rescale(image): return image.divide(10000) # SE2Baseline SE2Baseline = ee.ImageCollection('COPERNICUS/S2') \ .filterDate(BaseImageStartDate, BaseImageEndDate) \ .filterBounds(AOI) \ .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20)) \ .map(rescale) \ .sort('CLOUDY_PIXEL_PERCENTAGE') \ .first() \ .updateMask(nonWater) # SE2CollectionPostEvent SE2CollectionPostEvent = ee.ImageCollection('COPERNICUS/S2') \ .filterDate(PostEventSearchStartDate, PostEventSearchEndDate) \ .filterBounds(AOI) \ .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30)) \ .map(rescale) \ .sort('CLOUDY_PIXEL_PERCENTAGE') \ .first() \ .updateMask(nonWater) # Get the image id used for post-event properties = ee.String(SE2CollectionPostEvent.get('system:index')) print('Image id used for post-event: ') print(properties.getInfo()) # CloudsPostEvent CloudsPostEvent = SE2CollectionPostEvent.select('QA60').eq(0) # Calculate AWEInshBase AWEInshBase = SE2Baseline.expression( '4 * ((GREEN - SWIR1) - (0.25 * NIR + 2.75 * SWIR2))', { 'GREEN': SE2Baseline.select('B3'), 'SWIR1': SE2Baseline.select('B11'), 'SWIR2': SE2Baseline.select('B12'), 'NIR': SE2Baseline.select('B8') }) # Calculate AWEInshEvent AWEInshEvent = SE2CollectionPostEvent.expression( '4 * ((GREEN - SWIR1) - (0.25 * NIR + 2.75 * SWIR2))', { 'GREEN': SE2CollectionPostEvent.select('B3'), 'SWIR1': SE2CollectionPostEvent.select('B11'), 'SWIR2': SE2CollectionPostEvent.select('B12'), 'NIR': SE2CollectionPostEvent.select('B8') }).updateMask(nonWater) # Compute the threshold for AWEInshBase threshAWEBase = ee.Number(AWEInshBase.reduceRegion( reducer=ee.Reducer.percentile([95]), scale=10, bestEffort=True, geometry=AOI ).get('constant')) # Compute the threshold for AWEInshEvent threshAWEEvent = ee.Number(AWEInshEvent.reduceRegion( reducer=ee.Reducer.percentile([95]), scale=10, bestEffort=True, geometry=AOI ).get('constant')) # Print the thresholds print("AWE base threshold") print(threshAWEBase.getInfo()) print("AWE event threshold") print(threshAWEEvent.getInfo()) # Define highAWEBase highAWEBase = AWEInshBase.gt(threshAWEBase) highAWEEvent = AWEInshEvent.gt(threshAWEEvent) # Define floodedAreas floodedAreas = highAWEBase.eq(0).And(highAWEEvent.eq(1)) floodedAreas2 = highAWEEvent.mask(highAWEEvent).eq(1).And(highAWEBase.mask(highAWEBase).eq(0)) # Define SE2PostEventExtent SE2PostEventExtent = ee.ImageCollection('COPERNICUS/S2') \ .filterDate(PostEventSearchStartDate, PostEventSearchEndDate) \ .filterBounds(AOI) \ .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30)) \ .map(rescale) \ .sort('CLOUDY_PIXEL_PERCENTAGE') \ .first() # Define postEventSceneExtent postEventSceneExtent = SE2PostEventExtent.select('B1').gt(0) # # Define sceneFootprint # sceneFootprint = postEventSceneExtent.addBands(postEventSceneExtent) \ # .reduceToVectors( # reducer=ee.Reducer.mean(), # geometry=AOI, # scale=1000 # ) # # Add sceneFootprint layer to the map # Map.addLayer(sceneFootprint) # In[4]: Map = geemap.Map() Map.addLayer(slope.clip(AOI), {}, 'Slope', False) Map.addLayer(slopeMask.clip(AOI), {}, 'Slope Mask', False) Map.addLayer(slope.clip(AOI), {}, 'Slope', False) Map.addLayer(gtGrade10.clip(AOI), {}, 'Less than grade 10 slope', False) Map.addLayer(AOI, {}, 'AOI', False) # Add layers to the map Map.addLayer(SE2CollectionPostEvent, {'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 0.3}, 'RGB Post event', False) Map.addLayer(SE2Baseline, {'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 0.3}, 'RGB base line', False) # Add AWEInshBase layer to the map Map.addLayer(AWEInshBase, {'min': -6, 'max': 0.4}, 'AWEInsh_Base', False) # Add AWEInshEvent layer to the map Map.addLayer(AWEInshEvent, {'min': -6, 'max': 0.4}, 'AWEInsh_Event', False) # Add AWEInshBase layer to the map Map.addLayer(AWEInshBase, {'min': -6, 'max': 0.4}, 'AWEInsh_Base', False) # Add AWEInshEvent layer to the map Map.addLayer(AWEInshEvent, {'min': -6, 'max': 0.4}, 'AWEInsh_Event', False) # Add clipped AWEInshBase layer to the map Map.addLayer(AWEInshBase.clip(AOI), {}, 'clipped', False) Map.addLayer(highAWEBase.mask(highAWEBase).updateMask(nonWater), {'palette': ['green']}, 'water base', False) Map.addLayer(highAWEEvent.mask(highAWEEvent).updateMask(nonWater), {'palette': ['blue']}, 'water event', False) Map.addLayer(floodedAreas.mask(floodedAreas).updateMask(nonWater), {'palette': ['red']}, 'potential floods', False) Map.addLayer(floodedAreas.mask(floodedAreas).updateMask(nonWater).updateMask(CloudsPostEvent).clip(AOI), {'palette':['red']}, 'Floods, cloud mask', False) Map.addLayer(floodedAreas.mask(floodedAreas).updateMask(nonWater).updateMask(CloudsPostEvent).clip(AOI).updateMask(slopeMask), {'palette':['red']}, 'floods masked slope', False) # center the map Map.centerObject(AOI, 10) Map # In[5]: # return a list of the coordinates of the bounding box of the aoi coords = aoi.bounds.values[0] # Reorder the coordinates to match Earth Engine's convention (west, south, east, north) coords_ee = [coords[0], coords[1], coords[2], coords[1], coords[2], coords[3], coords[0], coords[3], coords[0], coords[1]] # Create a polygon geometry using the reordered bounding box coordinates polygon = ee.Geometry.Polygon(coords_ee) # Convert the geometry to a feature aoi_feature = ee.Feature(polygon) # Set AOI as the geometry of the feature AOI = aoi_feature.geometry() # Define the rescale function def rescale(image): return image.divide(10000) # SE2Baseline SE2Baseline = ee.ImageCollection('COPERNICUS/S2') \ .filterDate(BaseImageStartDate, BaseImageEndDate) \ .filterBounds(AOI) \ .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20)) \ .map(rescale) \ .sort('CLOUDY_PIXEL_PERCENTAGE') \ .first() \ .clip(AOI) # SE2CollectionPostEvent SE2CollectionPostEvent = ee.ImageCollection('COPERNICUS/S2') \ .filterDate(PostEventSearchStartDate, PostEventSearchEndDate) \ .filterBounds(AOI) \ .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30)) \ .map(rescale) \ .sort('CLOUDY_PIXEL_PERCENTAGE') \ .first() \ .clip(AOI) # Calculate AWEInshBase AWEInshBase = SE2Baseline.expression( '4 * ((GREEN - SWIR1) - (0.25 * NIR + 2.75 * SWIR2))', { 'GREEN': SE2Baseline.select('B3'), 'SWIR1': SE2Baseline.select('B11'), 'SWIR2': SE2Baseline.select('B12'), 'NIR': SE2Baseline.select('B8') }) # Calculate AWEInshEvent AWEInshEvent = SE2CollectionPostEvent.expression( '4 * ((GREEN - SWIR1) - (0.25 * NIR + 2.75 * SWIR2))', { 'GREEN': SE2CollectionPostEvent.select('B3'), 'SWIR1': SE2CollectionPostEvent.select('B11'), 'SWIR2': SE2CollectionPostEvent.select('B12'), 'NIR': SE2CollectionPostEvent.select('B8') }) # Compute the threshold for AWEInshBase threshAWEBase = ee.Number(AWEInshBase.reduceRegion( reducer=ee.Reducer.percentile([95]), scale=10, bestEffort=True, geometry=AOI ).get('constant')) # Compute the threshold for AWEInshEvent threshAWEEvent = ee.Number(AWEInshEvent.reduceRegion( reducer=ee.Reducer.percentile([95]), scale=10, bestEffort=True, geometry=AOI ).get('constant')) # Define highAWEBase highAWEBase = AWEInshBase.gt(threshAWEBase) # Define highAWEEvent highAWEEvent = AWEInshEvent.gt(threshAWEEvent) # Define floodedAreas floodedAreas = highAWEBase.eq(0).And(highAWEEvent.eq(1)) # Add layers to the map Map.addLayer(SE2CollectionPostEvent, {'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 0.3}, 'RGB Post event') Map.addLayer(SE2Baseline, {'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 0.3}, 'RGB base line') Map.addLayer(AWEInshBase, {'min': -6, 'max': 0.4}, 'AWEInsh_Base') Map.addLayer(AWEInshEvent, {'min': -6, 'max': 0.4}, 'AWEInsh_Event') Map.addLayer(highAWEBase.mask(highAWEBase), {'palette': ['green']}, 'water base') Map.addLayer(highAWEEvent.mask(highAWEEvent), {'palette': ['blue']}, 'water event') Map.addLayer(floodedAreas.mask(floodedAreas), {'palette': ['red']}, 'potential floods') # Display the map Map # In[7]: coords # In[6]: AOI