From f4f20da035b6c7b47ef8fb4d22be6494c6411a97 Mon Sep 17 00:00:00 2001 From: emanuel-schmid Date: Tue, 24 Oct 2023 21:10:34 +0200 Subject: [PATCH] folium plot workshop --- .../Plotting Climada Data with Folium.ipynb | 273 ++++++++++++++++++ climada-days/2023-10/climada_folium.py | 227 +++++++++++++++ 2 files changed, 500 insertions(+) create mode 100644 climada-days/2023-10/Plotting Climada Data with Folium.ipynb create mode 100644 climada-days/2023-10/climada_folium.py diff --git a/climada-days/2023-10/Plotting Climada Data with Folium.ipynb b/climada-days/2023-10/Plotting Climada Data with Folium.ipynb new file mode 100644 index 000000000..f7513f0dc --- /dev/null +++ b/climada-days/2023-10/Plotting Climada Data with Folium.ipynb @@ -0,0 +1,273 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "b650547e-4c47-4df0-a3a7-047fefeceacc", + "metadata": {}, + "source": [ + "# Interactive maps with folium\n", + "\n", + "- home: [Folium](https://python-visualization.github.io/folium/latest/index.html)\n", + "- based on [Leaflet](https://leafletjs.com/)\n", + "- pretty cool library for creating interactive html/java-script maps from python" + ] + }, + { + "cell_type": "markdown", + "id": "1a44b15a-e7e1-4b5c-9b85-d36ca5d38627", + "metadata": {}, + "source": [ + "## Installation\n", + "\n", + "```bash\n", + "pip install folium\n", + "```\n", + "\n", + "## Basics" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16ca28c3-c6bb-4f53-a2a0-893369ff040c", + "metadata": {}, + "outputs": [], + "source": [ + "import folium\n", + "[x for x in dir(folium) if not x.startswith('_')]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c45a5d3b-2599-40af-aa3b-e22bd341477b", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# A map of the London Area (initially)\n", + "latlon = (51.5, 0.0)\n", + "m = folium.Map(location=latlon)\n", + "\n", + "# show coordingates\n", + "#folium.LatLngPopup().add_to(m)\n", + "\n", + "# a Marker\n", + "#folium.Marker(location=(51.3, -0.4), popup='just an ordinary marker').add_to(m)\n", + "\n", + "# a Circle with size in meters\n", + "#folium.Circle(location=(51.7, -0.4), radius=2000, color='red').add_to(m)\n", + "\n", + "# also a Circle but size in pixels\n", + "#folium.CircleMarker(location=(51.7, +0.4), radius=20, color='green').add_to(m)\n", + "\n", + "# a Polygon\n", + "#folium.Polygon(locations=[(51.31, 0.42),(51.31, 0.39),(51.29, 0.38),(51.29, 0.41)], fillColor='#3186cc', fillOpacity=0.4).add_to(m)\n", + "\n", + "# change the Tiles\n", + "#folium.TileLayer('stamenterrain').add_to(m)\n", + "\n", + "# switch between Tiles\n", + "#folium.LayerControl().add_to(m)\n", + "\n", + "m" + ] + }, + { + "cell_type": "markdown", + "id": "e0fef59f-1808-43a6-a616-d592a3a24957", + "metadata": {}, + "source": [ + "## Plotting Climada Data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ee4af038-dc2d-4433-9853-943cad93a888", + "metadata": {}, + "outputs": [], + "source": [ + "from climada.util.api_client import Client\n", + "c = Client()\n", + "\n", + "uk_lp_30 = c.get_litpop(country='GBR', exponents=(3,0))\n", + "uk_lp_01 = c.get_litpop(country='GBR', exponents=(0,1))\n", + "uk_lp_11 = c.get_litpop(country='GBR', exponents=(1,1))\n", + "\n", + "uk_lp_11.gdf.head()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "90923796-76ca-4389-990a-a25704d2279c", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "folium.GeoJson?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d427d2ae-25e0-411d-ae8e-4a86181fa271", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "latlon = (51.5, 0.0)\n", + "m = folium.Map(location=latlon)\n", + "maxv = uk_lp_30.gdf.value.max()\n", + "folium.GeoJson(\n", + " data=uk_lp_30.gdf[\n", + " uk_lp_30.gdf.latitude.between(51.3, 51.7)\n", + " & uk_lp_30.gdf.longitude.between(-0.4, 0.4)\n", + " ].loc[:,['geometry', 'value']],\n", + " #marker=folium.Circle(radius=1500, fillColor='blue'),\n", + " #style_function=lambda lp: {\n", + " # 'opacity': 0,\n", + " # 'fillColor': 'blue',\n", + " # 'fillOpacity': 0.6 * lp['properties']['value']/maxv,\n", + " #},\n", + " #popup=folium.GeoJsonPopup(fields=['value']),\n", + ").add_to(m)\n", + "m" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a2f18d1c-c304-4094-925e-f35d529da3c7", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "from shapely import Polygon\n", + "latlon = (51.5, 0.0)\n", + "m = folium.Map(location=latlon)\n", + "maxv = uk_lp_30.gdf.value.max()\n", + "\n", + "def raster_polygon(lat, lon, arc_res):\n", + " return Polygon([\n", + " (lon+arc_res/2, lat+arc_res/2),\n", + " (lon+arc_res/2, lat-arc_res/2),\n", + " (lon-arc_res/2, lat-arc_res/2),\n", + " (lon-arc_res/2, lat+arc_res/2)])\n", + " \n", + "data = uk_lp_30.gdf[\n", + " uk_lp_30.gdf.latitude.between(51.3, 51.7)\n", + " & uk_lp_30.gdf.longitude.between(-0.4, 0.4)\n", + " ]\n", + "\n", + "data['geometry'] = data.apply(lambda x: raster_polygon(x.latitude, x.longitude, 150/3600), axis=1)\n", + "folium.GeoJson(\n", + " data=data,\n", + " popup=folium.GeoJsonPopup(fields=['value']),\n", + " style_function=lambda lp: {\n", + " 'opacity': 0,\n", + " 'fillColor': 'blue',\n", + " 'fillOpacity':0.6 * lp['properties']['value']/maxv,\n", + " }\n", + ").add_to(m)\n", + "m" + ] + }, + { + "cell_type": "markdown", + "id": "5a51710f-42ae-46ca-8b25-3b698ea4f6b4", + "metadata": {}, + "source": [ + "## Annotating large areas\n", + "\n", + "`folium.Choropleth` is better suited for large areas then `folium.GeoJson` and provides a sophisticated interface for plotting choropleth maps." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e4d32459-9ea3-4d68-9609-c9d26f0fe3a5", + "metadata": {}, + "outputs": [], + "source": [ + "folium.Choropleth?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "caac0f3b-f975-4a5b-acf2-994c20889f50", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "from climada_folium import exposures_folium, raster_polygons\n", + "\n", + "latlon = (51.5, 0.0)\n", + "m = folium.Map(location=latlon)\n", + "\n", + "exposures_folium(\n", + " exposures=uk_lp_30,\n", + " polygon_maker=raster_polygons(150/3600),\n", + " add_to_map=m\n", + ")\n", + "m" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "41a1e692-f27e-4eaf-b120-11ed3591db82", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "exposures_folium(\n", + " exposures=uk_lp_01,\n", + " polygon_maker=raster_polygons(150/3600),\n", + " fill_color='YlGn',\n", + " add_to_map=m\n", + ")\n", + "\n", + "exposures_folium(\n", + " exposures=uk_lp_11,\n", + " polygon_maker=raster_polygons(150/3600),\n", + " fill_color='OrRd',\n", + " add_to_map=m\n", + ")\n", + "\n", + "folium.LayerControl().add_to(m)\n", + "\n", + "m" + ] + } + ], + "metadata": { + "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.9.18" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/climada-days/2023-10/climada_folium.py b/climada-days/2023-10/climada_folium.py new file mode 100644 index 000000000..ca806c8b0 --- /dev/null +++ b/climada-days/2023-10/climada_folium.py @@ -0,0 +1,227 @@ +import math + +import folium +import geopandas as gpd +import numpy as np +import pandas as pd +from shapely import Polygon + +import climada.util.coordinates as u_coord + + +def estimate_zoom_level(geodf): + """crude implementation for guessing the initial zoom level + based on the latitudinal expansion of the geo dataframe + TODO: find a better formula + """ + vdist = abs(geodf.latitude.max() - geodf.latitude.min()) + return 9 - int(math.log2(vdist) / 0.8) + + +def center(geodf): + cy = (geodf.latitude.max() + geodf.latitude.min()) / 2 + cx = (geodf.longitude.max() + geodf.longitude.min()) / 2 + return (cy, cx) + + +def voronoi_polygons(country): + from geovoronoi import voronoi_regions_from_coords + + iso3a = u_coord.country_to_iso(country, representation="alpha3") + country_shape = u_coord.get_land_geometry([iso3a]) + def _(geodf): + coords = geodf.loc[:,['longitude', 'latitude']].values + voronoi_polys, voronoi_idx = voronoi_regions_from_coords(coords, country_shape) + remap = {vi: k for k,v in voronoi_idx.items() for vi in v} + return pd.Series([voronoi_polys[remap.get(i, 0) # TODO: gross! fix it! + ] for i in range(geodf.shape[0])]) + return _ + + +def raster_polygons(width): + def _(geodf): + return geodf.apply(lambda row: Polygon(( + (row.longitude - width/2, row.latitude - width/2), + (row.longitude - width/2, row.latitude + width/2), + (row.longitude + width/2, row.latitude + width/2), + (row.longitude + width/2, row.latitude - width/2), + (row.longitude - width/2, row.latitude - width/2), + )), axis=1) + return _ + + +def make_bins(dseries, n): + distance = dseries.max() - dseries.min() + om = int(math.log(distance, 10)) + lower = int(dseries.min() / math.pow(10, om)) + upper = (int(dseries.max() / math.pow(10, om)) + 1) + step = (upper - lower) / n + return [lower + n*step for n in range(n+1)], om + + +def show_folium(mymap, + geodf, + polygon_maker, + description, + value_unit, + val_col, + fill_color, + threshold=0, + bins_no=5, + crs='EPSG:4326'): + + polygons = polygon_maker(geodf) + mygdf = gpd.GeoDataFrame({ + 'geometry': polygons, + val_col: geodf[val_col], + 'latitude': geodf.latitude, + 'longitude': geodf.longitude + }, geometry='geometry', crs=crs) + bins, om = make_bins(mygdf[val_col], bins_no) + + mygdf['evalue'] = mygdf[val_col] / math.pow(10, om) + + choropleth = folium.Choropleth( + geo_data=mygdf[mygdf[val_col]>threshold], + data=mygdf.evalue, + key_on="feature.id", # index of the geodataframe is transformed into `id` field + fill_color=fill_color, + fill_opacity=0.6, + bins=bins, + line_weight=0, + line_color="black", + nan_fill_opacity=0.0, + legend_name=f"{val_col} [1.0e{om:+03d} {value_unit}]", # title under the color scale + name=f"{description}", # name of thew layer, e.g. in the layer control + ).add_to(mymap) + + cols_xy = [c for c in ["latitude", "longitude"] if c in mygdf] + labels_xy = ["latitude:", "longitude:"] if cols_xy else [] + # add tooltip to appear, when pointing at a hectar + tooltip = folium.GeoJsonTooltip( + # column names with values to be displayed + fields=["evalue"] + cols_xy, + # text to be shown explaining each value + aliases=[f"{val_col} [1.0e{om:+03d}]"] + labels_xy, + localize=True, + max_width=800, + ) + tooltip.add_to(choropleth.geojson) + + return mymap + + +def map_maker(geodf, zoom_start=None): + (sy,sx) = center(geodf) + zoom_start = zoom_start or estimate_zoom_level(geodf) + mymap = folium.Map(location=(sy,sx), zoom_start=zoom_start) + return mymap + + +def exposures_folium( + exposures, + polygon_maker, + add_to_map=None, + val_col='value', + threshold=0, + bins_no=5, + fill_color="PuBu", + zoom_start=None): + """Creates a folium (leaflet.js) map from the exposures object + + Parameters + ---------- + exposures : climada.entity.Exposures + polygon_maker : function(geopandas.GeoDataFrame) -> geopandas.GeoDataFrame + applied to exposures.gdf + supposed to make a copy and add a geometry column with Polygons around lat/lon + add_to_map : folium.Map, optional + the map where this layer is being integrated, + if None a new map is created based on the exposures' coordinates + val_col : str, optional + the column name of the exposures value of interest, default 'value' + threshold : int, optional + below this level values are not showed in the map, default 0 + bins_no : int, optional + number of categories/bins to be shown, default 5 + fill_color : str, optional + must be a ColorBrewer's name (see https://colorbrewer2.org) + default: 'PuBu', Purple -> Blue + zoom_start : int, optional + starting zoom level, if None, the zoom level is guessed from the latitudinal range + + Returns + ------- + folium.folium.Map + A folium Map with the exposures values shading a polygon around their lat/lon coordinates + """ + return show_folium(mymap=add_to_map or map_maker(exposures.gdf, zoom_start=zoom_start), + geodf=exposures.gdf, + polygon_maker=polygon_maker, + description=exposures.description or val_col, + value_unit=exposures.value_unit, + val_col=val_col, + fill_color=fill_color, + threshold=threshold, + bins_no=bins_no) + + +def impact_folium( + impact, + polygon_maker, + add_to_map=None, + label=None, + threshold=0, + bins_no=5, + fill_color="YlOrBr", + zoom_start=None): + """Creates a folium (leaflet.js) map from the exposures object + + Parameters + ---------- + impact : climada.engine.Impact + polygon_maker : function(geopandas.GeoDataFrame) -> geopandas.GeoDataFrame + applied to exposures.gdf + supposed to make a copy and add a geometry column with Polygons around lat/lon + add_to_map : folium.Map, optional + the map where this layer is being integrated, + if None a new map is created based on the exposures' coordinates + label : str, optional + impact label, will be shown in the LayerControl of the folium.Map + default is 'Expected annual impact for hazard type [HAZ_TYPE]' + threshold : int, optional + below this level values are not showed in the map, default 0 + bins_no : int, optional + number of categories/bins to be shown, default 5 + fill_color : str, optional + must be a ColorBrewer's name (see https://colorbrewer2.org) + default 'YlOrBr', Yellow -> Orange -> Brown + zoom_start : int, optional + starting zoom level, if None, the zoom level is guessed from the latitudinal range + + Returns + ------- + folium.folium.Map + A folium Map with the exposures values shading a polygon around their lat/lon coordinates + """ + imp_df = pd.DataFrame({ + 'longitude': impact.coord_exp[:,1], + 'latitude': impact.coord_exp[:,0], + 'eai_exp': impact.eai_exp}) + label = label or f"Expected annual impact for hazard type {impact.haz_type}" + return show_folium(mymap=add_to_map or map_maker(imp_df, zoom_start=zoom_start), + geodf=imp_df, + polygon_maker=polygon_maker, + description=label, + value_unit=impact.unit, + val_col='eai_exp', + fill_color=fill_color, + threshold=threshold, + bins_no=bins_no) + + +def add_layer_control(mymap): + # https://github.com/python-visualization/folium/issues/816 + # "Both folium and Leaflet require the LayerControl to be added last: otherwise it doesn't know what layers there are." + folium.LayerControl().add_to(mymap) + return mymap