{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "#### post flood standing water" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "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.\n", " warnings.warn(\n", "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:\n", "\n", "import os\n", "os.environ['USE_PYGEOS'] = '0'\n", "import geopandas\n", "\n", "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).\n", " import geopandas as gpd\n" ] } ], "source": [ "import ee\n", "import geemap\n", "import geopandas as gpd\n", "\n", "ee.Authenticate()\n", "ee.Initialize()" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Event started 2017-03-13 and ended 2017-03-27\n", "Image id used for post-event: \n", "20170416T090551_20170416T093109_T33LTK\n", "AWE base threshold\n", "-1.719208923005133\n", "AWE event threshold\n", "-1.7345601910828032\n" ] } ], "source": [ "aoi = gpd.read_file('C:/Users/jtrum/Wash scan/data/aoiLuanda.geojson')\n", "luanda_geometry = ee.Geometry.Point([13.2344, -8.8115]) \n", "luanda_extent = ee.Geometry.Rectangle([12.8, -9.1, 13.6, -8.5])\n", "pt = luanda_geometry\n", "AOI = luanda_extent\n", "id = \"L1a\"\n", "\n", "'''\n", "# return a list of the coordinates of the bounding box of the aoi\n", "coords = aoi.bounds.values[0]\n", "# Reorder the coordinates to match Earth Engine's convention (west, south, east, north)\n", "coords_ee = [coords[0], coords[1], coords[2], coords[1], coords[2], coords[3], coords[0], coords[3], coords[0], coords[1]]\n", "# Create a polygon geometry using the reordered bounding box coordinates\n", "polygon = ee.Geometry.Polygon(coords_ee)\n", "# Convert the geometry to a feature\n", "aoi_feature = ee.Feature(polygon)\n", "# Set AOI as the geometry of the feature\n", "AOI = aoi_feature.geometry()\n", "id = 'L1a'\n", "'''\n", "\n", "srtm = ee.Image(\"USGS/SRTMGL1_003\")\n", "slope = ee.Terrain.slope(srtm)\n", "grade10 = slope.gt(15)\n", "gtGrade10 = grade10.updateMask(grade10.neq(0))\n", "slopeMask = gtGrade10.clip(AOI).unmask(0).subtract(1).multiply(-1)\n", "\n", "\n", "\n", "BaseImageStartDate = '2016-11-01'\n", "BaseImageEndDate = '2017-02-28'\n", "EventStartDate = '2017-03-13'\n", "EventEndDate = '2017-03-27'\n", "PostEventSearchStartDate = '2017-04-10'\n", "PostEventSearchEndDate = '2017-04-24'\n", "print('Event started ' + EventStartDate + ' and ended ' + EventEndDate)\n", "\n", "waterMask = ee.Image(\"JRC/GSW1_0/GlobalSurfaceWater\")\n", "\n", "# Select transition band from water mask\n", "waterMask = waterMask.select('transition') #not transition\n", "\n", "# Create a blank image\n", "blank = ee.Image(0)\n", "\n", "# Create non-water mask\n", "nonWater = blank.addBands(waterMask).unmask().select('transition').eq(0).rename('non_water') # 0 is no change \n", "# Define the rescale function\n", "def rescale(image):\n", " return image.divide(10000)\n", "\n", "# SE2Baseline\n", "SE2Baseline = ee.ImageCollection('COPERNICUS/S2') \\\n", " .filterDate(BaseImageStartDate, BaseImageEndDate) \\\n", " .filterBounds(AOI) \\\n", " .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20)) \\\n", " .map(rescale) \\\n", " .sort('CLOUDY_PIXEL_PERCENTAGE') \\\n", " .first() \\\n", " .updateMask(nonWater)\n", "\n", "# SE2CollectionPostEvent\n", "SE2CollectionPostEvent = ee.ImageCollection('COPERNICUS/S2') \\\n", " .filterDate(PostEventSearchStartDate, PostEventSearchEndDate) \\\n", " .filterBounds(AOI) \\\n", " .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30)) \\\n", " .map(rescale) \\\n", " .sort('CLOUDY_PIXEL_PERCENTAGE') \\\n", " .first() \\\n", " .updateMask(nonWater)\n", "\n", "# Get the image id used for post-event\n", "properties = ee.String(SE2CollectionPostEvent.get('system:index'))\n", "print('Image id used for post-event: ')\n", "print(properties.getInfo())\n", "\n", "# CloudsPostEvent\n", "CloudsPostEvent = SE2CollectionPostEvent.select('QA60').eq(0)\n", "\n", "# Calculate AWEInshBase\n", "AWEInshBase = SE2Baseline.expression(\n", " '4 * ((GREEN - SWIR1) - (0.25 * NIR + 2.75 * SWIR2))', {\n", " 'GREEN': SE2Baseline.select('B3'),\n", " 'SWIR1': SE2Baseline.select('B11'),\n", " 'SWIR2': SE2Baseline.select('B12'),\n", " 'NIR': SE2Baseline.select('B8')\n", "})\n", "\n", "# Calculate AWEInshEvent\n", "AWEInshEvent = SE2CollectionPostEvent.expression(\n", " '4 * ((GREEN - SWIR1) - (0.25 * NIR + 2.75 * SWIR2))', {\n", " 'GREEN': SE2CollectionPostEvent.select('B3'),\n", " 'SWIR1': SE2CollectionPostEvent.select('B11'),\n", " 'SWIR2': SE2CollectionPostEvent.select('B12'),\n", " 'NIR': SE2CollectionPostEvent.select('B8')\n", "}).updateMask(nonWater)\n", "\n", "\n", "# Compute the threshold for AWEInshBase\n", "threshAWEBase = ee.Number(AWEInshBase.reduceRegion(\n", " reducer=ee.Reducer.percentile([95]),\n", " scale=10,\n", " bestEffort=True,\n", " geometry=AOI\n", ").get('constant'))\n", "\n", "# Compute the threshold for AWEInshEvent\n", "threshAWEEvent = ee.Number(AWEInshEvent.reduceRegion(\n", " reducer=ee.Reducer.percentile([95]),\n", " scale=10,\n", " bestEffort=True,\n", " geometry=AOI\n", ").get('constant'))\n", "\n", "# Print the thresholds\n", "print(\"AWE base threshold\")\n", "print(threshAWEBase.getInfo())\n", "print(\"AWE event threshold\")\n", "print(threshAWEEvent.getInfo())\n", "\n", "# Define highAWEBase\n", "highAWEBase = AWEInshBase.gt(threshAWEBase)\n", "highAWEEvent = AWEInshEvent.gt(threshAWEEvent)\n", "\n", "# Define floodedAreas\n", "floodedAreas = highAWEBase.eq(0).And(highAWEEvent.eq(1))\n", "floodedAreas2 = highAWEEvent.mask(highAWEEvent).eq(1).And(highAWEBase.mask(highAWEBase).eq(0))\n", "\n", "# Define SE2PostEventExtent\n", "SE2PostEventExtent = ee.ImageCollection('COPERNICUS/S2') \\\n", " .filterDate(PostEventSearchStartDate, PostEventSearchEndDate) \\\n", " .filterBounds(AOI) \\\n", " .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30)) \\\n", " .map(rescale) \\\n", " .sort('CLOUDY_PIXEL_PERCENTAGE') \\\n", " .first()\n", "\n", "# Define postEventSceneExtent\n", "postEventSceneExtent = SE2PostEventExtent.select('B1').gt(0)\n", "\n", "# # Define sceneFootprint\n", "# sceneFootprint = postEventSceneExtent.addBands(postEventSceneExtent) \\\n", "# .reduceToVectors(\n", "# reducer=ee.Reducer.mean(),\n", "# geometry=AOI,\n", "# scale=1000\n", "# )\n", "\n", "# # Add sceneFootprint layer to the map\n", "# Map.addLayer(sceneFootprint)\n", "\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "2fce57f6974448b49f438857d3f4e65b", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Map(center=[-8.800130028331832, 13.200000000000014], controls=(WidgetControl(options=['position', 'transparent…" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Map = geemap.Map()\n", "Map.addLayer(slope.clip(AOI), {}, 'Slope', False)\n", "Map.addLayer(slopeMask.clip(AOI), {}, 'Slope Mask', False)\n", "Map.addLayer(slope.clip(AOI), {}, 'Slope', False)\n", "Map.addLayer(gtGrade10.clip(AOI), {}, 'Less than grade 10 slope', False)\n", "Map.addLayer(AOI, {}, 'AOI', False)\n", "# Add layers to the map\n", "Map.addLayer(SE2CollectionPostEvent, {'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 0.3}, 'RGB Post event', False)\n", "Map.addLayer(SE2Baseline, {'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 0.3}, 'RGB base line', False)\n", "# Add AWEInshBase layer to the map\n", "Map.addLayer(AWEInshBase, {'min': -6, 'max': 0.4}, 'AWEInsh_Base', False)\n", "# Add AWEInshEvent layer to the map\n", "Map.addLayer(AWEInshEvent, {'min': -6, 'max': 0.4}, 'AWEInsh_Event', False)\n", "# Add AWEInshBase layer to the map\n", "Map.addLayer(AWEInshBase, {'min': -6, 'max': 0.4}, 'AWEInsh_Base', False)\n", "# Add AWEInshEvent layer to the map\n", "Map.addLayer(AWEInshEvent, {'min': -6, 'max': 0.4}, 'AWEInsh_Event', False)\n", "# Add clipped AWEInshBase layer to the map\n", "Map.addLayer(AWEInshBase.clip(AOI), {}, 'clipped', False)\n", "Map.addLayer(highAWEBase.mask(highAWEBase).updateMask(nonWater), {'palette': ['green']}, 'water base', False)\n", "Map.addLayer(highAWEEvent.mask(highAWEEvent).updateMask(nonWater), {'palette': ['blue']}, 'water event', False)\n", "Map.addLayer(floodedAreas.mask(floodedAreas).updateMask(nonWater), {'palette': ['red']}, 'potential floods', False)\n", "Map.addLayer(floodedAreas.mask(floodedAreas).updateMask(nonWater).updateMask(CloudsPostEvent).clip(AOI), {'palette':['red']}, 'Floods, cloud mask', False)\n", "Map.addLayer(floodedAreas.mask(floodedAreas).updateMask(nonWater).updateMask(CloudsPostEvent).clip(AOI).updateMask(slopeMask), {'palette':['red']}, 'floods masked slope', False)\n", "# center the map\n", "Map.centerObject(AOI, 10)\n", "Map" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "f8866462a08d4ea6aa0cb9bfe1c47012", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Map(bottom=17531.0, center=[-9.259360434570048, 15.792257819652129], controls=(WidgetControl(options=['positio…" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# return a list of the coordinates of the bounding box of the aoi\n", "coords = aoi.bounds.values[0]\n", "\n", "# Reorder the coordinates to match Earth Engine's convention (west, south, east, north)\n", "coords_ee = [coords[0], coords[1], coords[2], coords[1], coords[2], coords[3], coords[0], coords[3], coords[0], coords[1]]\n", "\n", "# Create a polygon geometry using the reordered bounding box coordinates\n", "polygon = ee.Geometry.Polygon(coords_ee)\n", "\n", "# Convert the geometry to a feature\n", "aoi_feature = ee.Feature(polygon)\n", "\n", "# Set AOI as the geometry of the feature\n", "AOI = aoi_feature.geometry()\n", "\n", "# Define the rescale function\n", "def rescale(image):\n", " return image.divide(10000)\n", "\n", "# SE2Baseline\n", "SE2Baseline = ee.ImageCollection('COPERNICUS/S2') \\\n", " .filterDate(BaseImageStartDate, BaseImageEndDate) \\\n", " .filterBounds(AOI) \\\n", " .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20)) \\\n", " .map(rescale) \\\n", " .sort('CLOUDY_PIXEL_PERCENTAGE') \\\n", " .first() \\\n", " .clip(AOI)\n", "\n", "# SE2CollectionPostEvent\n", "SE2CollectionPostEvent = ee.ImageCollection('COPERNICUS/S2') \\\n", " .filterDate(PostEventSearchStartDate, PostEventSearchEndDate) \\\n", " .filterBounds(AOI) \\\n", " .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30)) \\\n", " .map(rescale) \\\n", " .sort('CLOUDY_PIXEL_PERCENTAGE') \\\n", " .first() \\\n", " .clip(AOI)\n", "\n", "# Calculate AWEInshBase\n", "AWEInshBase = SE2Baseline.expression(\n", " '4 * ((GREEN - SWIR1) - (0.25 * NIR + 2.75 * SWIR2))', {\n", " 'GREEN': SE2Baseline.select('B3'),\n", " 'SWIR1': SE2Baseline.select('B11'),\n", " 'SWIR2': SE2Baseline.select('B12'),\n", " 'NIR': SE2Baseline.select('B8')\n", "})\n", "\n", "# Calculate AWEInshEvent\n", "AWEInshEvent = SE2CollectionPostEvent.expression(\n", " '4 * ((GREEN - SWIR1) - (0.25 * NIR + 2.75 * SWIR2))', {\n", " 'GREEN': SE2CollectionPostEvent.select('B3'),\n", " 'SWIR1': SE2CollectionPostEvent.select('B11'),\n", " 'SWIR2': SE2CollectionPostEvent.select('B12'),\n", " 'NIR': SE2CollectionPostEvent.select('B8')\n", "})\n", "\n", "# Compute the threshold for AWEInshBase\n", "threshAWEBase = ee.Number(AWEInshBase.reduceRegion(\n", " reducer=ee.Reducer.percentile([95]),\n", " scale=10,\n", " bestEffort=True,\n", " geometry=AOI\n", ").get('constant'))\n", "\n", "# Compute the threshold for AWEInshEvent\n", "threshAWEEvent = ee.Number(AWEInshEvent.reduceRegion(\n", " reducer=ee.Reducer.percentile([95]),\n", " scale=10,\n", " bestEffort=True,\n", " geometry=AOI\n", ").get('constant'))\n", "\n", "# Define highAWEBase\n", "highAWEBase = AWEInshBase.gt(threshAWEBase)\n", "\n", "# Define highAWEEvent\n", "highAWEEvent = AWEInshEvent.gt(threshAWEEvent)\n", "\n", "# Define floodedAreas\n", "floodedAreas = highAWEBase.eq(0).And(highAWEEvent.eq(1))\n", "\n", "# Add layers to the map\n", "Map.addLayer(SE2CollectionPostEvent, {'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 0.3}, 'RGB Post event')\n", "Map.addLayer(SE2Baseline, {'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 0.3}, 'RGB base line')\n", "Map.addLayer(AWEInshBase, {'min': -6, 'max': 0.4}, 'AWEInsh_Base')\n", "Map.addLayer(AWEInshEvent, {'min': -6, 'max': 0.4}, 'AWEInsh_Event')\n", "Map.addLayer(highAWEBase.mask(highAWEBase), {'palette': ['green']}, 'water base')\n", "Map.addLayer(highAWEEvent.mask(highAWEEvent), {'palette': ['blue']}, 'water event')\n", "Map.addLayer(floodedAreas.mask(floodedAreas), {'palette': ['red']}, 'potential floods')\n", "\n", "# Display the map\n", "Map\n" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "array([12.9925 , -9.34425918, 13.63397804, -8.63636811])" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "coords" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
    • type:Polygon
          • 0:12.992500001000167
          • 1:-9.344259181999917
          • 0:13.633978038000066
          • 1:-9.344259181999917
          • 0:13.633978038000066
          • 1:-8.636368108999932
          • 0:12.992500001000167
          • 1:-8.636368108999932
          • 0:12.992500001000167
          • 1:-9.344259181999917
" ], "text/plain": [ "ee.Geometry({\n", " \"functionInvocationValue\": {\n", " \"functionName\": \"Feature.geometry\",\n", " \"arguments\": {\n", " \"feature\": {\n", " \"functionInvocationValue\": {\n", " \"functionName\": \"Feature\",\n", " \"arguments\": {\n", " \"geometry\": {\n", " \"functionInvocationValue\": {\n", " \"functionName\": \"GeometryConstructors.Polygon\",\n", " \"arguments\": {\n", " \"coordinates\": {\n", " \"constantValue\": [\n", " [\n", " [\n", " 12.992500001000167,\n", " -9.344259181999917\n", " ],\n", " [\n", " 13.633978038000066,\n", " -9.344259181999917\n", " ],\n", " [\n", " 13.633978038000066,\n", " -8.636368108999932\n", " ],\n", " [\n", " 12.992500001000167,\n", " -8.636368108999932\n", " ],\n", " [\n", " 12.992500001000167,\n", " -9.344259181999917\n", " ]\n", " ]\n", " ]\n", " },\n", " \"evenOdd\": {\n", " \"constantValue\": true\n", " }\n", " }\n", " }\n", " }\n", " }\n", " }\n", " }\n", " }\n", " }\n", "})" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "AOI" ] } ], "metadata": { "kernelspec": { "display_name": "wash", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.18" } }, "nbformat": 4, "nbformat_minor": 2 }