Notebook

post flood standing water¶

In [2]:

import ee
import geemap
import geopandas as gpd

ee.Authenticate()
ee.Initialize()

c:\Users\jtrum\miniconda3\envs\wash\lib\site-packages\geopandas\_compat.py:124: UserWarning: The Shapely GEOS version (3.11.2-CAPI-1.17.2) is incompatible with the GEOS version PyGEOS was compiled with (3.10.4-CAPI-1.16.2). Conversions between both will be slow.
  warnings.warn(
C:\Users\jtrum\AppData\Local\Temp\ipykernel_26400\430109623.py:3: DeprecationWarning: Shapely 2.0 is installed, but because PyGEOS is also installed, GeoPandas still uses PyGEOS by default. However, starting with version 0.14, the default will switch to Shapely. To force to use Shapely 2.0 now, you can either uninstall PyGEOS or set the environment variable USE_PYGEOS=0. You can do this before starting the Python process, or in your code before importing geopandas:

import os
os.environ['USE_PYGEOS'] = '0'
import geopandas

In the next release, GeoPandas will switch to using Shapely by default, even if PyGEOS is installed. If you only have PyGEOS installed to get speed-ups, this switch should be smooth. However, if you are using PyGEOS directly (calling PyGEOS functions on geometries from GeoPandas), this will then stop working and you are encouraged to migrate from PyGEOS to Shapely 2.0 (https://shapely.readthedocs.io/en/latest/migration_pygeos.html).
  import geopandas as gpd

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)

Event started 2017-03-13 and ended 2017-03-27
Image id used for post-event: 
20170416T090551_20170416T093109_T33LTK
AWE base threshold
-1.719208923005133
AWE event threshold
-1.7345601910828032

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

Out[4]:

Map(center=[-8.800130028331832, 13.200000000000014], controls=(WidgetControl(options=['position', 'transparent…

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

Out[5]:

Map(bottom=17531.0, center=[-9.259360434570048, 15.792257819652129], controls=(WidgetControl(options=['positio…

In [7]:

coords

Out[7]:

array([12.9925    , -9.34425918, 13.63397804, -8.63636811])

In [6]:

AOI

Out[6]:

Polygon (5 vertices)
- type:Polygon
- coordinates: List (1 element)
  - 0: List (5 elements)
    - 0: [12.992500001000167, -9.344259181999917]
      - 0:12.992500001000167
      - 1:-9.344259181999917
    - 1: [13.633978038000066, -9.344259181999917]
      - 0:13.633978038000066
      - 1:-9.344259181999917
    - 2: [13.633978038000066, -8.636368108999932]
      - 0:13.633978038000066
      - 1:-8.636368108999932
    - 3: [12.992500001000167, -8.636368108999932]
      - 0:12.992500001000167
      - 1:-8.636368108999932
    - 4: [12.992500001000167, -9.344259181999917]
      - 0:12.992500001000167
      - 1:-9.344259181999917