{ "cells": [ { "cell_type": "markdown", "metadata": { "id": "s1nFKY6h-g35", "tags": [] }, "source": [ "# Tutorial 2: Natural Hazard Risk Assessment" ] }, { "cell_type": "markdown", "metadata": { "id": "F5jQ4EU--g38" }, "source": [ "In the second tutorial of this week, we are going to use publicly available hazard data and exposure data to do a risk assessment for an area of choice within Europe. More specifically we will look at damage due to wind storms and flooding. \n", "\n", "We will use both Copernicus Land Cover data and OpenStreetMap to estimate the potential damage of natural hazards to the built environment. We will use Copernicus Land Cover data to estimate the damage to specific land-uses, whereas we will use OpenStreetMap to assess the potential damage to the road system." ] }, { "cell_type": "markdown", "metadata": { "id": "K66fASIyP--9" }, "source": [ "### Important before we start\n", "---\n", "Make sure that you save this file before you continue, else you will lose everything. To do so, go to **Bestand/File** and click on **Een kopie opslaan in Drive/Save a Copy on Drive**!\n", "\n", "Now, rename the file into Week5_Tutorial2.ipynb. You can do so by clicking on the name in the top of this screen." ] }, { "cell_type": "markdown", "metadata": { "id": "hkXRFEqH-g38" }, "source": [ "## Learning Objectives\n", "
" ] }, { "cell_type": "markdown", "metadata": { "id": "WDQJZqvL-g39" }, "source": [ "- To know how to download data from the Copernicus Climate Data Store using the `cdsapi` and access it through Python.\n", "- To be able to open and visualize this hazard data.\n", "- To know how to access and open information from the Copernicus Land Monitoring System. Specifically the Corine Land Cover data.\n", "- To understand the basic approach of a natural hazard risk assessment.\n", "- To be able to use the `DamageScanner` to do a damage assessment.\n", "- To interpret and compare the damage estimates." ] }, { "cell_type": "markdown", "metadata": { "id": "4KseyvEh-g39" }, "source": [ "

Tutorial Outline

\n", "
\n", "
" ] }, { "cell_type": "markdown", "metadata": { "id": "ysRSLb-6-g3-" }, "source": [ "## 1.Introducing the packages\n", "
" ] }, { "cell_type": "markdown", "metadata": { "id": "h3aMPYWo-g3-" }, "source": [ "Within this tutorial, we are going to make use of the following packages: \n", "\n", "[**GeoPandas**](https://geopandas.org/) is a Python packagee that extends the datatypes used by pandas to allow spatial operations on geometric types.\n", "\n", "[**OSMnx**](https://osmnx.readthedocs.io/) is a Python package that lets you download geospatial data from OpenStreetMap and model, project, visualize, and analyze real-world street networks and any other geospatial geometries. You can download and model walkable, drivable, or bikeable urban networks with a single line of Python code then easily analyze and visualize them. You can just as easily download and work with other infrastructure types, amenities/points of interest, building footprints, elevation data, street bearings/orientations, and speed/travel time.\n", "\n", "[**xarray**](https://docs.xarray.dev/) is a Python package that allows for easy and efficient use of multi-dimensional arrays.\n", "\n", "[**Matplotlib**](https://matplotlib.org/) is a comprehensive Python package for creating static, animated, and interactive visualizations in Python. Matplotlib makes easy things easy and hard things possible.\n", "\n", "*We will first need to install the missing packages in the cell below. Uncomment them to make sure we can pip install them*" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "xAk93LFz-g3-", "outputId": "b0dca5b3-3894-457e-8f2e-e0188b7b8706" }, "outputs": [], "source": [ "!pip install osmnx\n", "!pip install rioxarray\n", "!pip install cdsapi" ] }, { "cell_type": "markdown", "metadata": { "id": "T28OUcKK-g3_" }, "source": [ "Now we will import these packages in the cell below:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "HNZcbQQ_-g4A", "outputId": "2e8b1214-69e9-40f5-bc84-d02672bf3a5d" }, "outputs": [], "source": [ "import os\n", "import cdsapi\n", "import shapely \n", "import matplotlib\n", "import urllib3\n", "import pyproj\n", "\n", "import osmnx as ox\n", "import numpy as np\n", "import xarray as xr\n", "import geopandas as gpd\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import networkx as nx\n", "\n", "from matplotlib.colors import ListedColormap\n", "from zipfile import ZipFile\n", "from io import BytesIO\n", "from urllib.request import urlopen\n", "from zipfile import ZipFile\n", "from tqdm import tqdm\n", "\n", "urllib3.disable_warnings()" ] }, { "cell_type": "markdown", "metadata": { "id": "TJr_yUdjFMKz" }, "source": [ "### Connect to google drive\n", "---\n", "To be able to read the data from Google Drive, we need to *mount* our Drive to this notebook.\n", "\n", "As you can see in the cell below, make sure that in your **My Drive** folder, you have created a **BigData** folder and within that folder, you have created a **Week5_Data** folder in which you can store the files that are required to run this analysis.\n", "\n", "Please go the URL when its prompted in the box underneath the following cell, and copy the authorization code in that box." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "-Xqbo0nMFTkd", "outputId": "8ce3d2fd-fb29-44f4-b874-a77a9199a41f" }, "outputs": [], "source": [ "from google.colab import drive\n", "drive.mount('/content/gdrive/')\n", "\n", "import sys\n", "sys.path.append(\"/content/gdrive/My Drive/BigData/\")\n", "\n", "data_path = os.path.join('/content/gdrive/My Drive/BigData/','Week5_Data')" ] }, { "cell_type": "markdown", "metadata": { "id": "fqDSHMah-g4A", "tags": [] }, "source": [ "## 2. Downloading and accessing natural hazard data\n", "
" ] }, { "cell_type": "markdown", "metadata": { "id": "1YWXcMRV-g4A" }, "source": [ "We are going to perform a damage assessment using both windstorm data and flood data for Europe.\n", "\n", "### Windstorm Data\n", "\n", "The windstorm data will be downloaded from the [Copernicus Climate Data Store](https://cds.climate.copernicus.eu/). As we have seen during the lecture, and as you can also see by browsing on this website, there is an awful lot of climate data available through this Data Store. As such, it is very valuable to understand how to access and download this information to use within an analysis. To keep things simple, we only download one dataset today: [A winter windstorm](https://cds.climate.copernicus.eu/cdsapp#!/dataset/sis-european-wind-storm-indicators?tab=overview). \n", "\n", "We will do so using an **API**, which is the acronym for application programming interface. It is a software intermediary that allows two applications to talk to each other. APIs are an accessible way to extract and share data within and across organizations. APIs are all around us. Every time you use a rideshare app, send a mobile payment, or change the thermostat temperature from your phone, you’re using an API.\n", "\n", "However, before we can access this **API**, we need to take a few steps. Most importantly, we need to register ourselves on the [Copernicus Climate Data Store](https://cds.climate.copernicus.eu/) portal. To do so, we need to register, as explained in the video clip below:\n", "\n", "\n", "
\n", "\n", "Now, the next step is to access the API. You can now login on the website of the [Copernicus Climate Data Store](https://cds.climate.copernicus.eu/). After you login, you can click on your name in the top right corner of the webpage (next to the login button). On the personal page that has just opened, you will find your user ID (**uid**) and your personal **API**. You need to add those in the cell below to be able to download the windstorm.\n", "\n", "As you can see in the cell below, we download a specific windstorm that has occured on the 28th of October in 2013. This is storm [Carmen (also called St Jude)](https://en.wikipedia.org/wiki/St._Jude_storm). " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "42ejGQJF-g4B", "outputId": "2936be1c-6e11-4d67-91cc-e339af3e661b" }, "outputs": [], "source": [ "uid = XXX\n", "apikey = 'XXX'\n", "\n", "c = cdsapi.Client(key=f\"{uid}:{apikey}\", url=\"https://cds.climate.copernicus.eu/api/v2\")\n", "\n", "c.retrieve(\n", " 'sis-european-wind-storm-indicators',\n", " {\n", " 'variable': 'all',\n", " 'format': 'zip',\n", " 'product': 'windstorm_footprints',\n", " 'year': '2013',\n", " 'month': '10',\n", " 'day': '28',\n", " },\n", " 'Carmen.zip')" ] }, { "cell_type": "markdown", "metadata": { "id": "0RKepWbV-g4B", "tags": [] }, "source": [ "### Flood Data\n", "\n", "The flood data we will extract from a repository maintained by the European Commission Joint Research Centre. We will download river flood hazard maps from their [Flood Data Collection](https://data.jrc.ec.europa.eu/dataset/1d128b6c-a4ee-4858-9e34-6210707f3c81). \n", "\n", "Here we do not need to use an API and we also do not need to register ourselves, so we can download any of the files directly. To do so, we use the `urllib` package." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 363 }, "id": "gdS3cvIv-g4B", "outputId": "233e924d-c368-471b-9590-5b13c9e4e806" }, "outputs": [], "source": [ "## this is the link to the 1/100 flood map for Europe\n", "zipurl = 'https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/FLOODS/EuropeanMaps/floodMap_RP100.zip'\n", "\n", "# and now we open and extract the data\n", "with urlopen(zipurl) as zipresp:\n", " with ZipFile(BytesIO(zipresp.read())) as zfile:\n", " zfile.extractall(data_path)" ] }, { "cell_type": "markdown", "metadata": { "id": "IJrBotZevbZ6" }, "source": [ "The download and zip in the cell above sometimes does not work. If that is indeed the case (e.g., when it seems to remain stuck), download the files manually through the link and upload them in the data folder for this week (as explained at the start of this tutorial.)" ] }, { "cell_type": "markdown", "metadata": { "id": "GCBoCWv2-g4B" }, "source": [ "### Set location to explore\n", "---\n", "Before we continue, we need to specify our location of interest. This should be a province that will have some flooding and relative high wind speeds occuring (else we will find zero damage).\n", "\n", "Specify the region in the cell below by using the `geocode_to_gdf()` function again. I have chosen Gelderland, but feel free to choose a different region. It would be good to double check later on whether there is actually some high wind speeds and flooding in your chosen area." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "bQfZLu6T-g4B" }, "outputs": [], "source": [ "place_name = \"Gelderland, The Netherlands\"\n", "area = ox.geocode_to_gdf(place_name)" ] }, { "cell_type": "markdown", "metadata": { "id": "qES1iP8N-g4C" }, "source": [ "## 3. Explore the natural hazard data\n", "
" ] }, { "cell_type": "markdown", "metadata": { "id": "v89JIHRz-g4C" }, "source": [ "As you can see in the section above, we have downloaded the storm footprint in a zipfile. Let's open the zipfile and load the dataset using the `xarray` package through the `open_dataset()` function." ] }, { "cell_type": "markdown", "metadata": { "id": "I8FlFSFw-g4C" }, "source": [ "### Windstorm Data\n", "---" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "UFcP1DTs-g4C" }, "outputs": [], "source": [ "with ZipFile('Carmen.zip') as zf:\n", " \n", " # Let's get the filename first\n", " file = zf.namelist()[0]\n", " \n", " # And now we can open and select the file within Python\n", " with zf.open(file) as f:\n", " windstorm_europe = xr.open_dataset(f)" ] }, { "cell_type": "markdown", "metadata": { "id": "KFZEcL_C-g4C" }, "source": [ "Let's have a look at the storm we have downloaded!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "windstorm_europe['FX'].plot()" ] }, { "cell_type": "markdown", "metadata": { "id": "Xz30IPy7QAiY" }, "source": [ "
\n", "Question 1: Describe windstorm Carmen. When did this event happen, which areas were most affected? Can you say something about the maximum wind speeds in different areas, based on the plot? And what does FX mean?\n", "
" ] }, { "cell_type": "markdown", "metadata": { "id": "zLb6mFA9-g4D" }, "source": [ "Unfortunately, our data does not have a proper coordinate system defined yet. As such, we will need to use the `rio.write_crs()` function to set the coordinate system to **EPSG:4326** (the standard global coordinate reference system). \n", "\n", "We also need to make sure that the functions will know what the exact parameters are that we have to use for our spatial dimenions (e.g. longitude and latitude). It prefers to be named `x` and `y`. So we use the `rename()` function before we use the `set_spatial_dims()` function." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 291 }, "id": "WqAmuTVo-g4D", "outputId": "0337b286-bc80-4d33-f182-f35d0458686f" }, "outputs": [], "source": [ "windstorm_europe.rio.write_crs(4326, inplace=True)\n", "windstorm_europe = windstorm_europe.rename({'Latitude': 'y','Longitude': 'x'})\n", "windstorm_europe.rio.set_spatial_dims(x_dim=\"x\",y_dim=\"y\", inplace=True)" ] }, { "cell_type": "markdown", "metadata": { "id": "sHxrXrOkRNqA" }, "source": [ "
\n", "Question 2: Climate data is often stored as a netCDF file. Please describe what a netCDF file is. Which information is stored in the netCDF file we have downloaded for the windstorm? What type of metadata does it contain?\n", "
" ] }, { "cell_type": "markdown", "metadata": { "id": "E2cHd5hV-g4D" }, "source": [ "Following, we also make sure it will be in the European coordinate system **EPSG:3035** to ensure we can easily use it together with the other data. To do so, we use the `reproject()` function. You can simple add the number of the coordinate system." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "AqZvdyGk-g4D" }, "outputs": [], "source": [ "windstorm_europe = windstorm_europe.rio.reproject(XXXX)" ] }, { "cell_type": "markdown", "metadata": { "id": "6o_aa1sp-g4D" }, "source": [ "Now we have all the information to clip the windstorm data to our area of interest:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "NRQV7Io3-g4E" }, "outputs": [], "source": [ "windstorm_map = windstorm_europe.rio.clip(area.envelope.values, area.crs)" ] }, { "cell_type": "markdown", "metadata": { "id": "WlWwtOKH-g4E" }, "source": [ "And let's have a look as well by using the `plot()` function. Please note that the legend is in meters per second." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 324 }, "id": "AzolOfUC-g4E", "outputId": "3de9b3d2-8907-48ee-b828-4b03a097a80d" }, "outputs": [], "source": [ "windstorm_map['FX']. XXXX" ] }, { "cell_type": "markdown", "metadata": { "id": "5jQ2_51Z-g4E" }, "source": [ "### Flood data\n", "---" ] }, { "cell_type": "markdown", "metadata": { "id": "btJy-1aR-g4E" }, "source": [ "And similarly, we want to open the flood map. But now we do not have to unzip the file anymore and we can directly open it through using `xarray`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "oRvMEvkm-g4E" }, "outputs": [], "source": [ "flood_map_path = os.path.join(data_path,'floodmap_EFAS_RP100_C.tif')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 328 }, "id": "eEA1EKlt-g4E", "outputId": "1968dec7-b036-4df4-bd9f-b16c41102b33" }, "outputs": [], "source": [ "flood_map = xr.open_dataset(flood_map_path, engine=\"rasterio\")\n", "flood_map" ] }, { "cell_type": "markdown", "metadata": { "id": "O4EIxBK_-g4F" }, "source": [ "And let's make sure we set all the variables and the CRS correctly again to be able to open the data properly. Note that we now use **EPSG:3035**. This is the standard coordinate system for Europe, in meters (instead of degrees)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 291 }, "id": "LxTaPlzW-g4F", "outputId": "8c9d5bcb-3ceb-4628-88fb-0aabe161365f", "tags": [] }, "outputs": [], "source": [ "flood_map.rio.write_crs(3035, inplace=True)\n", "flood_map.rio.set_spatial_dims(x_dim=\"x\",y_dim=\"y\", inplace=True)" ] }, { "cell_type": "markdown", "metadata": { "id": "utiDVMpJ-g4F" }, "source": [ "Now it is pretty difficult to explore the data for our area of interest, so let's clip the flood data. \n", "\n", "We want to clip our flood data to our chosen area. The code, however, is very inefficient and will run into memories issues on Google Colab. As such, we first need to clip it by using a bounding box, followed by the actual clip.\n", "\n", "
\n", "Question 3: Please provide the lines of code below in which you show how you have clipped the flood map to your area.\n", "
\n", "\n", "*A few hints*:\n", "\n", "* carefully read the documentation of the `.clip_box()` function of rioxarray. Which information do you need? \n", "* is the GeoDataFrame of your region (the area GeoDataframe) in the same coordinate system? Perhaps you need to convert it using the `.to_crs()` function. \n", "* how do you get the bounds from your area GeoDataFrame? \n", "* The final step of the clip would be to use the `.rio.clip()` function, using the actual area file and the flood map clipped to the bounding box. Please note that you should **not** use the envelope here, like we did in the previous clip. Here we really want to use the exact geometry values.\n", "\n", "As you will see, we first clip it very efficiently using the bounding box. After that, we do an exact clip." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "PuJ290bGISlL" }, "outputs": [], "source": [ "min_lon = area.to_crs(epsg=3035).bounds.minx.values[0]\n", "min_lat = area.to_crs(epsg=3035).bounds.miny\n", "max_lon = area.to_crs(epsg=3035).bounds\n", "max_lat = area.to_crs(epsg=3035).\n", "\n", "flood_map_area = flood_map.rio.clip_box(minx=.... )\n", "flood_map_area = flood_map_area.rio.clip(area.XXXX.values, area.crs)" ] }, { "cell_type": "markdown", "metadata": { "id": "xNjj8RT--g4F" }, "source": [ "And let's have a look as well. Please note that the legend is in meters." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 324 }, "id": "v_wldK5x-g4F", "outputId": "15c3b526-d094-40e1-dae1-9dbfe57339ea" }, "outputs": [], "source": [ "flood_map_area['band_data'].plot(cmap='Blues',vmax=10)" ] }, { "cell_type": "markdown", "metadata": { "id": "iiOgZMi0-g4F", "tags": [] }, "source": [ "## 4. Download and access Copernicus Land Cover data\n", "
" ] }, { "cell_type": "markdown", "metadata": { "id": "JTp2SMuK-g4F" }, "source": [ "Unfortunately, there is no API option to download the [Corine Land Cover](https://land.copernicus.eu/pan-european/corine-land-cover) data. We will have to download the data from the website first.\n", "\n", "To do so, we will first have to register ourselves again on the website. Please find in the video clip below how to register yourself on the website of the [Copernicus Land Monitoring Service](https://land.copernicus.eu/):\n", "\n", "\n", "\n", "Now click on the Login button in the top right corner to login on the website. There are many interesting datasets on this website, but we just want to download the Corine Land Cover data, and specifically the latest version: [Corine Land Cover 2018](https://land.copernicus.eu/pan-european/corine-land-cover/clc2018?tab=download). To do so, please select the **Corine Land Cover - 100 meter**. Now click on the large green Download button. Your download should start any minute.\n", "\n", "Slightly annoying, the file you have downloaded is double zipped. Its slightly inconvenient to open this through Python and within Google Drive. So let's unzip it twice outside of Python (on your local machine) and then direct yourself to the `DATA` directory within the unzipped file. Here you can find a file called `U2018_CLC2018_V2020_20u1.tif`. Drop this file into this week's data directory, as specified at the start of this tutorial when we mounted our Google Drive." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "e3OV0J1N-g4G" }, "outputs": [], "source": [ "CLC_location = os.path.join(data_path,'U2018_CLC2018_V2020_20u1.tif')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "stGRY3U5-g4G", "outputId": "afaec16f-d144-46ba-af4c-6cc03df68763" }, "outputs": [], "source": [ "CLC = xr.open_dataset(CLC_location, engine=\"rasterio\")" ] }, { "cell_type": "markdown", "metadata": { "id": "EF_na6jnKAvZ" }, "source": [ "Similarly to the flood map data, we need to do a two-stage clip again (like we did before in this tutorial to ensure we get only our area of interest without exceeding our RAM." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "LHBLJONMJ_Zf" }, "outputs": [], "source": [ "CLC_region = CLC.rio.clip_box(minx=.....,)\n", "CLC_region = CLC_region.rio.clip(area.geometry.values,area.crs)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 291 }, "id": "CE5AhRYt-g4G", "outputId": "861397ac-9256-4e1d-a876-c17f0d18cfca" }, "outputs": [], "source": [ "CLC_region = CLC_region.rename({'x': 'lat','y': 'lon'})\n", "CLC_region.rio.set_spatial_dims(x_dim=\"lat\",y_dim=\"lon\", inplace=True)" ] }, { "cell_type": "markdown", "metadata": { "id": "dGDKag4cKTnO" }, "source": [ "And now we create a *color_dict* again, similarly as we did for the raster data in the previous tutorial " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "g7I4fbKs-g4G" }, "outputs": [], "source": [ "CLC_values = [111, 112, 121, 122, 123, 124, 131, 132, 133, 141, 142, 211, 212, 213, 221, 222, 223, 231, 241, 242,\n", " 243, 244, 311, 312, 313, 321, 322, 323, 324, 331, 332, 333, 334, 335, 411, 412, 421, 422, 423, 511, 512, 521, 522, 523]\n", "\n", "CLC_colors = ['#E6004D', '#FF0000', '#CC4DF2', '#CC0000', '#E6CCCC', '#E6CCE6', '#A600CC', '#A64DCC', '#FF4DFF', '#FFA6FF', '#FFE6FF', '#FFFFA8', '#FFFF00', '#E6E600',\n", " '#E68000', '#F2A64D', '#E6A600', '#E6E64D', '#FFE6A6', '#FFE64D', '#E6CC4D', '#F2CCA6', '#80FF00', '#00A600',\n", " '#4DFF00', '#CCF24D', '#A6FF80', '#A6E64D', '#A6F200', '#E6E6E6', '#CCCCCC', '#CCFFCC', '#000000', '#A6E6CC',\n", " '#A6A6FF', '#4D4DFF', '#CCCCFF', '#E6E6FF', '#A6A6E6', '#00CCF2', '#80F2E6', '#00FFA6', '#A6FFE6', '#E6F2FF']" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "Y9Ye19av-g4G" }, "outputs": [], "source": [ "color_dict_raster = dict(zip(CLC_values,CLC_colors))\n", "\n", "# We create a colormar from our list of colors\n", "cm = ListedColormap(CLC_colors)\n", "\n", "# Let's also define the description of each category : 1 (blue) is Sea; 2 (red) is burnt, etc... Order should be respected here ! Or using another dict maybe could help.\n", "labels = np.array(CLC_values)\n", "len_lab = len(labels)\n", "\n", "# prepare normalizer\n", "## Prepare bins for the normalizer\n", "norm_bins = np.sort([*color_dict_raster.keys()]) + 0.5\n", "norm_bins = np.insert(norm_bins, 0, np.min(norm_bins) - 1.0)\n", "\n", "## Make normalizer and formatter\n", "norm = matplotlib.colors.BoundaryNorm(norm_bins, len_lab, clip=True)\n", "fmt = matplotlib.ticker.FuncFormatter(lambda x, pos: labels[norm(x)])" ] }, { "cell_type": "markdown", "metadata": { "id": "5CQoFmdqKcMe" }, "source": [ "And let's plot the Corine Land Cover data for our area of interest" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 650 }, "id": "JIPpIZRh-g4G", "outputId": "59ecd822-c679-42a5-eb5f-ead83aee5102" }, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1,figsize=(14,10))\n", "\n", "CLC_region[\"band_data\"].plot(ax=ax,levels=len(CLC_colors),colors=CLC_colors)" ] }, { "cell_type": "markdown", "metadata": { "id": "AdSaPHzIfNbi" }, "source": [ "
\n", "Question 4: Describe the different land-use classes within your region that you see on the Corine Land Cover map.\n", "
" ] }, { "cell_type": "markdown", "metadata": { "id": "9zbsc7d_-g4G", "tags": [] }, "source": [ "## 5. Perform a damage assessment using Coring Land Cover\n", "
" ] }, { "cell_type": "markdown", "metadata": { "id": "Agxq2HqY-g4H" }, "source": [ "To calculate the potential damage to both windstorms and floods, we use stage-damage curves, which relate the intensity of the hazard to the fraction of maximum damage that can be sustained by a certain land use. As you can see on the Corine Land Cover map that we just plotted, there are a lot of land use classes (44), though not all will suffer damage from either the windstorm or the flood event. For each of the land-use classes a curve and a maximum damage number are assigned.\n", "\n", "To Assess the damage for both the flood and windstorm event, we are going to make use of the [DamageScanner](https://damagescanner.readthedocs.io/en/latest/), which is a tool to calculate potential flood damages based on inundation depth and land use using depth-damage curves in the Netherlands. The DamageScanner was originally developed for the 'Netherlands Later' project [(Klijn et al., 2007)](https://www.rivm.nl/bibliotheek/digitaaldepot/WL_rapport_Overstromingsrisicos_Nederland.pdf). The original land-use classes were based on the Land-Use Scanner in order to evaluate the effect of future land-use change on flood damages. We have tailored the input of the DamageScanner to make sure it can estimate the damages using Corine Land Cover.\n", "\n", "Because the simplicity of the model, we can use this for any raster-based hazard map with some level of intensity. Hence, we can use it for both hazards." ] }, { "cell_type": "markdown", "metadata": { "id": "5m_RAcp_fraF" }, "source": [ "
\n", "Question 5: Describe in your own words what the `DamageScanner()` function does. Please walk us through the different steps. Which inputs do you need to be able to run this damage assessment?\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "jDrTp44Q-g4H" }, "outputs": [], "source": [ "def DamageScanner(landuse_map,inun_map,curve_path,maxdam_path,cellsize=100):\n", " \n", " # load land-use map\n", " landuse = landuse_map.copy()\n", " \n", " # Load inundation map\n", " inundation = inun_map.copy()\n", " \n", " inundation = np.nan_to_num(inundation) \n", "\n", " # Load curves\n", " if isinstance(curve_path, pd.DataFrame):\n", " curves = curve_path.values \n", " elif isinstance(curve_path, np.ndarray):\n", " curves = curve_path\n", "\n", " #Load maximum damages\n", " if isinstance(maxdam_path, pd.DataFrame):\n", " maxdam = maxdam_path.values \n", " elif isinstance(maxdam_path, np.ndarray):\n", " maxdam = maxdam_path\n", " \n", " # Speed up calculation by only considering feasible points\n", " inun = inundation * (inundation>=0) + 0\n", " inun[inun>=curves[:,0].max()] = curves[:,0].max()\n", " waterdepth = inun[inun>0]\n", " landuse = landuse[inun>0]\n", "\n", " # Calculate damage per land-use class for structures\n", " numberofclasses = len(maxdam)\n", " alldamage = np.zeros(landuse.shape[0])\n", " damagebin = np.zeros((numberofclasses, 4,))\n", " for i in range(0,numberofclasses):\n", " n = maxdam[i,0]\n", " damagebin[i,0] = n\n", " wd = waterdepth[landuse==n]\n", " alpha = np.interp(wd,((curves[:,0])),curves[:,i+1])\n", " damage = alpha*(maxdam[i,1]*cellsize)\n", " damagebin[i,1] = sum(damage)\n", " damagebin[i,2] = len(wd)\n", " if len(wd) == 0:\n", " damagebin[i,3] = 0\n", " else:\n", " damagebin[i,3] = np.mean(wd)\n", " alldamage[landuse==n] = damage\n", "\n", " # create pandas dataframe with output\n", " loss_df = pd.DataFrame(damagebin.astype(float),columns=['landuse','losses','area','avg_depth']).groupby('landuse').sum()\n", " \n", " # return output\n", " return loss_df.sum().values[0],loss_df" ] }, { "cell_type": "markdown", "metadata": { "id": "Y7PB8oJz-g4H" }, "source": [ "### Windstorm Damage\n", "---\n", "To estimate the potential damage of our windstorm, we use the vulnerability curves developed by [Yamin et al. (2014)](https://www.sciencedirect.com/science/article/pii/S2212420914000466). Following [Yamin et al. (2014)](https://www.sciencedirect.com/science/article/pii/S2212420914000466), we will apply a sigmoidal vulnerability function satisfying two constraints: (i) a minimum threshold for the occurrence of damage with an upper bound of 100% direct damage; (ii) a high power-law function for the slope, describing an increase in damage with increasing wind speeds. Due to the limited amount of vulnerability curves available for windstorm damage, we will use the damage curve that represents low-rise *reinforced masonry* buildings for all land-use classes that may contain buildings. Obviously, this is a large oversimplification of the real world, but this should be sufficient for this exercise. When doing a proper stand-alone windstorm risk assessment, one should take more effort in collecting the right vulnerability curves for different building types. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "-RxvAEQh-g4H", "tags": [] }, "outputs": [], "source": [ "wind_curves = pd.read_excel(\"https://github.com/ElcoK/BigData_AED/raw/main/week5/damage_curves.xlsx\",sheet_name='wind_curves')\n", "maxdam = pd.read_excel(\"https://github.com/ElcoK/BigData_AED/raw/main/week5/damage_curves.xlsx\",sheet_name='maxdam')" ] }, { "cell_type": "markdown", "metadata": { "id": "uLZ7vl1w-g4H" }, "source": [ "Unfortunately, we run into a *classic* problem when we want to overlay the windstorm data with the Corine Land Cover data. The windstorm data is not only stored in a different coordinate system (we had to convert it from **EPSG:4326** to **EPSG:3035**), it is in a different resolution (**1km** instead of the **100m** of Corine Land Cover). \n", "\n", "Let's first have a look how our clipped data look's like. If you have decided to use Gelderland, you will see that we have 102 columns (our Lattitude/lat) and 74 rows (our Longitude/lon). If you scroll above to our Corine Land Cover data, you see that dimensions are different: 1270 columns (Lattitude/lat/x) and 870 rows (Longitude/lon/y). " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 291 }, "id": "gG2OXOySj8Ra", "outputId": "67135491-52de-4571-f8c9-7b6a74e19b66" }, "outputs": [], "source": [ "windstorm_map" ] }, { "cell_type": "markdown", "metadata": { "id": "igfFBqcK-g4H" }, "source": [ "The first thing we are going to do is try to make sure our data will be in the correct resolution (moving from **1km** to **100m**). To do so, we will use the `rio.reproject()` function. You will see that specify the resolution as **100**. Because **EPSG:3035** is a coordinate system in meters, we can simply use meters to define the resolution. We use the `rio.clip()` function to make sure we clip it again to our area of interest. The function below (`match_rasters`) will do the hard work for us. Please note all the input variables to understand what's happening." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "Kud2CWEDhz1O" }, "outputs": [], "source": [ "def match_rasters(hazard,landuse,haz_crs=3035,lu_crs=3035,resolution=100,hazard_col=['FX']):\n", " \"\"\"\n", " Clips, reprojections, and matches the resolutions of two rasters, `hazard` and `landuse`,\n", " to prepare them for further analysis.\n", "\n", " Parameters\n", " ----------\n", " hazard : xarray.DataArray\n", " A 2D or 3D array containing hazard data.\n", " landuse : xarray.DataArray\n", " A 2D array containing land use data.\n", " haz_crs : int, optional\n", " The CRS of `hazard`. Default is EPSG:3035.\n", " lu_crs : int, optional\n", " The CRS of `landuse`. Default is EPSG:3035.\n", " resolution : float, optional\n", " The desired resolution in meters for both `hazard` and `landuse` after reprojection. Default is 100.\n", " hazard_col : list, optional\n", " A list of column names or indices for the hazard variable. Default is ['FX'].\n", "\n", " Returns\n", " -------\n", " tuple\n", " A tuple containing two xarray.DataArray objects:\n", " - The land use variable with matching resolution and dimensions to the hazard variable.\n", " - The hazard variable clipped to the extent of the land use variable, with matching resolution and dimensions.\n", " \"\"\"\n", " \n", " # Set the crs of the hazard variable to haz_crs\n", " hazard.rio.write_crs(haz_crs, inplace=True)\n", "\n", " # Set the x and y dimensions in the hazard variable to 'x' and 'y' respectively\n", " hazard.rio.set_spatial_dims(x_dim=\"x\",y_dim=\"y\", inplace=True)\n", "\n", " # Reproject the landuse variable from EPSG:4326 to EPSG:3857\n", " landuse = CLC_region.rio.reproject(\"EPSG:3857\",resolution=resolution)\n", "\n", " # Get the minimum longitude and latitude values in the landuse variable\n", " min_lon = landuse.x.min().to_dict()['data']\n", " min_lat = landuse.y.min().to_dict()['data']\n", "\n", " # Get the maximum longitude and latitude values in the landuse variable\n", " max_lon = landuse.x.max().to_dict()['data']\n", " max_lat = landuse.y.max().to_dict()['data']\n", "\n", " # Create a bounding box using the minimum and maximum latitude and longitude values\n", " area = gpd.GeoDataFrame([shapely.box(min_lon,min_lat,max_lon, max_lat)],columns=['geometry'])\n", "\n", " # Set the crs of the bounding box to EPSG:3857\n", " area.crs = 'epsg:3857'\n", "\n", " # Convert the crs of the bounding box to EPSG:4326\n", " area = area.to_crs(f'epsg:{haz_crs}')\n", "\n", " # Clip the hazard variable to the extent of the bounding box\n", " hazard = hazard.rio.clip(area.geometry.values, area.crs)\n", "\n", " # Reproject the hazard variable to EPSG:3857 with the desired resolution\n", " hazard = hazard.rio.reproject(\"EPSG:3857\",resolution=resolution)\n", "\n", " # Clip the hazard variable again to the extent of the bounding box\n", " hazard = hazard.rio.clip(area.geometry.values, area.crs)\n", "\n", " # If the hazard variable has fewer columns and rows than the landuse variable, reproject the landuse variable to match the hazard variable\n", " if (len(hazard.x)len(landuse.x)) & (len(hazard.y)>len(landuse.y)):\n", " hazard = hazard.rio.reproject_match(landuse)\n", "\n", " # return the new landuse and hazard map\n", " return landuse,hazard" ] }, { "cell_type": "markdown", "metadata": { "id": "Vkf6YKPZ-g4I" }, "source": [ "Now let's run the `match_rasters` function and let it do its magic." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "v8NW3c1Q-g4I" }, "outputs": [], "source": [ "CLC_region_wind, windstorm = match_rasters(windstorm_europe,\n", " CLC_region,\n", " haz_crs=3035,\n", " lu_crs=3035,\n", " resolution=100,\n", " hazard_col=['FX'])" ] }, { "cell_type": "markdown", "metadata": { "id": "GgcwJe_6nJip" }, "source": [ "And let's have a look if the two rasters are now the same extend:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 291 }, "id": "vzMbkiSLldlQ", "outputId": "6e73f8b1-33ad-4a7c-b95c-d320b6c75439" }, "outputs": [], "source": [ "CLC_region_wind" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 291 }, "id": "DXnxCBS_ldWg", "outputId": "ec49b756-0fd9-4d49-f9ff-9f68a52bf83c" }, "outputs": [], "source": [ "windstorm" ] }, { "cell_type": "markdown", "metadata": { "id": "6123eX9C-g4J" }, "source": [ "It worked! And to double check, let's also plot it:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 324 }, "id": "Aeay_slW-g4J", "outputId": "11424ad3-2a00-49db-db13-b4db8e73671e" }, "outputs": [], "source": [ "windstorm.FX.plot()" ] }, { "cell_type": "markdown", "metadata": { "id": "JlZF-cs4gSuu" }, "source": [ "
\n", "Question 6: Describe the various steps you have taken to make sure that the windstorm map is now exactly the same extent as the corine land cover map. Feel free to include lines of code in your answer and also describe the different functions you have used along the way.\n", "
" ] }, { "cell_type": "markdown", "metadata": { "id": "LW158xPh-g4J" }, "source": [ "Now its finally time to do our damage assessment! To do so, we need to convert our data to `numpy.arrays()` to do our calculation:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "QZIzWIeP-g4J" }, "outputs": [], "source": [ "landuse_map = CLC_region_wind['band_data'].to_numpy()[0,:,:]\n", "wind_map = windstorm['FX'].to_numpy()[0,:,:]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "aBqqRqbkmA1Y", "outputId": "709d6c91-4ad9-4e27-e6a9-e6201df32dc7" }, "outputs": [], "source": [ "wind_map.shape" ] }, { "cell_type": "markdown", "metadata": { "id": "J9QHyhSU-g4J" }, "source": [ "And remember that our windstorm data was stored in **m/s**. Hence, we need to convert it to **km/h**:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "GqdUCXD_-g4J" }, "outputs": [], "source": [ "wind_map_kmh = wind_map*XXX" ] }, { "cell_type": "markdown", "metadata": { "id": "ln7NqRB1-g4J" }, "source": [ "And now let's run the DamageScanner to obtain the damage results" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "y_g0pj1h-g4J", "tags": [] }, "outputs": [], "source": [ "wind_damage_CLC = DamageScanner(landuse_map,wind_map_kmh,wind_curves,maxdam)[1]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 1000 }, "id": "-6DbFD_JeFA1", "outputId": "fb251350-8885-4dde-e665-893fa04cde6e" }, "outputs": [], "source": [ "wind_damage_CLC" ] }, { "cell_type": "markdown", "metadata": { "id": "5UNySYvk-g4J", "tags": [] }, "source": [ "### Flood Damage\n", "---\n", "To Assess the flood damage, we are again going to make use of the [DamageScanner](https://damagescanner.readthedocs.io/en/latest/). The Corine Land Cover data is widely used in European flood risk assessments. As such, we can simply make use of pre-developed curves. We are using the damage curves as developed by Huizinga et al. (2007). Again, let's first load the maximum damages and the depth-damage curves:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "ua2xyAGW-g4J" }, "outputs": [], "source": [ "flood_curves = pd.read_excel(\"https://github.com/ElcoK/BigData_AED/raw/main/week5/damage_curves.xlsx\",sheet_name='flood_curves',engine='openpyxl')\n", "maxdam = pd.read_excel(\"https://github.com/ElcoK/BigData_AED/raw/main/week5/damage_curves.xlsx\",sheet_name='maxdam')" ] }, { "cell_type": "markdown", "metadata": { "id": "HT54wRvs-g4K" }, "source": [ "And convert our data to `numpy.arrays()` to do our calculation:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "qzXKNmg2-g4K" }, "outputs": [], "source": [ "landuse_map = CLC_region['band_data'].to_numpy()\n", "flood_map = flood_map_area['band_data'].to_numpy()" ] }, { "cell_type": "markdown", "metadata": { "id": "ttGra99k-g4K" }, "source": [ "And now let's run the DamageScanner to obtain the damage results" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "2qL8UATu-g4K" }, "outputs": [], "source": [ "flood_damage_CLC = DamageScanner(landuse_map,flood_map,flood_curves,maxdam)[1]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 1000 }, "id": "WnWu6AMUeFA2", "outputId": "a16faddd-7a66-40e4-f103-9a1c1f723735" }, "outputs": [], "source": [ "flood_damage_CLC" ] }, { "cell_type": "markdown", "metadata": { "id": "hGWhg9fogvM3" }, "source": [ "
\n", "Question 7: Describe the results of the flood and wind damage assessments. Do you notice any differences between the outcomes? Do you observe specific land-use classes that are severely damaged?\n", "
" ] }, { "cell_type": "markdown", "metadata": { "id": "E0ohOKwd-g4K" }, "source": [ "## 6. Perform a damage assessment of the road network using OpenStreetMap\n", "
" ] }, { "cell_type": "markdown", "metadata": { "id": "bKFmTKpj-g4K" }, "source": [ "Generally, wind damage does not cause much damage to roads. There will be clean-up cost of the trees that will fall on the roads, but structural damage is rare. As such, we will only do a flood damage assessment for the road network of our region.\n", "\n", "To do so, we first need to extract the roads again. We will use the `graph_from_place()` function again to do so. However, the area will be to large to extract roads, so we will focus our analysis on the main network." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "CUqFG7AD-g4K" }, "outputs": [], "source": [ "cf = '[\"highway\"~\"trunk|motorway|primary|secondary\"]'\n", "G = ox.graph_from_place(place_name, network_type=\"drive\", custom_filter=cf)" ] }, { "cell_type": "markdown", "metadata": { "id": "modIJTEz-g4K" }, "source": [ "And convert the road network to a `geodataframe`, as done in the previous tutorial as well." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "sxOBriok-g4K", "outputId": "6ddc37c2-d8d1-4222-f989-08a2f9d08d43" }, "outputs": [], "source": [ "roads = gpd.GeoDataFrame(nx.to_pandas_edgelist(G))\n", "roads.highway = roads.highway.astype('str')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In case the above does not work, you can continue the assignment by using the code below (make sure you remove the hashtags to run it). " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#from urllib import request\n", "# remote_url = 'https://github.com/ElcoK/BigData_AED/raw/main/week5/kampen_roads.gpkg'\n", "# file = 'kampen_roads.gpkg'\n", "# \n", "# #request.urlretrieve(remote_url, file)\n", "# roads = gpd.GeoDataFrame.from_file('kampen_roads.gpkg')" ] }, { "cell_type": "markdown", "metadata": { "id": "VIaMGLxA-g4K" }, "source": [ "And lets have a look at the data:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 495 }, "id": "dx_299FS-g4L", "outputId": "ffdab479-6a25-4794-da3a-cbacf2accaa4" }, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1,figsize=(12,10))\n", "\n", "\n", "roads.plot(column='highway',legend=True,ax=ax,legend_kwds={'loc': 'lower right'});\n", "\n", "\n", "# remove the ax labels\n", "ax.set_xticks([])\n", "ax.set_yticks([])\n", "ax.set_axis_off()" ] }, { "cell_type": "markdown", "metadata": { "id": "PSGo7dC3-g4L" }, "source": [ "It is actually quite inconvenient to have all these lists in the data for when we want to do the damage assessment. Let's clean this up a bit. To do so, we first make sure that all the lists are represented as actual lists, and not lists wrapped within a string." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "S_LZSRI6-g4L" }, "outputs": [], "source": [ "roads.highway = roads.highway.apply(lambda x: x.strip('][').split(', '))" ] }, { "cell_type": "markdown", "metadata": { "id": "cwQKiRDd-g4L" }, "source": [ "Now we just need to grab the first element of each of the lists." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "qe86tcET-g4L" }, "outputs": [], "source": [ "roads.highway = roads.highway.apply(lambda x: x[0] if isinstance(x, list) else x)\n", "roads.highway = roads.highway.str.replace(\"'\",\"\")" ] }, { "cell_type": "markdown", "metadata": { "id": "TkyDDDIP-g4L" }, "source": [ "And let's have a look whether this worked:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 495 }, "id": "f5qPSBmq-g4L", "outputId": "6ac152f2-14f1-4e6c-c53a-9d68db347ce6" }, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1,figsize=(12,10))\n", "\n", "roads.plot(column='highway',legend=True,ax=ax,legend_kwds={'loc': 'upper left','ncol':1});\n", "\n", "# remove the ax labels\n", "ax.set_xticks([])\n", "ax.set_yticks([])\n", "ax.set_axis_off()" ] }, { "cell_type": "markdown", "metadata": { "id": "u0crazq8iQjQ" }, "source": [ "
\n", "Question 8: Upload a figure of the cleaned road network (e.g. in which you do not see any of the listed road types anymore)\n", "
" ] }, { "cell_type": "markdown", "metadata": { "id": "J63sExRp-g4L" }, "source": [ "Nice! now let's start with the damage calculation. As you already have may have noticed, our data is now not stored in raster format, but in vector format. One way to deal with this issue is to convert our vector data to raster data, but we will lose a lot of information and detail. As such, we will perform the damage assessment on the road elements, using the xarray flood map.\n", "\n", "Let's start with preparing the flood data into vector format:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "uHmaZFXV-g4L" }, "outputs": [], "source": [ "# get the mean values\n", "flood_map_vector = flood_map_area['band_data'].to_dataframe().reset_index()\n", "\n", "# create geometry values and drop lat lon columns\n", "flood_map_vector['geometry'] = [shapely.points(x) for x in list(zip(flood_map_vector['x'],flood_map_vector['y']))]\n", "flood_map_vector = flood_map_vector.drop(['x','y','band','spatial_ref'],axis=1)\n", "\n", "# drop all non values to reduce size\n", "flood_map_vector = flood_map_vector.loc[~flood_map_vector['band_data'].isna()].reset_index(drop=True)\n", "\n", "# and turn them into squares again:\n", "flood_map_vector.geometry= shapely.buffer(flood_map_vector.geometry,distance=100/2,cap_style='square').values" ] }, { "cell_type": "markdown", "metadata": { "id": "sKb-ig4Q-g4M" }, "source": [ "And let's plot the results:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 306 }, "id": "LL3YU6r1-g4M", "outputId": "e0a9e61f-e376-436c-a11a-fa58651bf15a" }, "outputs": [], "source": [ "gpd.GeoDataFrame(flood_map_vector.copy()).plot(column='band_data',cmap='Blues',vmax=5,linewidth=0)" ] }, { "cell_type": "markdown", "metadata": { "id": "XBsxnhjN-g4M" }, "source": [ "We will need a bunch of functions to make sure we can do our calculations. They are specified below. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "T-XfgGLB-g4M" }, "outputs": [], "source": [ "def reproject(df_ds,current_crs=\"epsg:4326\",approximate_crs = \"epsg:3035\"):\n", " geometries = df_ds['geometry']\n", " coords = shapely.get_coordinates(geometries)\n", " transformer=pyproj.Transformer.from_crs(current_crs, approximate_crs,always_xy=True)\n", " new_coords = transformer.transform(coords[:, 0], coords[:, 1])\n", " \n", " return shapely.set_coordinates(geometries.copy(), np.array(new_coords).T) \n", "\n", "def buffer_assets(assets,buffer_size=100):\n", " assets['buffered'] = shapely.buffer(assets.geometry.values,buffer_size)\n", " return assets\n", "\n", "def overlay_hazard_assets(df_ds,assets):\n", "\n", " #overlay \n", " hazard_tree = shapely.STRtree(df_ds.geometry.values)\n", " if (shapely.get_type_id(assets.iloc[0].geometry) == 3) | (shapely.get_type_id(assets.iloc[0].geometry) == 6):\n", " return hazard_tree.query(assets.geometry,predicate='intersects') \n", " else:\n", " return hazard_tree.query(assets.buffered,predicate='intersects')\n", " \n", "def get_damage_per_asset(asset,df_ds,assets):\n", " # find the exact hazard overlays:\n", " get_hazard_points = df_ds.iloc[asset[1]['hazard_point'].values].reset_index()\n", " get_hazard_points = get_hazard_points.loc[shapely.intersects(get_hazard_points.geometry.values,assets.iloc[asset[0]].geometry)]\n", "\n", " asset_geom = assets.iloc[asset[0]].geometry\n", "\n", " maxdam_asset = 100\n", " hazard_intensity = np.arange(0,10,0.1) \n", " fragility_values = np.arange(0,1,0.01) \n", " \n", " if len(get_hazard_points) == 0:\n", " return asset[0],0\n", " else:\n", " get_hazard_points['overlay_meters'] = shapely.length(shapely.intersection(get_hazard_points.geometry.values,asset_geom))\n", " return asset[0],np.sum((np.interp(get_hazard_points.band_data.values,hazard_intensity,fragility_values))*get_hazard_points.overlay_meters*maxdam_asset)" ] }, { "cell_type": "markdown", "metadata": { "id": "og2Bkcv--g4M" }, "source": [ "Now we need to make sure that the road data is the same coordinate system. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "NvtDuspN-g4M", "outputId": "8d920938-ab07-42f9-bb88-34483e751c3f" }, "outputs": [], "source": [ "roads.geometry = reproject(roads)" ] }, { "cell_type": "markdown", "metadata": { "id": "4JT25WTv-g4M" }, "source": [ "And we can now overlay the roads with the flood data" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "8rtBYbX_-g4M" }, "outputs": [], "source": [ "overlay_roads = pd.DataFrame(overlay_hazard_assets(flood_map_vector,buffer_assets(roads)).T,columns=['asset','hazard_point'])" ] }, { "cell_type": "markdown", "metadata": { "id": "s82DyD_y-g4M" }, "source": [ "And estimate the damages" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "LGqPFklh-g4N", "outputId": "207937c9-dcc5-41a5-ebc4-76aa03003022" }, "outputs": [], "source": [ "collect_output = []\n", "for asset in tqdm(overlay_roads.groupby('asset'),total=len(overlay_roads.asset.unique()),\n", " desc='polyline damage calculation for'):\n", " collect_output.append(get_damage_per_asset(asset,flood_map_vector,roads))\n", " \n", "damaged_roads = roads.merge(pd.DataFrame(collect_output,columns=['index','damage']),\n", " left_index=True,right_on='index')[['highway','geometry','damage']]" ] }, { "cell_type": "markdown", "metadata": { "id": "RFGZxWl7i7pQ" }, "source": [ "
\n", "Question 9: Describe the various steps we have taken to perform the damage assessment on the road network. How is this approach different compared to the raster-based approach? Highlight the differences you find most important. Include any line of code you may want to include to make your story clear.\n", "
" ] }, { "cell_type": "markdown", "metadata": { "id": "B5jpsbyC-g4N" }, "source": [ "And let's plot the results" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 514 }, "id": "n25j-3wG-g4N", "outputId": "b926a8e7-7e51-4434-f3b8-d61e6278bc24" }, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1,figsize=(12,10))\n", "\n", "damaged_roads.plot(column='damage',cmap='Reds',ax=ax);" ] }, { "cell_type": "markdown", "metadata": { "id": "bTfmvwW2jchB" }, "source": [ "
\n", "Question 10: Describe the most severely damaged parts of the road network. Use Google Maps to identify these roads. Are you surprised by the results?\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "SG3FSqsLeFA5" }, "outputs": [], "source": [] } ], "metadata": { "colab": { "provenance": [], "toc_visible": true }, "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.10.6" }, "vscode": { "interpreter": { "hash": "f323064ae63d54ed8d769390a968e914fbf7abacffc63e116cd2e04a08ed2d24" } } }, "nbformat": 4, "nbformat_minor": 4 }