How to use polygons or lines as exposure¶
Exposure in CLIMADA are usually represented as individual points or a raster of points. See Exposures tutorial to learn how to fill and use exposures. In this tutorial we show you how to use CLIMADA Impf you have your exposure in the form of shapes/polygons or in the form of lines.
The approach follows three steps: 1. transform your polygon or line in a set of points 2. do the impact calculation in CLIMADA with that set of points 3. transform the calculated Impact back to your polygon or line
Polygons¶
Polygons or shapes are a common geographical representation of countries, states etc. as for example in NaturalEarth. Here we want to show you how to deal with exposure information as polygons.
Lets assume we have the following data given. The polygons of the admin-1 regions of the netherlands and an exposure value each. We want to know the Impact of Lothar on each admin-1 region.
[1]:
from cartopy.io import shapereader
from climada_petals.entity.exposures.black_marble import country_iso_geom
# open the file containing the Netherlands admin-1 polygons
shp_file = shapereader.natural_earth(resolution='10m',
category='cultural',
name='admin_0_countries')
shp_file = shapereader.Reader(shp_file)
# extract the NL polygons
prov_names = {'Netherlands': ['Groningen', 'Drenthe',
'Overijssel', 'Gelderland',
'Limburg', 'Zeeland',
'Noord-Brabant', 'Zuid-Holland',
'Noord-Holland', 'Friesland',
'Flevoland', 'Utrecht']
}
polygon_Netherlands, polygons_prov_NL = country_iso_geom(prov_names,
shp_file)
# assign a value to each admin-1 area (assumption 100'000 USD per inhabitant)
population_prov_NL = {'Drenthe':493449, 'Flevoland':422202,
'Friesland':649988, 'Gelderland':2084478,
'Groningen':585881, 'Limburg':1118223,
'Noord-Brabant':2562566, 'Noord-Holland':2877909,
'Overijssel':1162215, 'Zuid-Holland':3705625,
'Utrecht':1353596, 'Zeeland':383689}
value_prov_NL = {n: 100000 * population_prov_NL[n] for n in population_prov_NL.keys()}
Assume a unImpform distribution of values within your polygons¶
This helps you in the case you have a given total exposure value per polygon and we assume this value is distributed evenly within the polygon.
We can now perform the three steps for this example:
[2]:
import numpy as np
from pandas import DataFrame
from climada.entity import Exposures
from climada.util.coordinates import coord_on_land
### 1. transform your polygon or line in a set of points
# create exposure with points
exp_df = DataFrame()
n_exp = 200*200
lat, lon = np.mgrid[50 : 54 : complex(0, np.sqrt(n_exp)),
3 : 8 : complex(0, np.sqrt(n_exp))]
exp_df['latitude'] = lat.flatten() # provide latitude
exp_df['longitude'] = lon.flatten() # provide longitude
exp_df['impf_WS'] = np.ones(n_exp, int) # provide impact functions
# now we assign each point a province and a value, Impf the points are within one of the polygons defined above
exp_df['province'] = ''
exp_df['value'] = np.ones((exp_df.shape[0],))*np.nan
for prov_name_i, prob_polygon_i in zip(prov_names['Netherlands'],polygons_prov_NL['NLD']):
in_geom = coord_on_land(lat=exp_df['latitude'],
lon=exp_df['longitude'],
land_geom=prob_polygon_i)
np.put(exp_df['province'].values, np.where(in_geom)[0], prov_name_i)
np.put(exp_df['value'].values, np.where(in_geom)[0], value_prov_NL[prov_name_i]/sum(in_geom))
exp_df = Exposures(exp_df)
exp_df.set_geometry_points() # set geometry attribute (shapely Points) from GeoDataFrame from latitude and longitude
exp_df.check()
exp_df.plot_hexbin()
[2]:
<GeoAxesSubplot:>
[3]:
from climada.hazard.storm_europe import StormEurope
from climada.util.constants import WS_DEMO_NC
from climada.entity.impact_funcs.storm_europe import ImpfStormEurope
from climada.entity.impact_funcs import ImpactFuncSet
from climada.engine import Impact
### 2. do the impact calculation in CLIMADA with that set of points
# define hazard
storms = StormEurope.from_footprints(WS_DEMO_NC, description='test_description')
# define impact function
impact_func = ImpfStormEurope.from_welker()
impact_function_set = ImpactFuncSet()
impact_function_set.append(impact_func)
# calculate hazard
impact_NL = Impact()
impact_NL.calc(exp_df, impact_function_set, storms, save_mat=True)
impact_NL.plot_hexbin_impact_exposure()
$CLIMADA_SRC/climada/hazard/centroids/centr.py:822: UserWarning: Geometry is in a geographic CRS. Results from 'buffer' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.
xy_pixels = self.geometry.buffer(res / 2).envelope
[3]:
<GeoAxesSubplot:>
[4]:
import pandas as pd
### 3. transform the calculated Impact back to your polygon or line
impact_at_province_raw = pd.DataFrame(np.mean(impact_NL.imp_mat.todense().transpose(), axis=1),
index=exp_df.gdf['province'])
impact_at_province = impact_at_province_raw.groupby(impact_at_province_raw.index).sum()
print(impact_at_province)
0
province
0.000000e+00
Drenthe 6.852537e+04
Flevoland 2.716078e+05
Friesland 7.782136e+05
Gelderland 4.056456e+05
Groningen 1.150089e+05
Limburg 2.559739e+05
Noord-Brabant 7.391625e+05
Noord-Holland 5.908293e+06
Overijssel 1.588620e+05
Utrecht 4.846288e+05
Zeeland 4.069767e+05
Zuid-Holland 2.893966e+06
Use the LitPop module to disaggregate your exposure values¶
Instead of a unImpform distribution, another geographical distribution can be chosen to disaggregate within the value within each polygon. We show here the example of LitPop. The same three steps apply:
[5]:
import numpy as np
from climada.entity import LitPop
from climada.util.coordinates import coord_on_land
### 1. transform your polygon or line in a set of points
# create exposure with points
exp_df_lp = LitPop.from_countries('Netherlands',res_arcsec = 60, fin_mode = 'none')
exp_df_lp.gdf['impf_WS'] = np.ones(exp_df_lp.gdf.shape[0], int) # provide impact functions
# now we assign each point a province and a value, Impf the points are within one of the polygons defined above
exp_df_lp.gdf['province'] = ''
for prov_name_i, prob_polygon_i in zip(prov_names['Netherlands'],polygons_prov_NL['NLD']):
in_geom = coord_on_land(lat=exp_df_lp.gdf['latitude'],
lon=exp_df_lp.gdf['longitude'],
land_geom=prob_polygon_i)
np.put(exp_df_lp.gdf['province'].values, np.where(in_geom)[0], prov_name_i)
exp_df_lp.gdf['value'][np.where(in_geom)[0]] = \
exp_df_lp.gdf['value'][np.where(in_geom)[0]] * value_prov_NL[prov_name_i]/sum(exp_df_lp.gdf['value'][np.where(in_geom)[0]])
exp_df_lp.gdf = exp_df_lp.gdf.drop(np.where(exp_df_lp.gdf['province']=='')[0]) #drop carribean islands for this example
exp_df_lp.set_geometry_points()
exp_df_lp.check()
exp_df_lp.plot_hexbin()
2021-10-19 16:54:04,710 - climada.entity.exposures.litpop.gpw_population - WARNING - Reference year: 2018. Using nearest available year for GPW data: 2020
[5]:
<GeoAxesSubplot:title={'center':"LitPop Exposure for ['Netherlands'] at 60 as, year: 2018, financial\nmode: none, exp: (1, 1), admin1_calc: False"}>
[6]:
from climada.hazard.storm_europe import StormEurope
from climada.util.constants import WS_DEMO_NC
from climada.entity.impact_funcs.storm_europe import ImpfStormEurope
from climada.entity.impact_funcs import ImpactFuncSet
from climada.engine import Impact
### 2. do the impact calculation in CLIMADA with that set of points
# define hazard
storms = StormEurope.from_footprints(WS_DEMO_NC, description='test_description')
# define impact function
impact_func = ImpfStormEurope.from_welker()
impact_function_set = ImpactFuncSet()
impact_function_set.append(impact_func)
# calculate hazard
impact_NL = Impact()
impact_NL.calc(exp_df_lp, impact_function_set, storms, save_mat=True)
impact_NL.plot_hexbin_impact_exposure()
$CLIMADA_SRC/climada/hazard/centroids/centr.py:822: UserWarning: Geometry is in a geographic CRS. Results from 'buffer' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.
xy_pixels = self.geometry.buffer(res / 2).envelope
[6]:
<GeoAxesSubplot:>
[7]:
import pandas as pd
### 3. transform the calculated Impact back to your polygon or line
impact_at_province_raw = pd.DataFrame(np.mean(impact_NL.imp_mat.todense().transpose(), axis=1),
index=exp_df_lp.gdf['province'])
impact_at_province_lp = impact_at_province_raw.groupby(impact_at_province_raw.index).sum()
print(impact_at_province_lp)
0
province
Drenthe 6.500784e+04
Flevoland 1.684639e+05
Friesland 4.575450e+05
Gelderland 3.345681e+05
Groningen 1.604080e+05
Limburg 3.647827e+05
Noord-Brabant 6.253040e+05
Noord-Holland 1.819285e+06
Overijssel 1.094167e+05
Utrecht 4.495409e+05
Zeeland 7.701477e+05
Zuid-Holland 2.823286e+06
Comparison of both modelled impacts:
[8]:
from matplotlib import pyplot as plt
plt.plot(impact_at_province[impact_at_province.index!=''],impact_at_province_lp, '.k')
plt.xlabel('unImpform disaggregation')
plt.ylabel('litpop disaggregation')
[8]:
Text(0, 0.5, 'litpop disaggregation')
Further statistical analysis of hazard on polygon level¶
imagine that you need access to the hazard centroids in oder to provide some statistical analysis on the province level
[9]:
# this provides the wind speed value for each event at the corresponding exposure
import scipy
exp_df_lp.gdf[:5]
l1, l2, vals = scipy.sparse.find(storms.intensity)
# provide columns for both events
exp_df_lp.gdf['wind_0'] = 0
exp_df_lp.gdf['wind_1'] = 0
for evt, idx, val in zip(l1, l2, vals):
if evt == 0:
exp_df_lp.gdf['wind_0'].values[exp_df_lp.gdf['centr_WS'] == idx] = val
else:
exp_df_lp.gdf['wind_1'].values[exp_df_lp.gdf['centr_WS'] == idx] = val
[10]:
# now you can perform additional statistical analysis and aggregate it to the province level
import pandas as pd
import geopandas as gpd
exp_province_raw = exp_df_lp.copy()
def f(x): # define function for statistical aggregation with pandas
d = {}
d['value'] = x['value'].sum()
d['wind_0'] = x['wind_0'].max()
d['wind_1'] = x['wind_1'].mean()
# one could also be interested in centroid of max wind with respect to province
#d['centr_WS'] = x.loc[x.index[x['wind_0'].max()],'centr_WS']
return pd.Series(d, index=['value', 'wind_0', 'wind_1'])
exp_province = exp_province_raw.gdf.groupby('province').apply(f).reset_index() # Result is not a GeoDataFrame anymore
# add geometries to DataFrame and plot results
exp_province=gpd.GeoDataFrame(exp_province, geometry=None)
for prov,poly in zip(list(prov_names.values())[0],polygons_prov_NL['NLD']):
exp_province.loc[exp_province.index[exp_province['province']== prov],'geometry']= gpd.GeoDataFrame(geometry=[poly]).geometry.values
exp_province
print('Plot maximum wind per province for first event')
exp_province.plot(column='wind_0', cmap='Reds', legend=True)
plt.show()
print('Plot mean wind per province for second event')
exp_province.plot(column='wind_1', cmap='Reds', legend=True)
plt.show()
Plot maximum wind per province for first event
Plot mean wind per province for second event
Lines¶
Lines are common geographical representation of transport infrastructure like streets, train tracks or powerlines etc.
[11]:
# under construction. here follows an example on how to deal with lines