diff --git a/TAA2/tutorial.ipynb b/TAA2/tutorial.ipynb index d4e2abe..ffe4bc0 100644 --- a/TAA2/tutorial.ipynb +++ b/TAA2/tutorial.ipynb @@ -476,7 +476,7 @@ "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 4: Please provide the lines of code below in which you show how you have clipped the flood map to your area.\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", @@ -701,7 +701,7 @@ "metadata": {}, "source": [ "
\n", - "Question 5: Describe the different land-use classes within your region that you see on the Corine Land Cover map. Do you see any dominant land-use classes? \n", + "Question 4: Describe the different land-use classes within your region that you see on the Corine Land Cover map. Do you see any dominant land-use classes? \n", "
" ] }, @@ -768,7 +768,7 @@ "metadata": {}, "source": [ "
\n", - "Question 1: To make sure we understand which area you focus on, please submit the figure that outlines your area.\n", + "Question 5: To make sure we understand which area you focus on, please submit the figure that outlines your area.\n", "
" ] }, @@ -871,7 +871,7 @@ "metadata": {}, "source": [ "
\n", - "Question 2: Please provide in the answer box in Canvas the code that you used to make sure that all land uses are now registered within the landuse column.\n", + "Question 6: Please provide in the answer box in Canvas the code that you used to make sure that all land uses are now registered within the landuse column.\n", "
" ] }, @@ -953,7 +953,7 @@ "In case anything is missing, add them to the color_dict dictionairy and re-run that cell. \n", "\n", "
\n", - "Question 3: Show us in Canvas (i) which land-use categories you had to add, and (ii) how your final color dictionary looks like.\n", + "Question 7: Show us in Canvas (i) which land-use categories you had to add, and (ii) how your final color dictionary looks like.\n", "
\n", "\n", "```{tip}\n", @@ -1028,7 +1028,7 @@ "metadata": {}, "source": [ "
\n", - "Question 4: Please upload a figure of your land-use map, using OpenStreetMap. \n", + "Question 8: Please upload a figure of your land-use map, using OpenStreetMap. \n", "
" ] }, @@ -1268,7 +1268,7 @@ "metadata": {}, "source": [ "
\n", - "Question 5: In the rasterization process, we use the `.make_geocube()` function. Please elaborate on the following: i)why is it important to specify the right coordinate system? What could happen if you choose the wrong coordinate system? ii) which resolution did you choose and why? iii)Why did the first result did not give us the right output with the correct colors? How did you solve this? \n", + "Question 9: In the rasterization process, we use the `.make_geocube()` function. Please elaborate on the following: i)why is it important to specify the right coordinate system? What could happen if you choose the wrong coordinate system? ii) which resolution did you choose and why? iii)Why did the first result did not give us the right output with the correct colors? How did you solve this? \n", "
" ] }, @@ -1344,46 +1344,468 @@ "
" ] }, + { + "cell_type": "markdown", + "id": "5dfabe80-3c8f-4e3a-a411-a29349c2d75b", + "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", + "id": "ed5b5f52-bbad-42a6-acd5-1ea0248591f8", + "metadata": { + "id": "5m_RAcp_fraF" + }, + "source": [ + "
\n", + "Question 10: 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, - "id": "d965589e-5317-4111-a9ad-06094f76569f", - "metadata": {}, + "id": "94d33d66-6bf6-43b2-8c24-3e1a1069af2b", + "metadata": { + "id": "jDrTp44Q-g4H" + }, "outputs": [], "source": [ - "! @TASK: Copy paste from Tutorial 2" + "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", + "id": "146c6b79-ce68-436b-a7fb-8a889528f513", + "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, - "id": "e4c4ad75-a555-414d-bb64-6af5f10a2eb1", - "metadata": {}, + "id": "e0e63daa-636b-4a1e-ba6c-5bb78e56f290", + "metadata": { + "id": "-RxvAEQh-g4H", + "tags": [] + }, "outputs": [], "source": [ - "! @TASK: Let students also perform risk assessment based on rasterized OSM land use data?? Code needs to be written for that. Or think about a really good question that we can ask here about this?" + "wind_curves = pd.read_excel(\"https://github.com/VU-IVM/UNIGIS_ProgrammingGIS/tree/main/TAA2/damage_curves.xlsx\",sheet_name='wind_curves')\n", + "maxdam = pd.read_excel(\"https://github.com/VU-IVM/UNIGIS_ProgrammingGIS/tree/main/TAA2/damage_curves.xlsx\",sheet_name='maxdam')" + ] + }, + { + "cell_type": "markdown", + "id": "beb4f312-c4fb-4169-b563-e5878dfd95e2", + "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, - "id": "38ac963e-2f8a-4e84-a313-96c3d00f2b1c", - "metadata": {}, + "id": "57bc4345-b1e5-45c7-8aa2-95ef39b1ea78", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 291 + }, + "id": "gG2OXOySj8Ra", + "outputId": "67135491-52de-4571-f8c9-7b6a74e19b66" + }, "outputs": [], - "source": [] + "source": [ + "windstorm_map" + ] + }, + { + "cell_type": "markdown", + "id": "e1688219-1647-442e-8ed7-e4277ca75a16", + "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, - "id": "38d46699-d759-48e3-a5ef-09d45ec8e41a", + "id": "d960af61-b324-459e-9e4c-da2f70bf7ecf", + "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", + "id": "9f9ccb0d-2644-45ea-b435-79d83cd12a8e", + "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, + "id": "d5c3409f-5c41-4ffa-b881-fa166eec2521", + "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", + "id": "aa2fdc04-2327-47d2-b43a-521dc95b6882", + "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, + "id": "7ca69a83-350b-4d0e-adea-d6a00b9dfc38", + "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, + "id": "0b0216bd-6497-4a0a-8903-31d46599a854", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 291 + }, + "id": "DXnxCBS_ldWg", + "outputId": "ec49b756-0fd9-4d49-f9ff-9f68a52bf83c" + }, + "outputs": [], + "source": [ + "windstorm" + ] + }, + { + "cell_type": "markdown", + "id": "50f3a8c5-aac5-4def-b79f-8e20be1cc822", + "metadata": { + "id": "6123eX9C-g4J" + }, + "source": [ + "It worked! And to double check, let's also plot it:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "58536a6a-6e5b-40cd-8d14-dc453d781405", + "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", + "id": "8e445352-2695-4d54-9abe-e3215817f355", + "metadata": { + "id": "JlZF-cs4gSuu" + }, + "source": [ + "
\n", + "Question 11: 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", + "id": "4959f306-df08-4548-8470-9083c814b203", + "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, + "id": "f0f82a7c-f2b1-43dd-bcb6-4a9d24738d75", + "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, + "id": "53d83a65-c63b-427b-b950-2f4674fe8b82", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "aBqqRqbkmA1Y", + "outputId": "709d6c91-4ad9-4e27-e6a9-e6201df32dc7" + }, + "outputs": [], + "source": [ + "wind_map.shape" + ] + }, + { + "cell_type": "markdown", + "id": "aca9bbe2-95f2-4b0d-9fb7-bb0ffc94ecef", + "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, + "id": "11ddca33-7360-4efb-84f2-7613fd360ca3", + "metadata": { + "id": "GqdUCXD_-g4J" + }, + "outputs": [], + "source": [ + "wind_map_kmh = wind_map*XXX" + ] + }, + { + "cell_type": "markdown", + "id": "53b5ac7e-e817-49e2-b588-614e997adf8a", + "metadata": { + "id": "ln7NqRB1-g4J" + }, + "source": [ + "And now let's run the DamageScanner to obtain the damage results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "31da4031-fdd4-485f-af99-30fa49fdabe0", + "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, + "id": "a717a4a2-8d6b-436b-b855-bebc7835d0b3", + "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": "code", + "execution_count": null, + "id": "e4c4ad75-a555-414d-bb64-6af5f10a2eb1", "metadata": {}, "outputs": [], - "source": [] + "source": [ + "! @TASK: Let students also perform risk assessment based on rasterized OSM land use data?? Code needs to be written for that. Or think about a really good question that we can ask here about this?" + ] }, { "cell_type": "code", "execution_count": null, - "id": "786c9da2-4257-4749-9fbf-fe0f94ef61a0", + "id": "38d46699-d759-48e3-a5ef-09d45ec8e41a", "metadata": {}, "outputs": [], "source": [] @@ -1406,6 +1828,16 @@ "---" ] }, + { + "cell_type": "markdown", + "id": "4f42ba03-9c17-465d-97eb-77df502dd8ee", + "metadata": { + "id": "WMOSa6JF47DS" + }, + "source": [ + "There is a lot more data to extract from OpenStreetMap besides land-use information. Let's extract some building data. To do so, we use the *\"building\"* tag." + ] + }, { "cell_type": "code", "execution_count": null, @@ -1419,115 +1851,1173 @@ }, { "cell_type": "markdown", - "id": "68ae51a7-1c32-4d96-bff0-cb67ce796fc8", + "id": "772d95a6-c2ce-48a6-a7d0-bc1f4b77c323", "metadata": {}, "source": [ - "There is a lot more data to extract from OpenStreetMap besides land-use information. Let's extract some building data. To do so, we use the *\"building\"* tag." + "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). If you decide to use the data as specified below, also change the map at the start to 'Kampen'." ] }, { "cell_type": "code", - "execution_count": null, - "id": "26a60806-c927-4d0c-ab04-dc07c9c7a688", + "execution_count": 64, + "id": "9a7e508a-b20a-416c-8f2a-c8dd36a60389", "metadata": {}, "outputs": [], "source": [ - "buildings.head()" + "# remote_url = https://github.com/VU-IVM/UNIGIS_ProgrammingGIS/tree/main/TAA2'\n", + "# file = 'kampen_buildings.gpkg'\n", + "# \n", + "# #request.urlretrieve(remote_url, file)\n", + "# buildings = gpd.GeoDataFrame.from_file('kampen_buildings.gpkg')" + ] + }, + { + "cell_type": "markdown", + "id": "53e45bdb-36a3-4efd-b0e4-111f5fd858fe", + "metadata": { + "id": "bon7osXA47DT" + }, + "source": [ + "Now let's see what information is actually extracted:" ] }, { "cell_type": "code", "execution_count": null, - "id": "f9f8340c-21e7-4e9b-996d-e2234ced09d9", + "id": "26a60806-c927-4d0c-ab04-dc07c9c7a688", "metadata": {}, "outputs": [], "source": [ - "! @TASK: copy-paste cells from tutorial 1" + "buildings.head()" ] }, { "cell_type": "markdown", - "id": "90a6ba47-e148-48e1-a0a2-e7e2c75c4257", - "metadata": {}, + "id": "2f7f7ff1-ed38-4dca-b2ad-f09472566879", + "metadata": { + "id": "p8dhZx1n47DT" + }, "source": [ - "### Analyze and visualize building stock\n", - "---" + "As you notice in the output of the cell above, there are many columns which just contain \"NaN\". And there even seem to be to many columns to even visualize properly in one view.\n", + "\n", + "Let's check what information is collected for the different buildings:" ] }, { "cell_type": "code", "execution_count": null, - "id": "75442fe8-0e84-4315-a591-3c705d280357", - "metadata": {}, + "id": "b78b1253-16b0-4513-8903-966b1816715f", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "executionInfo": { + "elapsed": 35, + "status": "ok", + "timestamp": 1675087861529, + "user": { + "displayName": "RA Odongo", + "userId": "17326618845752559881" + }, + "user_tz": -60 + }, + "id": "W3ZpsIkO47DT", + "outputId": "385616f8-25fe-4675-f5f0-16d7fa9b0df7" + }, "outputs": [], "source": [ - "fig,ax = plt.subplots(1,1,figsize=(5,18))\n", - "\n", - "building_year.plot(kind='barh',ax=ax)\n", - "\n", - "ax.tick_params(axis='y', which='major', labelsize=7)" + "buildings.columns" ] }, { - "cell_type": "code", - "execution_count": null, - "id": "73f5051e-797f-4a42-8680-cd9767086d08", - "metadata": {}, - "outputs": [], + "cell_type": "markdown", + "id": "772be8a9-7ab6-4996-9311-e165fcc6c8ed", + "metadata": { + "id": "sH7NTKENA4WO" + }, "source": [ - "! @TASK: copy-paste cells from tutorial 1" + "
\n", + "Question 12: Let's have a look at the extracted building information. Please describe in your own words the information it contains. Is there specific information that suprises you to see, and do you think anything is missing that you expected? \n", + "
" ] }, { "cell_type": "markdown", - "id": "c59fa9dc-68f1-4d39-8adc-f917a0af3a04", - "metadata": {}, + "id": "d7077034-4cfe-436a-b5b2-64fc77f24ae5", + "metadata": { + "id": "Z37JRRc747DT" + }, "source": [ - "### Extracting roads from OpenStreetMap\n", - "---" + "One interesting column is called `start_date`. This shows the building year per building. \n", + "\n", + "Let's explore this year of building a bit more.\n", + "\n", + "First, it would be interesting to get an idea how many buildings are build in each year through using the `value_counts()` function. Normally, that functions ranks the values in descending order (high to low). We are more interested in how this has developed over time. So we use the `sort_index()` function to sort the values by year. Add these two functions in the cell below." ] }, { "cell_type": "code", "execution_count": null, - "id": "d34bcba3-d295-49d8-9948-054723c2e58f", - "metadata": {}, + "id": "e27c9a71-68b4-4adc-a1b4-52ff94275a6f", + "metadata": { + "executionInfo": { + "elapsed": 924, + "status": "ok", + "timestamp": 1675087884735, + "user": { + "displayName": "RA Odongo", + "userId": "17326618845752559881" + }, + "user_tz": -60 + }, + "id": "lNSe826i47DT" + }, "outputs": [], "source": [ - "! @TASK: copy-paste cells from tutorial 1" + "building_year = buildings.start_date. XXXX" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "9b68047e-2370-4584-b6fb-d97013ce2e9c", - "metadata": {}, - "outputs": [], - "source": [] - }, { "cell_type": "markdown", - "id": "508022f5-46d7-4c69-91d2-46711d6a514c", - "metadata": {}, + "id": "e9e9c436-fd01-46f7-8f73-44eccafe7816", + "metadata": { + "id": "X6b_T5xs47DU" + }, "source": [ - "## 7. Perform a damage assessment of the road network using OpenStreetMap\n", - "
" + "There is not better way to further explore this years than through plotting it. Don't forget to add things such as a x label, y label and title. Have a look at some of the matplotlib [tutorials](https://matplotlib.org/stable/tutorials/introductory/quick_start.html). Note that you need to look at the code that also uses subplots and where they use the `ax` option." ] }, { "cell_type": "code", "execution_count": null, - "id": "6be31e8c-34d3-42c7-a4bd-76d681702f83", - "metadata": {}, + "id": "35a6ee4b-f664-4955-9781-86d56920b3e4", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 1000 + }, + "executionInfo": { + "elapsed": 1636, + "status": "ok", + "timestamp": 1675087889714, + "user": { + "displayName": "RA Odongo", + "userId": "17326618845752559881" + }, + "user_tz": -60 + }, + "id": "fmk6m78I47DU", + "outputId": "78f3ea6a-284d-495f-b596-12db32305128" + }, + "outputs": [], + "source": [ + "fig,ax = plt.subplots(1,1,figsize=(5,18))\n", + "\n", + "building_year.plot(kind='barh',ax=ax)\n", + "\n", + "ax.tick_params(axis='y', which='major', labelsize=7)" + ] + }, + { + "cell_type": "markdown", + "id": "7e5ce5b6-f7b6-47a4-b4b7-2cac7cc933bf", + "metadata": { + "id": "AwRyhvNeBn0G" + }, + "source": [ + "
\n", + "Question 13: Please upload a figure that shows the development of building stock over the years in your region of interest. Make sure it contains all the necessary elements (labels on the axis, title, etc.)\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "089f6011-87fd-464a-ad8b-52ad810d8b75", + "metadata": { + "id": "wrcM1p-m47DU" + }, + "source": [ + "What we also noticed is that quite some buildings are identified as 'yes'. This is not very useful as it does not really say much about the use of the building. \n", + "\n", + "Let's see for how many buildings this is the case: " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ca163a97-54bd-4050-bea6-2721d024a1eb", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "executionInfo": { + "elapsed": 205, + "status": "ok", + "timestamp": 1675087945407, + "user": { + "displayName": "RA Odongo", + "userId": "17326618845752559881" + }, + "user_tz": -60 + }, + "id": "_4eKR2Bo47DU", + "outputId": "54e29c44-9717-4eee-984f-9555659317be" + }, + "outputs": [], + "source": [ + "buildings.building.value_counts()" + ] + }, + { + "cell_type": "markdown", + "id": "502fcb4a-cc93-46fa-bc1a-5dad7f4438bc", + "metadata": {}, + "source": [ + "As you have seen from the `value_counts` function, there are quite a few buildings with only very few tags. You could either consider to not include them in your plot at all (for example by using the `isin` function or the `query` function, see also [here](https://stackoverflow.com/questions/12096252/use-a-list-of-values-to-select-rows-from-a-pandas-dataframe)), or rename them, similar to how you named the natural land cover classes for the land-use map. " + ] + }, + { + "cell_type": "markdown", + "id": "22598d8d-a73c-41b7-befd-f8abafc86162", + "metadata": { + "id": "vF6WJ9_k47DU" + }, + "source": [ + "Now let's visualize the buildings again. We need to create a similar color dictionary as we did for the land-use categories. Now its up to you to make it!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "739a2389-8302-4e4f-8624-2beb0225ecfd", + "metadata": { + "executionInfo": { + "elapsed": 226, + "status": "ok", + "timestamp": 1675087956546, + "user": { + "displayName": "RA Odongo", + "userId": "17326618845752559881" + }, + "user_tz": -60 + }, + "id": "RYRxKXxd47DU" + }, + "outputs": [], + "source": [ + "color_dict = { 'yes' : \"#f1134b\", \n", + " 'house':'#f13013', \n", + " 'industrial':'#0f045c',\n", + " 'farm':'#fcfcb9', \n", + " 'bungalow':'#f13013',\n", + " 'service':'#CB8DDB' }" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "54672469-3961-4b67-860e-17f9320d1a62", + "metadata": { + "executionInfo": { + "elapsed": 232, + "status": "ok", + "timestamp": 1675087958726, + "user": { + "displayName": "RA Odongo", + "userId": "17326618845752559881" + }, + "user_tz": -60 + }, + "id": "YjRwC-K-47DU" + }, + "outputs": [], + "source": [ + "# Remove multiple keys from dictionary\n", + "color_dict = {key: color_dict[key]\n", + " for key in color_dict if key not in list(set(color_dict)-set(buildings.building.unique()))}\n", + "\n", + "map_dict = dict(zip(color_dict.keys(),[x for x in range(len(color_dict))]))\n", + "buildings['col_landuse'] =buildings.building.apply(lambda x: color_dict[x])" + ] + }, + { + "cell_type": "markdown", + "id": "0bbaef6d-6b28-45f4-ace9-6a5c9dec0f2f", + "metadata": { + "id": "wGh7rbnB47DU" + }, + "source": [ + "And plot the figure in the same manner!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7dad1b9f-855a-468e-9090-25d7ace2990d", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 616 + }, + "executionInfo": { + "elapsed": 3651, + "status": "ok", + "timestamp": 1675087966347, + "user": { + "displayName": "RA Odongo", + "userId": "17326618845752559881" + }, + "user_tz": -60 + }, + "id": "oDfamvIV47DV", + "outputId": "dbcee552-e955-4d2b-ddda-909ab4265105" + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(1, 1,figsize=(12,10))\n", + "\n", + "# add color scheme\n", + "color_scheme_map = list(color_dict.values())\n", + "cmap = LinearSegmentedColormap.from_list(name='landuse',\n", + " colors=color_scheme_map) \n", + "\n", + "# and plot the land-use map.\n", + "buildings.plot(color=buildings['col_landuse'],ax=ax,linewidth=0)\n", + "\n", + "# remove the ax labels\n", + "ax.set_xticks([])\n", + "ax.set_yticks([])\n", + "ax.set_axis_off()\n", + "\n", + "# add a legend:\n", + "legend_elements = []\n", + "for iter_,item in enumerate(color_dict):\n", + " legend_elements.append(Patch(facecolor=color_scheme_map[iter_],label=item)) \n", + "\n", + "ax.legend(handles=legend_elements,edgecolor='black',facecolor='#fefdfd',prop={'size':12},loc=(1.02,0.2)) \n", + "\n", + "# add a title\n", + "ax.set_title(place_name,fontweight='bold')" + ] + }, + { + "cell_type": "markdown", + "id": "7fc1437d-8deb-4569-a3ee-9270676f0e41", + "metadata": { + "id": "IJalDDJvB6Yd" + }, + "source": [ + "
\n", + "Question 14: Please upload a figure of your building stock map of your region of interest. Make sure that the interpretation is clear. If necessary, merge multiple categories into one (i.e., when some categories only contain 1 or 2 buildings).\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "c59fa9dc-68f1-4d39-8adc-f917a0af3a04", + "metadata": {}, + "source": [ + "### Extracting roads from OpenStreetMap\n", + "---" + ] + }, + { + "cell_type": "markdown", + "id": "83eb8f23-6b4e-47ce-a1df-12ca10e90897", + "metadata": { + "id": "HpWnKfIk47DV" + }, + "source": [ + "Let's continue (and end) this tutorial with the core data in OpenStreetMap (it is even in the name): roads!\n", + "\n", + "Now, instead of using tags, we want to identify what type of roads we would like to extract. Let's first only extract roads that can be used to drive.\n", + "\n", + "The `graph_from_place()` function returns a `NetworkX` Graph element. You can read more about these graph elements in the introduction page of [NetworkX](https://networkx.org/documentation/stable/reference/introduction.html)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "73036bbd-08b1-4d2d-bc14-9efd48d4ec76", + "metadata": { + "executionInfo": { + "elapsed": 13403, + "status": "ok", + "timestamp": 1675088179196, + "user": { + "displayName": "RA Odongo", + "userId": "17326618845752559881" + }, + "user_tz": -60 + }, + "id": "MpjxUadp47DV" + }, + "outputs": [], + "source": [ + "G = ox.graph_from_place(place_name, network_type=\"drive\")" + ] + }, + { + "cell_type": "markdown", + "id": "6c1596cb-9a1b-4a67-87cd-cee4f687ad11", + "metadata": { + "id": "NA_HAmHx47DV" + }, + "source": [ + "Unfortunately, it is bit difficult to easily view all the roads within such a Graph element. To be able to explore the data, we are going to convert it to a `Geopandas GeoDataFrame`, using the `to_pandas_edgelist()` function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1b52effc-4eaf-44b5-b031-2481d2125197", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "executionInfo": { + "elapsed": 15, + "status": "ok", + "timestamp": 1675088179197, + "user": { + "displayName": "RA Odongo", + "userId": "17326618845752559881" + }, + "user_tz": -60 + }, + "id": "4dMXZb_v47DV", + "outputId": "e615d5f8-6b96-4a32-f114-ce38b1481b70" + }, + "outputs": [], + "source": [ + "roads = gpd.GeoDataFrame(nx.to_pandas_edgelist(G))" + ] + }, + { + "cell_type": "markdown", + "id": "7a3a191c-0d36-4bb5-8ed7-4f12bb05e1f5", + "metadata": { + "id": "yds5NBlF47DV" + }, + "source": [ + "In some cases, roads are classified with more than one category. If that is the case, they are captured within a `list`. To overcome this issue, we specify that we want the entire `highway` column as a `string` dtype." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d5f76461-aac1-4534-8a97-2abc841107b4", + "metadata": { + "executionInfo": { + "elapsed": 294, + "status": "ok", + "timestamp": 1675088191403, + "user": { + "displayName": "RA Odongo", + "userId": "17326618845752559881" + }, + "user_tz": -60 + }, + "id": "NxTg1jb847DV" + }, + "outputs": [], + "source": [ + "roads.highway = roads.highway.astype('str')" + ] + }, + { + "cell_type": "markdown", + "id": "5080124a-c7f4-45f6-aced-a27f8793ef59", + "metadata": { + "id": "rbzl5JPR47DW" + }, + "source": [ + "Now we can create a plot to see how the road network is configured." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b0c2f581-4182-471e-823b-0280181494e0", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 575 + }, + "executionInfo": { + "elapsed": 888, + "status": "ok", + "timestamp": 1675088196301, + "user": { + "displayName": "RA Odongo", + "userId": "17326618845752559881" + }, + "user_tz": -60 + }, + "id": "uyo3Ggpc47DW", + "outputId": "e58c915a-b085-498c-e7ea-8d62f237d417" + }, + "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", + "id": "82a98a1c-e8d3-4f84-80df-5e8577a3458b", + "metadata": { + "id": "M18fuTWM47DW", + "tags": [] + }, + "source": [ + "It would also be interesting to explore the network a little but more interactively. **OSMnx** has a function called `plot_graph_folium()`, which allow us to use the [folium](https://python-visualization.github.io/folium/quickstart.html#Getting-Started) package to plot data interactively on a map. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dcc7ab55-f4ce-4706-83a0-5b944c2bfd7d", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 866 + }, + "executionInfo": { + "elapsed": 1720, + "status": "ok", + "timestamp": 1675088204394, + "user": { + "displayName": "RA Odongo", + "userId": "17326618845752559881" + }, + "user_tz": -60 + }, + "id": "lpWHR0Ez47DW", + "outputId": "38be3bea-b399-4e8f-a081-210baa559ef7" + }, + "outputs": [], + "source": [ + "m1 = ox.plot_graph_folium(G, popup_attribute=\"highway\", weight=2, color=\"#8b0000\")\n", + "m1" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "274d172a-62df-4095-987d-1de387104e8d", + "metadata": { + "id": "gm8NIAR947DW" + }, + "source": [ + "One of the exiting things we can do with this data is that we can compute and plot routes between two points on a map.\n", + "\n", + "Let's first select two random start and end points from the graph and compute the shortest route between them through using the `shortest_path()` function of the `NetworkX` package.\n", + "\n", + "The function `ox.nearest_nodes()` looks for the nearest point in your network based on a `X` and `Y` coordinate. For example, in the code below, the origin node is based on the northwestern corner of your bounding box, whereas the destination node is based on the coordinates of the southeastern corner of your bounding box. \n", + "\n", + "So this can also be rewritten as:\n", + "\n", + "```\n", + "origin_node = ox.nearest_nodes(G,4.65465, 56.6778) \n", + "destination_node = ox.nearest_nodes(G,4.61055, 59.5487) \n", + "route = nx.shortest_path(G, origin_node, destination_node) \n", + "```" + ] + }, + { + "cell_type": "markdown", + "id": "6e1d7552-5ec6-42f0-90a7-7e70b5c435ca", + "metadata": { + "id": "Zi_xhqXTCi3h" + }, + "source": [ + "
\n", + "Question 15: The last element of this tutorial is to play around with routing. Please explain in your own words what the .shortest_path() algorithm does. Include the term 'Dijkstra algorithm' in your answer.\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a5e7faf0-43f4-4947-a5b7-0d5f2c1632ea", + "metadata": { + "executionInfo": { + "elapsed": 837, + "status": "ok", + "timestamp": 1675088236066, + "user": { + "displayName": "RA Odongo", + "userId": "17326618845752559881" + }, + "user_tz": -60 + }, + "id": "Dr7Ftbbh47DW" + }, + "outputs": [], + "source": [ + "origin_node = ox.nearest_nodes(G,area['bbox_west'].values[0], area['bbox_north'].values[0])\n", + "destination_node = ox.nearest_nodes(G,area['bbox_east'].values[0], area['bbox_south'].values[0])\n", + "route = nx.shortest_path(G, origin_node, destination_node)" + ] + }, + { + "cell_type": "markdown", + "id": "b3724e1e-0cfa-4e7c-8d94-1d76e10c06d1", + "metadata": { + "id": "4UuCS8G247DW" + }, + "source": [ + "We can plot the route with folium. Like above, you can pass keyword args along to folium PolyLine to style the lines." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ac7785a7-c4c3-491b-8c65-44353cdddc43", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 866 + }, + "executionInfo": { + "elapsed": 240, + "status": "ok", + "timestamp": 1675088397348, + "user": { + "displayName": "RA Odongo", + "userId": "17326618845752559881" + }, + "user_tz": -60 + }, + "id": "7QrZ6Gpi47DX", + "outputId": "33fecead-a6d7-4090-99da-e09d0c34081a" + }, + "outputs": [], + "source": [ + "m2 = ox.plot_route_folium(G, route, weight=10)\n", + "m2" + ] + }, + { + "cell_type": "markdown", + "id": "9df6de86-02e3-4a06-b99f-ec39c71385e8", + "metadata": { + "id": "eedHqPli47DX" + }, + "source": [ + "Plot the route with folium on top of the previously created map\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d0813cf0-0fb3-4ba7-9012-dc5c3f362743", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 866 + }, + "executionInfo": { + "elapsed": 2518, + "status": "ok", + "timestamp": 1675088413924, + "user": { + "displayName": "RA Odongo", + "userId": "17326618845752559881" + }, + "user_tz": -60 + }, + "id": "aHzxZAbZ47DX", + "outputId": "01ee7d97-4792-41fb-925c-c44fb9b15a3f" + }, + "outputs": [], + "source": [ + "m3 = ox.plot_route_folium(G, route, route_map=m1, popup_attribute=\"name\", weight=7)\n", + "m3" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "d0eda49e-6b12-42b9-b88d-ba0938f2f507", + "metadata": {}, + "source": [ + "
\n", + "Question 16: Please add one more routes on a map and upload the resulting figure here.\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "508022f5-46d7-4c69-91d2-46711d6a514c", + "metadata": {}, + "source": [ + "## 7. Perform a damage assessment of the road network using OpenStreetMap\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "b8d7623d-40ee-4262-9781-0c2ed3ea1889", + "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, + "id": "abacf8e1-b27c-4156-8521-87f38004ac2e", + "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", + "id": "87de9a96-e3fc-4ab5-b090-be6412a0d0a1", + "metadata": { + "id": "modIJTEz-g4K" + }, + "source": [ + "Now we convert the road network to a `geodataframe`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "90acac15-3b1e-4c0a-b0a4-6033d243f9ea", + "metadata": {}, + "outputs": [], + "source": [ + "roads = gpd.GeoDataFrame(nx.to_pandas_edgelist(G))\n", + "roads.highway = roads.highway.astype('str')" + ] + }, + { + "cell_type": "markdown", + "id": "56f95cc0-2611-472f-92fa-c256f368c33c", + "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, + "id": "55752d12-81de-4eb5-b849-90473f88478a", + "metadata": {}, + "outputs": [], + "source": [ + "#from urllib import request\n", + "# remote_url = 'https://github.com/VU-IVM/UNIGIS_ProgrammingGIS/tree/main/TAA2'\n", + "# file = 'kampen_roads.gpkg'\n", + "# \n", + "# #request.urlretrieve(remote_url, file)\n", + "# roads = gpd.GeoDataFrame.from_file('kampen_roads.gpkg')" + ] + }, + { + "cell_type": "markdown", + "id": "73f4a201-87b1-4cd8-a206-6b965c9b4107", + "metadata": { + "id": "VIaMGLxA-g4K" + }, + "source": [ + "And lets have a look at the data:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "91a864f8-6af7-454f-869c-197e5e47151b", + "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", + "id": "4789ae59-e6c9-4a1c-a1ad-44e94c9172ca", + "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, + "id": "51c4746a-6efd-4b59-8b6b-5fe9a8c42372", + "metadata": { + "id": "S_LZSRI6-g4L" + }, + "outputs": [], + "source": [ + "roads.highway = roads.highway.apply(lambda x: x.strip('][').split(', '))" + ] + }, + { + "cell_type": "markdown", + "id": "eb311b10-0a1c-4e0a-ad1f-cef37aa40ad7", + "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, + "id": "9e449847-c279-479e-bf3b-382b8c50e018", + "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", + "id": "5b3f95ad-9a1f-43e1-b098-45255c9d5dd2", + "metadata": { + "id": "TkyDDDIP-g4L" + }, + "source": [ + "And let's have a look whether this worked:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21dc3a62-c266-4d0e-bd03-ad1d1b4b22e8", + "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", + "id": "a3049d42-da1f-4dc2-a56f-1ce6b3080fd0", + "metadata": { + "id": "u0crazq8iQjQ" + }, + "source": [ + "
\n", + "Question 17: 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", + "id": "d784f6dc-7091-4e60-b7df-58682efb66ef", + "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, + "id": "af68dc43-fd65-416c-bc7a-a1b1231187d3", + "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", + "id": "9894f1ba-e637-4ede-b46c-5aa4747eaea0", + "metadata": { + "id": "sKb-ig4Q-g4M" + }, + "source": [ + "And let's plot the results:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5f6ab83a-62af-4142-bb26-157e077d868b", + "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", + "id": "1cdbc0d8-beba-4b37-be49-f5226953f19b", + "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, + "id": "d470c747-42aa-4f20-bbf3-90773c3d4cc8", + "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", + "id": "622a243c-dd99-44c1-975c-7cb4eab99359", + "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, + "id": "484c68a5-90ae-4394-9e4f-4ed5ba0c665d", + "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", + "id": "efebb735-7955-4a75-bedf-1b213067fe20", + "metadata": { + "id": "4JT25WTv-g4M" + }, + "source": [ + "And we can now overlay the roads with the flood data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "60de298f-9e06-40c3-b3a2-eecf2aa7205d", + "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", + "id": "53aeaec5-4ea3-48a5-93b5-109a1578c77e", + "metadata": { + "id": "s82DyD_y-g4M" + }, + "source": [ + "And estimate the damages" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "daa56503-90e8-4b7f-a5dc-82a4b235e95b", + "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", + "id": "580663ca-b83f-4724-a699-b6a4790678ad", + "metadata": { + "id": "RFGZxWl7i7pQ" + }, + "source": [ + "
\n", + "Question 18: 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", + "id": "2e42c7f3-c868-400c-858e-7ac0c848fd4c", + "metadata": { + "id": "B5jpsbyC-g4N" + }, + "source": [ + "And let's plot the results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ef910b4e-201e-444b-92b6-826be82a5164", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 514 + }, + "id": "n25j-3wG-g4N", + "outputId": "b926a8e7-7e51-4434-f3b8-d61e6278bc24" + }, "outputs": [], "source": [ - "! @TASK: copy-paste cells from tutorial 2" + "fig, ax = plt.subplots(1, 1,figsize=(12,10))\n", + "\n", + "damaged_roads.plot(column='damage',cmap='Reds',ax=ax);" + ] + }, + { + "cell_type": "markdown", + "id": "49d5033f-ac25-4ae1-9991-ca0899c1b8d4", + "metadata": { + "id": "bTfmvwW2jchB" + }, + "source": [ + "
\n", + "Question 19: 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, - "id": "f1693da4-89a6-449d-bc3c-2ae39fd998d7", + "id": "fe92fb12-9049-4e15-8992-b3d0cbd459b5", "metadata": {}, "outputs": [], "source": [] @@ -1549,7 +3039,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.9" + "version": "3.12.3" } }, "nbformat": 4,