This notebook illustrates two blogposts (in French) :
Requires same packages as in "Spatial_Join" notebook, plus :
- conda install -c conda-forge dash
- conda install -c conda-forge -c plotly jupyter-dash
- conda install descartes
- conda install openpyxl
- conda install -c conda-forge flexpolyline
- conda install -c conda-forge python-kaleido
%matplotlib inline
import urllib
import time
import pandas as pd
import numpy as np
import geopandas as gpd
from shapely.geometry import Polygon
import contextily as ctx
import shapely
from shapely.geometry import shape
from pyproj import CRS
import os
import requests
import matplotlib.pyplot as plt
from tqdm.auto import tqdm, trange
import json
tqdm.pandas()
import shapely
import plotly.express as px
import dash
from jupyter_dash import JupyterDash
import dash_core_components as dcc
import dash_html_components as html
from dash.dependencies import Input, Output
import plotly.graph_objects as go
import flexpolyline as fp
from credentials import here_api_key, here_app_code, here_app_id
crs = 'epsg:3857'
osm_crs= 'epsg:4326'
try : os.mkdir("output")
except FileExistsError: pass
from_files = True
# If False, Here API is called for geocoding and isolines, and result is stored in a file
# If True, API is not called, data is retreived from the file and plots are produced
add_options={}
option=None
dataset= "smur"
# dataset = "maternity"
# dataset = "vaccin"
# Uncomment to two following lines to add traffic
# add_options = {"departureTime" : "2021-05-28T14:00:00+01:00"} # with heavy traffic
# option="traffic"
# # Could also be {"transportMode":"pedestrian"}
def get_geocode(addr):
"""
Geocode an address with Here API
Parameters
----------
addr: str
address to geocode
filename: str
local file to save
Returns
-------
(lat, lon): (float, float)
geographical coordinates (in epsg:4326)
"""
host_name = "geocoder.api.here.com"
params = urllib.parse.urlencode({"app_id": here_app_id,
"app_code": here_app_code,
"searchtext":addr,
})
url = f"https://{host_name}/6.2/geocode.json?{params}"
# print(url)
with urllib.request.urlopen(url) as response:
res = response.read()
res = json.loads(res)
r = res["Response"]["View"][0]["Result"][0]["Location"]["DisplayPosition"]
return r["Latitude"], r["Longitude"]
def geocode_points(starts):
"""
Geocode all addresses in a geopandas GeoDataFrame. Add columns "lat", "lon" (in epsg:4326),
as well as "geometry" (in epsg:3857)
Parameters
----------
starts: geopandas.GeoDataFrame
a dataframe which contains a column "address"
Returns
-------
geopandas.GeoDataFrame, with 3 additionnal columns : "lat", "lon", and "geometry"
"""
starts[["lat", "lon"]] = starts.address.progress_apply(get_geocode).apply(pd.Series)
# Convert lat/long in xy points
starts= starts.set_geometry(gpd.points_from_xy(x = starts["lon"], y=starts["lat"]))
# Set CRS
starts.crs = CRS(osm_crs)
starts = starts.to_crs(CRS(crs))
return starts
def get_isoline(lat, lon, range_list, add_options={}):
"""
Get isolines 'lat, lon' as origin points, for all delays given in 'range_list' in minutes.
Call Here isoline API, with some default parameters that can be overwritten or enriched with 'add_options'
If 'Too Many Requests' error is received from Here, wait for 15 seconds and retries
Parameters
----------
lat: float
Latitude (in epsg:4326)
lon: float
Longitude (in epsg:4326)
range_list: list
List of range values (in minutes)
add_options : dict
Example: {"departureTime" : "2021-05-28T14:00:00+01:00", "transportMode":"pedestrian"}
Returns
-------
A dict, with keys from range_list, and data is a list of coordinates tuple.
Example :
{10.0: [(50.818634, 4.257374),
(50.818291, 4.258232),
(50.817604, 4.258919),
(...)],
15.0: (50.821381, 4.266815),
(50.821037, 4.267845),
(50.820351, 4.268532),
(...)]}
"""
host_name = "isoline.router.hereapi.com"
params = {"apiKey": here_api_key,
"transportMode":f"car",
"departureTime":"any", # No traffic dependent
"origin" : f"{lat},{lon}",
"range[values]" :",".join([str(rng*60) for rng in range_list]),
"range[type]":"time"# "distance"
}
for k in add_options:
params[k] = add_options[k]
params = urllib.parse.urlencode(params)
url = f"https://{host_name}/v8/isolines?{params}"
# print(url)
try:
with urllib.request.urlopen(url) as response:
res = response.read()
res = json.loads(res)
# print(res)
# return res
for it in res["isolines"]:
if len(it["polygons"]) != 1:
print("More than 1 polygon!!")
return {it["range"]["value"]/60: fp.decode(it["polygons"][0]["outer"]) for it in res["isolines"]} #["isoline"][0]["component"][0]["shape"]
except urllib.error.HTTPError as ex:
if ex.getcode() == 429: #Too Many Requests
print("Too many requests, wait for 15 seconds...")
time.sleep(15)
return get_isoline(lat, lon, range_list)
else:
print(ex)
print(type(ex))
print(url)
return None
# get_isoline(50.83458, 4.33639, [10])
def get_isolines(starts, range_list, add_options={}):
"""
From the output of geocode_points(starts), call get_isoline for each record, with range_list and add_options as parameter.
Note: does NOT support multi-component isolines (https://developer.here.com/documentation/isoline-routing-api/8.4.0/dev_guide/topics/multi-component-isoline.html)
Parameters
----------
start: geopandas.GeoDataFrame
output of geocode_points(starts)
range_list: list
list of range values (in minutes)
add_options : dict
Example: {"departureTime" : "2021-05-28T14:00:00+01:00", "transportMode":"pedestrian"}
Returns
-------
geopandas.GeoDataFrame
A geopandas.GeoDataFrame, with starts.shape[0] * len(range_list) record. For each starting point and each range,
"geometry" columns is a polygon with the corresponding isoline
"""
isolines= starts.copy()
isolines["isolines"] = isolines.progress_apply(lambda row: get_isoline(row["lat"], row["lon"], range_list, add_options), axis=1)
for rng in range_list:
isolines[rng] = isolines["isolines"].apply(lambda x: x[rng])
isolines = isolines.drop(["geometry", "isolines"], axis=1)
cols= isolines.drop(range_list, axis=1).columns
isolines = isolines.set_index(list(cols)).stack().rename("isoline").reset_index().rename(columns={"level_7": "range",
"level_8": "range",
"level_9": "range",
})
# Convert isolines to polygons
isolines["geometry"] = isolines["isoline"].apply(lambda row:Polygon(map(lambda tpl: (tpl[1], tpl[0]), row)))
isolines = isolines.set_geometry("geometry")
isolines.crs = CRS(osm_crs)
isolines = isolines.to_crs(CRS(crs))
return isolines.drop("isoline", axis=1)
def get_isolines_grouped(isolines, subsets):
"""
From the output of get_isolines(starts), group all isolines for the same range, making a (multi)polygon from all the
(single starting points) polygons
Parameters
----------
isolines: geopandas.GeoDataFrame
output of get_isolines
subsets: list
list of scenarii (corrsponding to boolean columns in isolines dataframe)
Returns
-------
geopandas.GeoDataFrame
A geopandas.GeoDataFrame, with len(subsets) * len(range_list) record. For each scenario point and each range,
"geometry" columns is a (mutli)polygon with the corresponding reachable region
"""
isolines_grouped = []
for subset in subsets:
gr = isolines[isolines[subset]].groupby("range").apply(lambda bloc: bloc.unary_union).rename("geometry").reset_index()
gr["subset"]=subset
isolines_grouped.append(gr)
isolines_grouped = gpd.GeoDataFrame(pd.concat(isolines_grouped))
isolines_grouped.crs = CRS(crs)
isolines_grouped["label"] = isolines_grouped.apply(lambda x: f"{x['subset']} ({x['range']})", axis=1)
return isolines_grouped
# get_isolines_grouped(isolines, subsets)
def add_basemap(ax, alpha=None):
"""
Add a basemap on a plot. Tries first default (Stamen) basemap. If errors, tries OpenStreetMap.Mapnik
Parameters
----------
ax: matplotlib axes
Returns
-------
None
"""
try:
ctx.add_basemap(ax, alpha=alpha)
except requests.HTTPError:
print("Default basemap doesn't work...")
ctx.add_basemap(ax, source=ctx.providers.OpenStreetMap.Mapnik, alpha=alpha)
def disjoin_labels(geo_df):
"""
From the output of get_isolines_grouped, filtred on a unique range (we should then have len(subsets) records,
extract the overlapping regions.
Example :
geo_df= [["subs1", geom1],
["subs2", geom2]]
if geom2 in full included in geom1, output will be:
[["subs1", geom1-geom2],
["all scenarii", geom2]]
if there is an ovelapping region "common":
[["subs1", geom1-common],
["subs2", geom2-common],
["all scenarii", common]]
If more than 2 subsets, extracts also overlappings two by two
Parameters
----------
geo_df: geopandas.GeoDataFrame
Returns
-------
None
"""
## Find global common zone (common to all zones)
common_zone = geo_df.iloc[[0]]
for i in range(1, geo_df.shape[0]):
common_zone = gpd.overlay(common_zone,
geo_df.iloc[[i]],
how="intersection",
keep_geom_type=False)
common_zone = common_zone.geometry
if common_zone.shape[0]>0:
common_zone_s = gpd.GeoSeries([common_zone.iloc[0] for it in geo_df.iterrows()])
common_zone_s.crs= geo_df.crs
geo_df = geo_df.reset_index(drop=True).copy()
geo_df["common_zone"] = common_zone_s
try:
geo_df["geometry"]= geo_df.geometry.difference(geo_df["common_zone"])#.plot()
except shapely.errors.TopologicalError as e:
print(e)
print("Trying with buffer(0)...")
geo_df["geometry"]= geo_df.geometry.buffer(0).difference(geo_df["common_zone"])
#return geo_df#.drop("common_zone", axis=1)
m_range = geo_df.range.iloc[0]
geo_df = pd.concat([geo_df, gpd.GeoDataFrame([{"subset": "all scenarii",
"label" : f"all scenarii ({m_range})",
"geometry":common_zone.iloc[0],
"range": m_range,
"reach_ratio":0 # not correct!
}])]).drop("common_zone", axis=1).copy()
## Find one2one (remaining) overlaps (only if len(subsets)>2)
geo_df = gpd.GeoDataFrame(geo_df)
sh = geo_df.shape[0]
for i in range(sh):
for j in range(i+1, sh):
common_zone = geo_df.iloc[i].geometry.intersection(geo_df.iloc[j].geometry)
if geo_df.iloc[i].geometry.area > 0 and common_zone.area/geo_df.area.max() > 0.01:
geo_df["geometry"].iloc[i] = geo_df.iloc[i].geometry.difference(common_zone)
geo_df["geometry"].iloc[j] = geo_df.iloc[j].geometry.difference(common_zone)
geo_df = pd.concat([geo_df, gpd.GeoDataFrame([{"subset": f"{geo_df.iloc[i].subset} + {geo_df.iloc[j].subset}",
"label" : f"{geo_df.iloc[i].subset} + {geo_df.iloc[j].subset} ({m_range})",
"geometry":common_zone,
"range": m_range,
"reach_ratio":0 # not correct!
}])]).copy()
#common_zone.plot()
geo_df = geo_df[~geo_df.geometry.is_empty & (geo_df.area/geo_df.area.max() > 0.01)]
return geo_df
# from matplotlib import cm
from matplotlib.colors import ListedColormap
def plot_isolines(isolines_grouped, starts, m_range, boundary):
cmap_st = ListedColormap( [f"C{i}" for i in range(starts.label.nunique())])
ax=starts.sort_values("label").plot("label", cmap=cmap_st, figsize=(10, 10), legend=(m_range is None))
boundaries.boundary.plot(ax=ax, color="black")
if m_range is not None:
geo_df = isolines_grouped[(isolines_grouped.range == m_range)]
geo_df = disjoin_labels(geo_df).sort_values("label")
cmap = ListedColormap( [f"C{i}" for i in range(geo_df.label.nunique())])
ax= geo_df.plot("label", ax=ax, alpha=0.2, cmap=cmap)
# ax = geo_df.plot("label", ax=ax, facecolor="none", edgecolor="k", cmap=cmap, legend=True)
ax = geo_df.plot("label", ax=ax, facecolor="none", cmap=cmap, legend=True)
# boundaries.boundary.plot(ax=ax, color="C2", linewidth=5)
minx, miny, maxx, maxy = boundaries.total_bounds
width = maxx-minx
height = maxy-miny
ax.set_xlim(minx - width/10, maxx + width/10)
ax.set_ylim(miny - height/10, maxy + height/10)
ax.set_title(f"Reached zone ({dataset}{', '+option if option else ''})" + ("" if m_range is None else f" - range: {m_range} minutes"))
add_basemap(ax)#, alpha=0.8)
return ax
# ax= isolines[isolines.range == minute_range[5]].plot("name", alpha=0.8, figsize=(10, 10), legend=True)
# boundaries.boundary.plot(ax=ax, color="C2", linewidth=5)
# add_basemap(ax)
# ax = plot_isolines(isolines_grouped, starts, minute_range[10], boundaries)
# plt.axis("off")
# ax = plot_isolines(isolines_grouped, starts, minute_range[1], boundaries)
def plot_isolines_px(isolines_grouped, starts, m_range, boundary):
geo_df = isolines_grouped[(isolines_grouped.range == m_range)]
geo_df.crs = CRS(crs)
geo_df = geo_df.to_crs(osm_crs)
bnd = boundaries.copy()
bnd.crs = CRS(crs)
bnd = bnd.to_crs(osm_crs)
# display(geo_df)
geo_df = disjoin_labels(geo_df)
#display(geo_df)
fig = px.choropleth_mapbox(geo_df,
geojson=geo_df.geometry,
locations=geo_df.index,
hover_name="label",
hover_data=["label"],
color=geo_df['label'],
labels={"label": "Scenario"},
center={"lat": bnd.bounds[["miny","maxy"]].mean(axis=1).iloc[0],
"lon": bnd.bounds[["minx","maxx"]].mean(axis=1).iloc[0]},
# mapbox_style="carto-positron",
mapbox_style="open-street-map",
zoom=zoom,
width=950, height=800,
opacity=0.5,
title= f"Reached zone ({dataset}{', '+option if option else ''}) - range: {m_range} minutes"
)
#fig.update_geos(fitbounds="locations")
starts_df = starts.to_crs(osm_crs)
fig.add_trace(go.Scattermapbox(
lat=starts_df.lat,
lon=starts_df.lon,
text=starts_df.name,
name="Starts")
)
if isinstance(bnd.geometry.iloc[0], shapely.geometry.MultiPolygon):
# Get the main polygon
bnd_coords = pd.DataFrame(max(bnd.geometry.iloc[0], key=lambda a: a.area).exterior.coords)
# print(bnd)
else:
bnd_coords = pd.DataFrame(bnd.geometry.iloc[0].exterior.coords)
#[(x, y) for (x, y) in bnd.geometry.iloc[0].exterior.coords])
fig.add_trace(go.Scattermapbox(lat=bnd_coords[1],
lon=bnd_coords[0],
mode='lines',
name="Boundary"))
fig.update_layout(legend=dict(
yanchor="top",
y=0.99,
xanchor="right",
x=0.99
))
return fig
import os
def download_if_nexist(url, filename):
"""
If the (local) file <filename> does not exists, download it from <url>
Parameters
----------
url: str
url to fetch
filename: str
local file to save
Returns
-------
None
"""
if not os.path.isfile(filename):
with urllib.request.urlopen(url) as response:
with open(filename, "wb") as f:
f.write(response.read())
dataset_label=f"{dataset}{'_'+option if option else ''}"
output_path=f"output/{dataset_label}"
try : os.mkdir(output_path)
except FileExistsError: pass
boundaries_filename = "data/boundaries.geojson"
if not os.path.isfile(boundaries_filename):
zipcode_boundaries_filename = "data/zipcode_boundaries_shapefile_3812.zip"
# https://www.geo.be/catalog/details/9738c7c0-5255-11ea-8895-34e12d0f0423?l=fr
download_if_nexist("https://bgu.bpost.be/assets/9738c7c0-5255-11ea-8895-34e12d0f0423_x-shapefile_3812.zip",
zipcode_boundaries_filename)
zipcodes_boundaries = gpd.read_file(f"zip://{zipcode_boundaries_filename}/3812")
bru = Polygon(zipcodes_boundaries[zipcodes_boundaries.nouveau_PO.between("1000", "1299")].unary_union.exterior)
bel = Polygon(max(zipcodes_boundaries.unary_union, key=lambda a: a.area).exterior)
all_boundaries= gpd.GeoDataFrame([{"name":"BRU", "geometry": bru},
{"name":"BEL", "geometry": bel}
])
all_boundaries = all_boundaries.set_geometry("geometry")
all_boundaries.crs = zipcodes_boundaries.crs
all_boundaries = all_boundaries.to_crs(CRS(crs))
all_boundaries.to_file(boundaries_filename, driver="GeoJSON")
else:
all_boundaries = gpd.read_file(boundaries_filename)
all_boundaries.plot("name", alpha=0.5)
<AxesSubplot:>
if dataset == "smur":
starts = gpd.GeoDataFrame(columns =
["name", "address", "wk_xl", "wk_dlt", "wk_stjean", "wk_none"], data=[
["Paul Brien", "Av. Britsiers 5, 1030 Schaerbeek", True, True, True, True ],
["Saint-Pierre", "Rue Haute 290, 1000 Bruxelles", True, True, True, True ],
["Sainte-Elisabeth", "Avenue De Fré 206, 1180 Uccle", True, True, True, True ],
["Saint-Luc", "Avenue Hippocrate, 1200 Woluwe-Saint-Lambert", True, True, True, True ],
["Ărasme (ULB)", "Rue Meylemeersch 56, 1070 Anderlecht", True, True, True, True ],
["HM", "Rue Bruyn 1, 1120 Bruxelles", True, True, True, True ],
["AZVUB", "Avenue du Laerbeek 101, 1090 Jette", True, True, True, True ],
["Ixelles", "Rue Jean Paquot 63, 1050 Ixelles", True, False, False, False],
["Delta", "Boulevard du Triomphe 201, 1160 Auderghem", False, True, False, False],
["Saint-Jean", "Rue du Marais 104, 1000 Bruxelles", False, False, True, False],
])
subsets = ["wk_xl", "wk_dlt", "wk_stjean"] #, "wk_none"
minute_range = list(range(16)) # list(range(21)) #[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 16, 18]
if dataset == "smur":
boundaries = all_boundaries[all_boundaries.name == "BRU"]
zoom = 10.5
if dataset == "maternity":
starts = pd.read_csv("data/maternities.csv", usecols=["name", "address"])
starts = gpd.GeoDataFrame(starts)
starts["current"] = True
# https://kce.fgov.be/fr/une-r%C3%A9duction-du-nombre-de-maternit%C3%A9s-est-possible-dans-notre-pays
# From https://www.google.com/maps/d/u/0/viewer?mid=1BuSrXrmRLzqMNVNLR-hSq0PL53o9GaB9&ll=50.83934807412515%2C4.471731899999895&z=9
# https://plus.lesoir.be/273299/article/2020-01-16/voici-les-17-maternites-belges-qui-pourraient-fermer-leurs-portes-carte
# AZ Sint-Jan (Ostende)
# AZ Zeno (Knokke-Heist)
# AZ Delta (Turnhout)
# Sint-Andries (Tielt)
# St. Vincentius (Deinze)
# A.Z. (Audenarde)
# AZ Delta (Menin)
# AZ Glorieux (Renaix)
# Site St-Joseph, du CHR de Mons
# CH Haute Senne, site le Tilleriau (Soignies)
# CHU André Vésale - site de Lobbes (Montigny-le-Tilleul)
# C.H. Jolimont (Lobbes)
# CHR Sambre et Meuse (Auvelais)
# Centre Hospitalier RĂ©gional de Huy
# CHC - Clinique Sainte-Elisabeth (Verviers)
# St. Nikolaus Hospital (Eupen)
# ASZ (Grammont) - déjà fermé depuis l'été
if dataset == "maternity":
should_go_names = ["ST. JAN BRUGGE", "Zeno", "Torhout", "Tielt", "Deinze", "Oudenaarde", "Menen", "Glorieux",
"CHR Mons", "Tilleriau", "Vesale", "Lobbes", "Sambre", "Huy", "VERVIERS, CENTRE HOSPITALIER CHRETIEN",
"Eupen"]
starts["future"] = starts.name.str.upper().apply(lambda x: sum([sgn.lower() in x.lower() for sgn in should_go_names]) == 0)
assert starts[~starts.future].shape[0] == len(should_go_names)
subsets = ["current", "future"]
minute_range = [0, 5, 10, 15, 20, 25, 30, 35, 40, 45]
if dataset == "maternity":
boundaries = all_boundaries[all_boundaries.name == "BEL"]
zoom = 7
if dataset == "vaccin":
starts = gpd.GeoDataFrame(columns =
["name", "address", "current", "with_WB"], data=[
["Pachéco", "Boulevard Pachéco, 42, 1000 Bruxelles", True, True],
["Heysel ", "Avenue Impératrice Charlotte, 6, 1020 Laeken", True, True],
["WSP", "DrÚve des Shetlands, 15, 1150 Woluwé-Saint-Pierre", True, True],
["Schaerbeek", "Avenue du Suffrage Universel, 1030 Schaerbeek", True, True],
["Molenbeek", "Chaussée de Gand, 696, 1080 Molenbeek", True, True],
["Forest", "Avenue Jupiter, 201, 1190 Forest", True, True],
["Anderlecht", "Avenue Théo Verbeeck, 2, 1070 Anderlecht", True, True],
["WSL", "Avenue des Vaillants, 2, 1200 Woluwé-Saint-Lambert", True, True],
["Uccle", "Rue Ăgide Van Ophem, 110, 1180 Uccle", True, True],
["HM", "Rue Bruyn, 200 Ă 1120 Neder-over-Heembeek", True, True],
["Watermael", "Boulevard du Souverain 25, 1170 Watermael-Boitsfort", False, True]
])
subsets = ["current", "with_WB"]
minute_range = list(range(21)) #[0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20]
if dataset == "vaccin":
boundaries = all_boundaries[all_boundaries.name == "BRU"]
zoom = 10.5
starts["label"] = starts.apply(lambda row: ";".join(sb for sb in subsets if row[sb]), axis=1)
starts["label"] = starts["label"].str.replace(";".join(subsets), "all scenarii")
starts
name | address | wk_xl | wk_dlt | wk_stjean | wk_none | label | |
---|---|---|---|---|---|---|---|
0 | Paul Brien | Av. Britsiers 5, 1030 Schaerbeek | True | True | True | True | all scenarii |
1 | Saint-Pierre | Rue Haute 290, 1000 Bruxelles | True | True | True | True | all scenarii |
2 | Sainte-Elisabeth | Avenue De Fré 206, 1180 Uccle | True | True | True | True | all scenarii |
3 | Saint-Luc | Avenue Hippocrate, 1200 Woluwe-Saint-Lambert | True | True | True | True | all scenarii |
4 | Ărasme (ULB) | Rue Meylemeersch 56, 1070 Anderlecht | True | True | True | True | all scenarii |
5 | HM | Rue Bruyn 1, 1120 Bruxelles | True | True | True | True | all scenarii |
6 | AZVUB | Avenue du Laerbeek 101, 1090 Jette | True | True | True | True | all scenarii |
7 | Ixelles | Rue Jean Paquot 63, 1050 Ixelles | True | False | False | False | wk_xl |
8 | Delta | Boulevard du Triomphe 201, 1160 Auderghem | False | True | False | False | wk_dlt |
9 | Saint-Jean | Rue du Marais 104, 1000 Bruxelles | False | False | True | False | wk_stjean |
starts.label.value_counts()
all scenarii 7 wk_xl 1 wk_dlt 1 wk_stjean 1 Name: label, dtype: int64
if not from_files:
starts = geocode_points(starts)
starts.to_pickle(f"data/starts_{dataset_label}.pkl.gz")
else:
starts = pd.read_pickle(f"data/starts_{dataset_label}.pkl.gz")
# fixing a bug that makes geodataframe being loaded from pickle to be not plottable
# after loading from pickle: starts.name.dtype is np.dtype("O") == False
n = starts.name.str.replace("willnotappear", "").copy()
starts = starts.drop("name", axis=1).assign(name=n)
# after those lines : # starts.name.dtype is np.dtype("O") == True
starts
address | wk_xl | wk_dlt | wk_stjean | wk_none | label | lat | lon | geometry | name | |
---|---|---|---|---|---|---|---|---|---|---|
0 | Av. Britsiers 5, 1030 Schaerbeek | True | True | True | True | all scenarii | 50.869380 | 4.385800 | POINT (488225.023 6598221.000) | Paul Brien |
1 | Rue Haute 290, 1000 Bruxelles | True | True | True | True | all scenarii | 50.835730 | 4.348100 | POINT (484028.278 6592287.540) | Saint-Pierre |
2 | Avenue De Fré 206, 1180 Uccle | True | True | True | True | all scenarii | 50.805300 | 4.370630 | POINT (486536.306 6586925.541) | Sainte-Elisabeth |
3 | Avenue Hippocrate, 1200 Woluwe-Saint-Lambert | True | True | True | True | all scenarii | 50.854470 | 4.447360 | POINT (495077.851 6595591.411) | Saint-Luc |
4 | Rue Meylemeersch 56, 1070 Anderlecht | True | True | True | True | all scenarii | 50.812280 | 4.264400 | POINT (474710.837 6588155.162) | Ărasme (ULB) |
5 | Rue Bruyn 1, 1120 Bruxelles | True | True | True | True | all scenarii | 50.907530 | 4.390220 | POINT (488717.055 6604953.121) | HM |
6 | Avenue du Laerbeek 101, 1090 Jette | True | True | True | True | all scenarii | 50.886120 | 4.309440 | POINT (479724.666 6601174.337) | AZVUB |
7 | Rue Jean Paquot 63, 1050 Ixelles | True | False | False | False | wk_xl | 50.824562 | 4.379648 | POINT (487540.174 6590319.230) | Ixelles |
8 | Boulevard du Triomphe 201, 1160 Auderghem | False | True | False | False | wk_dlt | 50.816680 | 4.399960 | POINT (489801.307 6588930.375) | Delta |
9 | Rue du Marais 104, 1000 Bruxelles | False | False | True | False | wk_stjean | 50.853920 | 4.360510 | POINT (485409.753 6595494.427) | Saint-Jean |
fig = px.scatter_mapbox(starts,
lat="lat",
lon="lon",
hover_name="name",
color='label',
size_max=120,
zoom=10,
mapbox_style="open-street-map",
#opacity=0.2
)
bd = boundaries.to_crs(osm_crs)
fig_bd = px.choropleth_mapbox(bd,
geojson=bd.geometry,
locations=bd.index,
#labels="boundary",
#hover_name="name",
#fill_color="none",
#zoom=10,
opacity=0.2
)
fig.add_trace(fig_bd.data[0])
fig.show()
ax = starts.plot("label", legend=True, figsize=(10,10), markersize=200)
boundaries.boundary.plot(ax=ax, color="C2", linewidth=5)
add_basemap(ax, alpha=0.5)
plt.axis("off")
(471086.2557179927, 500291.65850749245, 6578273.415749143, 6607401.994184209)
from_files=True
if not from_files:
isolines = get_isolines(starts, minute_range, add_options)
isolines.to_pickle(f"data/isolines_{dataset_label}.pkl.gz")
else :
isolines = pd.read_pickle(f"data/isolines_{dataset_label}.pkl.gz")
# fixing a bug that makes geodataframe being loaded from pickle to be not plottable
n = isolines.name.str.replace("willnotappear", "").copy()
isolines = isolines.drop("name", axis=1).assign(name=n)
isolines
address | wk_xl | wk_dlt | wk_stjean | wk_none | label | lat | lon | range | geometry | name | |
---|---|---|---|---|---|---|---|---|---|---|---|
0 | Av. Britsiers 5, 1030 Schaerbeek | True | True | True | True | all scenarii | 50.86938 | 4.38580 | 0 | POLYGON ((488317.974 6598717.030, 488508.999 6... | Paul Brien |
1 | Av. Britsiers 5, 1030 Schaerbeek | True | True | True | True | all scenarii | 50.86938 | 4.38580 | 1 | POLYGON ((487706.497 6599201.623, 487897.521 6... | Paul Brien |
2 | Av. Britsiers 5, 1030 Schaerbeek | True | True | True | True | all scenarii | 50.86938 | 4.38580 | 2 | POLYGON ((487094.907 6599686.245, 487362.519 6... | Paul Brien |
3 | Av. Britsiers 5, 1030 Schaerbeek | True | True | True | True | all scenarii | 50.86938 | 4.38580 | 3 | POLYGON ((486712.747 6598838.219, 486751.041 6... | Paul Brien |
4 | Av. Britsiers 5, 1030 Schaerbeek | True | True | True | True | all scenarii | 50.86938 | 4.38580 | 4 | POLYGON ((486177.746 6598717.030, 486368.770 6... | Paul Brien |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
155 | Rue du Marais 104, 1000 Bruxelles | False | False | True | False | wk_stjean | 50.85392 | 4.36051 | 11 | POLYGON ((476011.605 6599807.272, 476126.264 6... | Saint-Jean |
156 | Rue du Marais 104, 1000 Bruxelles | False | False | True | False | wk_stjean | 50.85392 | 4.36051 | 12 | POLYGON ((473030.580 6600655.397, 473909.559 6... | Saint-Jean |
157 | Rue du Marais 104, 1000 Bruxelles | False | False | True | False | wk_stjean | 50.85392 | 4.36051 | 13 | POLYGON ((471196.035 6601139.928, 471387.171 6... | Saint-Jean |
158 | Rue du Marais 104, 1000 Bruxelles | False | False | True | False | wk_stjean | 50.85392 | 4.36051 | 14 | POLYGON ((467221.373 6600655.397, 468100.351 6... | Saint-Jean |
159 | Rue du Marais 104, 1000 Bruxelles | False | False | True | False | wk_stjean | 50.85392 | 4.36051 | 15 | POLYGON ((464163.871 6600655.397, 468100.351 6... | Saint-Jean |
160 rows Ă 11 columns
isolines_grouped = get_isolines_grouped(isolines, subsets)
isolines_grouped
/home/vb/anaconda3/envs/gis/lib/python3.8/site-packages/pandas/core/dtypes/cast.py:1983: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. /home/vb/anaconda3/envs/gis/lib/python3.8/site-packages/pandas/core/dtypes/cast.py:1983: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. /home/vb/anaconda3/envs/gis/lib/python3.8/site-packages/pandas/core/dtypes/cast.py:1983: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry.
range | geometry | subset | label | |
---|---|---|---|---|
0 | 0 | MULTIPOLYGON (((486617.235 6587097.114, 486712... | wk_xl | wk_xl (0) |
1 | 1 | MULTIPOLYGON (((485699.963 6586855.258, 485986... | wk_xl | wk_xl (1) |
2 | 2 | MULTIPOLYGON (((484782.690 6586855.258, 484916... | wk_xl | wk_xl (2) |
3 | 3 | MULTIPOLYGON (((473947.741 6589516.591, 474138... | wk_xl | wk_xl (3) |
4 | 4 | MULTIPOLYGON (((483215.757 6586673.827, 483273... | wk_xl | wk_xl (4) |
5 | 5 | MULTIPOLYGON (((471196.035 6589516.591, 471387... | wk_xl | wk_xl (5) |
6 | 6 | MULTIPOLYGON (((482298.484 6585948.142, 482355... | wk_xl | wk_xl (6) |
7 | 7 | MULTIPOLYGON (((481572.347 6585403.920, 481667... | wk_xl | wk_xl (7) |
8 | 8 | MULTIPOLYGON (((480808.028 6585403.920, 480903... | wk_xl | wk_xl (8) |
9 | 9 | POLYGON ((492082.418 6591361.934, 492063.319 6... | wk_xl | wk_xl (9) |
10 | 10 | POLYGON ((493209.910 6588367.281, 493286.386 6... | wk_xl | wk_xl (10) |
11 | 11 | POLYGON ((497681.502 6585887.730, 497719.685 6... | wk_xl | wk_xl (11) |
12 | 12 | POLYGON ((494776.876 6587218.115, 494700.478 6... | wk_xl | wk_xl (12) |
13 | 13 | POLYGON ((501095.696 6581656.110, 501121.163 6... | wk_xl | wk_xl (13) |
14 | 14 | POLYGON ((500356.738 6584194.782, 500471.391 6... | wk_xl | wk_xl (14) |
15 | 15 | POLYGON ((499745.254 6585162.114, 500203.890 6... | wk_xl | wk_xl (15) |
0 | 0 | MULTIPOLYGON (((486617.235 6587097.114, 486712... | wk_dlt | wk_dlt (0) |
1 | 1 | MULTIPOLYGON (((485699.963 6586855.258, 485986... | wk_dlt | wk_dlt (1) |
2 | 2 | MULTIPOLYGON (((474253.536 6589032.744, 474444... | wk_dlt | wk_dlt (2) |
3 | 3 | MULTIPOLYGON (((473947.741 6589516.591, 474138... | wk_dlt | wk_dlt (3) |
4 | 4 | MULTIPOLYGON (((472724.786 6589516.591, 473603... | wk_dlt | wk_dlt (4) |
5 | 5 | MULTIPOLYGON (((471196.035 6589516.591, 471387... | wk_dlt | wk_dlt (5) |
6 | 6 | MULTIPOLYGON (((469973.079 6589516.591, 470852... | wk_dlt | wk_dlt (6) |
7 | 7 | MULTIPOLYGON (((468750.123 6589032.744, 469017... | wk_dlt | wk_dlt (7) |
8 | 8 | MULTIPOLYGON (((480808.028 6585403.920, 480903... | wk_dlt | wk_dlt (8) |
9 | 9 | POLYGON ((495006.218 6591452.791, 495082.637 6... | wk_dlt | wk_dlt (9) |
10 | 10 | POLYGON ((498063.661 6588064.964, 498063.662 6... | wk_dlt | wk_dlt (10) |
11 | 11 | POLYGON ((499821.733 6588669.801, 499821.731 6... | wk_dlt | wk_dlt (11) |
12 | 12 | POLYGON ((505401.620 6597748.106, 505095.937 6... | wk_dlt | wk_dlt (12) |
13 | 13 | POLYGON ((505707.415 6598717.030, 506930.371 6... | wk_dlt | wk_dlt (13) |
14 | 14 | POLYGON ((507312.531 6598353.648, 507541.849 6... | wk_dlt | wk_dlt (14) |
15 | 15 | POLYGON ((510599.350 6597748.106, 511210.828 6... | wk_dlt | wk_dlt (15) |
0 | 0 | MULTIPOLYGON (((486617.235 6587097.114, 486712... | wk_stjean | wk_stjean (0) |
1 | 1 | MULTIPOLYGON (((485699.963 6586855.258, 485986... | wk_stjean | wk_stjean (1) |
2 | 2 | MULTIPOLYGON (((484782.690 6586855.258, 484916... | wk_stjean | wk_stjean (2) |
3 | 3 | MULTIPOLYGON (((484018.370 6586855.258, 484457... | wk_stjean | wk_stjean (3) |
4 | 4 | MULTIPOLYGON (((483215.757 6586673.827, 483273... | wk_stjean | wk_stjean (4) |
5 | 5 | MULTIPOLYGON (((471196.035 6589516.591, 471387... | wk_stjean | wk_stjean (5) |
6 | 6 | MULTIPOLYGON (((482298.484 6585948.142, 482355... | wk_stjean | wk_stjean (6) |
7 | 7 | MULTIPOLYGON (((481572.347 6585403.920, 481667... | wk_stjean | wk_stjean (7) |
8 | 8 | MULTIPOLYGON (((480808.028 6585403.920, 480903... | wk_stjean | wk_stjean (8) |
9 | 9 | POLYGON ((490343.488 6590484.546, 490305.361 6... | wk_stjean | wk_stjean (9) |
10 | 10 | POLYGON ((492923.293 6588820.915, 493209.910 6... | wk_stjean | wk_stjean (10) |
11 | 11 | POLYGON ((493477.410 6588790.922, 493515.704 6... | wk_stjean | wk_stjean (11) |
12 | 12 | POLYGON ((494471.048 6587540.720, 494432.865 6... | wk_stjean | wk_stjean (12) |
13 | 13 | POLYGON ((495044.453 6580628.651, 495044.455 6... | wk_stjean | wk_stjean (13) |
14 | 14 | POLYGON ((496114.568 6583287.967, 496114.569 6... | wk_stjean | wk_stjean (14) |
15 | 15 | POLYGON ((496420.362 6583771.647, 496496.729 6... | wk_stjean | wk_stjean (15) |
ax= isolines[isolines.range == minute_range[5]].plot("name", alpha=0.8, figsize=(10, 10), legend=True)
boundaries.boundary.plot(ax=ax, color="C2", linewidth=5)
add_basemap(ax, alpha=0.5)
plt.axis("off")
(469486.66218560725, 500367.8296280822, 6578162.735096075, 6609726.2878986355)
start_name= starts.name.iloc[-2]
sel = isolines[(isolines.name==start_name) & (isolines.range==minute_range[5])]
ax = sel.plot(alpha=0.2, figsize=(10, 10))
ax = sel.boundary.plot(ax=ax)
starts[(starts.name==start_name)].plot(ax=ax)
add_basemap(ax, alpha=0.5)#, source=ctx.providers.OpenStreetMap.Mapnik)
plt.axis("off")
(486965.04749862105, 495709.42727136327, 6586401.4932128545, 6593722.021967109)
sel = isolines_grouped[(isolines_grouped.range == minute_range[5]) &
(isolines_grouped.subset == subsets[1])]
ax= sel.plot("label", alpha=0.2, figsize=(10, 10), legend=False)
ax = sel.boundary.plot(ax=ax)
starts[starts[subsets[1]]].plot(ax=ax)
boundaries.boundary.plot(ax=ax, color="C1")
add_basemap(ax)
# ax = plot_isolines(isolines_grouped, starts, minute_range[6], boundaries)
plt.rcParams.update({'figure.max_open_warning': len(minute_range)})
for r in tqdm([None] + minute_range):
f = plot_isolines(isolines_grouped, starts, r, boundaries)
plt.axis("off")
r_lab = "none" if r is None else f"{r:02}"
plt.savefig(f"{output_path}/reached_zone_{dataset_label}_{r_lab}.png")
plt.close() # Remove this line to see all plots --> might be long!
0%| | 0/17 [00:00<?, ?it/s]
output_path
'output/smur'
plot_isolines(isolines_grouped, starts, minute_range[11], boundaries)
plt.axis("off")
(469758.7374093791, 501619.1768161061, 6576949.38945664, 6608726.020476712)
plot_isolines(isolines_grouped, starts, None, boundaries)
plt.axis("off")
m_range: None
(469758.7374093791, 501619.1768161061, 6576949.38945664, 6608726.020476712)
# Make a video from PNG images - Assumes ffmpeg is installed on the machine
print(f"Making {output_path}/reached_zone_{dataset_label}.mp4")
os.system(f"rm {output_path}/reached_zone_{dataset_label}.mp4")
os.system(f"cat {output_path}/reached_zone*.png | ffmpeg -f image2pipe -framerate 2 -i - -c:v libx264 -vf format=yuv420p -r 25 {output_path}/reached_zone_{dataset_label}.mp4")
Making output/smur/reached_zone_smur.mp4
ffmpeg version 4.4.2-0ubuntu0.22.04.1 Copyright (c) 2000-2021 the FFmpeg developers built with gcc 11 (Ubuntu 11.2.0-19ubuntu1) configuration: --prefix=/usr --extra-version=0ubuntu0.22.04.1 --toolchain=hardened --libdir=/usr/lib/x86_64-linux-gnu --incdir=/usr/include/x86_64-linux-gnu --arch=amd64 --enable-gpl --disable-stripping --enable-gnutls --enable-ladspa --enable-libaom --enable-libass --enable-libbluray --enable-libbs2b --enable-libcaca --enable-libcdio --enable-libcodec2 --enable-libdav1d --enable-libflite --enable-libfontconfig --enable-libfreetype --enable-libfribidi --enable-libgme --enable-libgsm --enable-libjack --enable-libmp3lame --enable-libmysofa --enable-libopenjpeg --enable-libopenmpt --enable-libopus --enable-libpulse --enable-librabbitmq --enable-librubberband --enable-libshine --enable-libsnappy --enable-libsoxr --enable-libspeex --enable-libsrt --enable-libssh --enable-libtheora --enable-libtwolame --enable-libvidstab --enable-libvorbis --enable-libvpx --enable-libwebp --enable-libx265 --enable-libxml2 --enable-libxvid --enable-libzimg --enable-libzmq --enable-libzvbi --enable-lv2 --enable-omx --enable-openal --enable-opencl --enable-opengl --enable-sdl2 --enable-pocketsphinx --enable-librsvg --enable-libmfx --enable-libdc1394 --enable-libdrm --enable-libiec61883 --enable-chromaprint --enable-frei0r --enable-libx264 --enable-shared libavutil 56. 70.100 / 56. 70.100 libavcodec 58.134.100 / 58.134.100 libavformat 58. 76.100 / 58. 76.100 libavdevice 58. 13.100 / 58. 13.100 libavfilter 7.110.100 / 7.110.100 libswscale 5. 9.100 / 5. 9.100 libswresample 3. 9.100 / 3. 9.100 libpostproc 55. 9.100 / 55. 9.100 Input #0, image2pipe, from 'pipe:': Duration: N/A, bitrate: N/A Stream #0:0: Video: png, rgba(pc), 720x720 [SAR 2835:2835 DAR 1:1], 2 fps, 2 tbr, 2 tbn, 2 tbc Stream mapping: Stream #0:0 -> #0:0 (png (native) -> h264 (libx264)) [libx264 @ 0x55642e36e780] using SAR=1/1 [libx264 @ 0x55642e36e780] using cpu capabilities: MMX2 SSE2 SSE3 Cache64 [libx264 @ 0x55642e36e780] profile High, level 3.1, 4:2:0, 8-bit [libx264 @ 0x55642e36e780] 264 - core 163 r3060 5db6aa6 - H.264/MPEG-4 AVC codec - Copyleft 2003-2021 - http://www.videolan.org/x264.html - options: cabac=1 ref=3 deblock=1:0:0 analyse=0x3:0x113 me=hex subme=7 psy=1 psy_rd=1.00:0.00 mixed_ref=1 me_range=16 chroma_me=1 trellis=1 8x8dct=1 cqm=0 deadzone=21,11 fast_pskip=1 chroma_qp_offset=-2 threads=12 lookahead_threads=2 sliced_threads=0 nr=0 decimate=1 interlaced=0 bluray_compat=0 constrained_intra=0 bframes=3 b_pyramid=2 b_adapt=1 b_bias=0 direct=1 weightb=1 open_gop=0 weightp=2 keyint=250 keyint_min=25 scenecut=40 intra_refresh=0 rc_lookahead=40 rc=crf mbtree=1 crf=23.0 qcomp=0.60 qpmin=0 qpmax=69 qpstep=4 ip_ratio=1.40 aq=1:1.00 Output #0, mp4, to 'output/smur/reached_zone_smur.mp4': Metadata: encoder : Lavf58.76.100 Stream #0:0: Video: h264 (avc1 / 0x31637661), yuv420p(tv, progressive), 720x720 [SAR 1:1 DAR 1:1], q=2-31, 25 fps, 12800 tbn Metadata: encoder : Lavc58.134.100 libx264 Side data: cpb: bitrate max/min/avg: 0/0/0 buffer size: 0 vbv_delay: N/A frame= 200 fps=0.0 q=-1.0 Lsize= 362kB time=00:00:07.88 bitrate= 376.2kbits/s dup=184 drop=0 speed=13.8x video:359kB audio:0kB subtitle:0kB other streams:0kB global headers:0kB muxing overhead: 0.889702% [libx264 @ 0x55642e36e780] frame I:1 Avg QP:18.19 size:104407 [libx264 @ 0x55642e36e780] frame P:50 Avg QP:18.69 size: 5087 [libx264 @ 0x55642e36e780] frame B:149 Avg QP:12.97 size: 53 [libx264 @ 0x55642e36e780] consecutive B-frames: 0.5% 0.0% 1.5% 98.0% [libx264 @ 0x55642e36e780] mb I I16..4: 15.4% 34.0% 50.6% [libx264 @ 0x55642e36e780] mb P I16..4: 0.0% 0.3% 1.3% P16..4: 3.4% 0.4% 0.9% 0.0% 0.0% skip:93.7% [libx264 @ 0x55642e36e780] mb B I16..4: 0.0% 0.0% 0.0% B16..8: 1.7% 0.1% 0.0% direct: 0.0% skip:98.2% L0:47.7% L1:52.2% BI: 0.1% [libx264 @ 0x55642e36e780] 8x8 transform intra:26.5% inter:13.5% [libx264 @ 0x55642e36e780] coded y,uvDC,uvAC intra: 76.0% 65.2% 57.9% inter: 0.5% 0.8% 0.2% [libx264 @ 0x55642e36e780] i16 v,h,dc,p: 76% 7% 17% 1% [libx264 @ 0x55642e36e780] i8 v,h,dc,ddl,ddr,vr,hd,vl,hu: 30% 14% 26% 5% 5% 5% 4% 5% 7% [libx264 @ 0x55642e36e780] i4 v,h,dc,ddl,ddr,vr,hd,vl,hu: 15% 14% 14% 8% 9% 10% 9% 9% 11% [libx264 @ 0x55642e36e780] i8c dc,h,v,p: 61% 15% 16% 7% [libx264 @ 0x55642e36e780] Weighted P-Frames: Y:0.0% UV:0.0% [libx264 @ 0x55642e36e780] ref P L0: 74.2% 16.6% 7.4% 1.8% [libx264 @ 0x55642e36e780] ref B L0: 90.1% 9.5% 0.4% [libx264 @ 0x55642e36e780] ref B L1: 96.4% 3.6% [libx264 @ 0x55642e36e780] kb/s:366.62
0
# Build App
app = JupyterDash(__name__)
app.layout = html.Div([
html.H1(f"Multipoint reachability ({dataset})"),
dcc.Graph(id='graph'),
html.Label([
"Range",
dcc.Slider(
id='range-slider',
min=minute_range[0],
max=minute_range[-1],
step=None,
marks={c: str(c) for c in minute_range},
value=minute_range[-3]
)
])
])
# Define callback to update graph
@app.callback(
Output('graph', 'figure'),
[Input("range-slider", "value")]#, Input('traffic-radio', 'value')]
)
def update_figure(m_range):#, traffic):
return plot_isolines_px(isolines_grouped, starts, m_range, boundaries)#, traffic)
ip = 'http://172.27.0.64/'
print(f"Watch outside Jupyter: http://{ip}:8050/")
app.run_server(mode='inline', debug=True, host=ip)
plot_isolines_px(isolines_grouped, starts, minute_range[5], boundaries)
isolines_grouped_bnd = gpd.overlay(isolines_grouped, boundaries[["geometry"]], how="intersection")#.plot()
isolines_grouped_bnd["reach_ratio"] = isolines_grouped_bnd.area/boundaries.area.iloc[0]
from matplotlib.colors import ListedColormap
pv=isolines_grouped_bnd.pivot(index="range", values="reach_ratio", columns= "subset") #["label", "traffic"])#
fig=pv.plot(figsize=(10, 6), title=f"Reach ratio ({dataset}{', '+option if option else ''})", colormap= ListedColormap(["C1", "C2", "C3"]))
plt.savefig(f"{output_path}/reach_ratio_{dataset_label}.png")
# With plotly
fig = px.line(pv, title=f"Reach ratio ({dataset}{', '+option if option else ''})",
color_discrete_sequence=px.colors.qualitative.D3)
fig.update_layout(
yaxis_title='Reached zone ratio',
xaxis_title='Range (minutes)',
title=f"Reach ratio ({dataset}{', '+option if option else ''})",
hovermode="x",
yaxis_tickformat= ',.2%'
)
fig.update_layout(legend=dict(
yanchor="top",
y=0.99,
xanchor="left",
x=0.01
))
fig.write_html(f"{output_path}/reach_ratio_{dataset_label}.html", include_plotlyjs="cdn")
fig
# https://statbel.fgov.be/fr/open-data/secteurs-statistiques-2020
download_if_nexist("https://statbel.fgov.be/sites/default/files/files/opendata/Statistische%20sectoren/sh_statbel_statistical_sectors_31370_20200101.shp.zip",
"data/stat_sectors_2020.zip")
statistical_sectors = gpd.read_file("zip://data/stat_sectors_2020.zip/sh_statbel_statistical_sectors_20200101.shp")
statistical_sectors = statistical_sectors.to_crs(crs)
#statistical_sectors.shape
# https://statbel.fgov.be/fr/nouvelles/la-population-par-secteur-statistique
download_if_nexist("https://statbel.fgov.be/sites/default/files/files/opendata/bevolking/sectoren/OPENDATA_SECTOREN_2020.xlsx",
"data/OPENDATA_SECTOREN_2020.xlsx")
population_per_sector = pd.read_excel("data/OPENDATA_SECTOREN_2020.xlsx")
pop_per_sector_geo = population_per_sector[["TX_DESCR_FR", "TX_DESCR_SECTOR_FR", "CD_SECTOR", "POPULATION"]]
pop_per_sector_geo = pop_per_sector_geo.merge(statistical_sectors[["CS01012020", "geometry"]].rename({"CS01012020":"CD_SECTOR"}, axis=1)).set_geometry("geometry")
pop_per_sector_geo.plot("POPULATION")
<AxesSubplot:>
It might be interessting to consider grid data, grouping population per 1km square. In this example, statistical sectors give a beter precision, so we do not use it.
# # https://statbel.fgov.be/fr/open-data/population-selon-la-grille-km2-grid-2020
# download_if_nexist("https://statbel.fgov.be/sites/default/files/files/opendata/Pop_GRID/TF_POPULATION_GRID_3035_20200101_geojson.zip",
# "data/pop_grid.zip")
# pop_per_sector_geo = gpd.read_file("zip://data/pop_grid.zip/TF_POPULATION_GRID_3035_20200101.geojson")
# pop_per_sector_geo = pop_per_sector_geo[["geometry", "ms_population_20200101", "grd_floaid"]]
# pop_per_sector_geo = pop_per_sector_geo.rename(columns={"ms_population_20200101": "POPULATION",
# "grd_floaid": "CD_SECTOR"})
# pop_per_sector_geo = pop_per_sector_geo.to_crs(crs)
# pop_per_sector_geo
# Keep only sector within boundaries
pop_per_sector_geo = gpd.sjoin(pop_per_sector_geo, boundaries.buffer(100).reset_index().set_geometry(0), op="within")
# pop_per_sector_geo.POPULATION.sort_values()
# pop_per_sector_geo.area.sort_values().astype(int)#.hist()
pop_per_sector_geo.plot("POPULATION")
<AxesSubplot:>
# isolines_grouped_bnd.assign(geometry=isolines_grouped_bnd.buffer(10))
security_buffer = 10
isolines_grouped_pop = isolines_grouped_bnd
for op in ['intersects', 'contains']:
isolines_grouped_op = gpd.sjoin(
#isolines_grouped_bnd[["subset", "range", "geometry"]],
isolines_grouped_bnd.assign(geometry=isolines_grouped_bnd.buffer(security_buffer))[["subset", "range", "geometry"]],
pop_per_sector_geo[["geometry", "POPULATION"]],
op=op)
isolines_grouped_op = isolines_grouped_op.groupby(["subset", "range"]).POPULATION.sum().rename("POP_"+op).reset_index()
isolines_grouped_pop = isolines_grouped_pop.merge(isolines_grouped_op, how="left")
isolines_grouped_pop
range | subset | label | geometry | reach_ratio | POP_intersects | POP_contains | |
---|---|---|---|---|---|---|---|
0 | 0 | wk_xl | wk_xl (0) | MULTIPOLYGON Z (((486617.235 6587097.114 0.000... | 0.004646 | 41477 | NaN |
1 | 1 | wk_xl | wk_xl (1) | MULTIPOLYGON Z (((485699.963 6586855.258 0.000... | 0.013392 | 72069 | NaN |
2 | 2 | wk_xl | wk_xl (2) | MULTIPOLYGON Z (((484782.690 6586855.258 0.000... | 0.038319 | 134367 | 8203.0 |
3 | 3 | wk_xl | wk_xl (3) | MULTIPOLYGON Z (((473947.741 6589516.591 0.000... | 0.084320 | 251772 | 50275.0 |
4 | 4 | wk_xl | wk_xl (4) | MULTIPOLYGON Z (((472724.786 6589516.591 0.000... | 0.164261 | 398218 | 113472.0 |
5 | 5 | wk_xl | wk_xl (5) | MULTIPOLYGON Z (((472494.982 6589516.591 0.000... | 0.263241 | 587283 | 208604.0 |
6 | 6 | wk_xl | wk_xl (6) | MULTIPOLYGON Z (((472564.698 6589580.943 0.000... | 0.391179 | 796436 | 400301.0 |
7 | 7 | wk_xl | wk_xl (7) | MULTIPOLYGON Z (((493553.819 6599807.593 0.000... | 0.521778 | 947700 | 570897.0 |
8 | 8 | wk_xl | wk_xl (8) | MULTIPOLYGON Z (((473671.769 6591103.149 0.000... | 0.657004 | 1083504 | 800134.0 |
9 | 9 | wk_xl | wk_xl (9) | MULTIPOLYGON Z (((498063.662 6589546.327 0.000... | 0.770608 | 1184193 | 936152.0 |
10 | 10 | wk_xl | wk_xl (10) | MULTIPOLYGON Z (((476482.114 6596873.951 0.000... | 0.847635 | 1210670 | 1077733.0 |
11 | 11 | wk_xl | wk_xl (11) | MULTIPOLYGON Z (((490695.520 6580882.981 0.000... | 0.901917 | 1213891 | 1175911.0 |
12 | 12 | wk_xl | wk_xl (12) | POLYGON Z ((494776.876 6587218.115 0.000, 4947... | 0.959866 | 1215012 | 1207596.0 |
13 | 13 | wk_xl | wk_xl (13) | POLYGON Z ((498075.863 6584378.605 0.000, 4980... | 0.960386 | 1215012 | 1207596.0 |
14 | 14 | wk_xl | wk_xl (14) | POLYGON Z ((497467.835 6584120.443 0.000, 4974... | 0.991851 | 1215012 | 1211867.0 |
15 | 15 | wk_xl | wk_xl (15) | POLYGON Z ((496917.071 6583885.914 0.000, 4969... | 0.992651 | 1215012 | 1214954.0 |
16 | 0 | wk_dlt | wk_dlt (0) | MULTIPOLYGON Z (((486617.235 6587097.114 0.000... | 0.005697 | 39278 | NaN |
17 | 1 | wk_dlt | wk_dlt (1) | MULTIPOLYGON Z (((485699.963 6586855.258 0.000... | 0.015875 | 72488 | NaN |
18 | 2 | wk_dlt | wk_dlt (2) | MULTIPOLYGON Z (((474253.536 6589032.744 0.000... | 0.043314 | 127170 | 7569.0 |
19 | 3 | wk_dlt | wk_dlt (3) | MULTIPOLYGON Z (((473947.741 6589516.591 0.000... | 0.088788 | 239388 | 45786.0 |
20 | 4 | wk_dlt | wk_dlt (4) | MULTIPOLYGON Z (((472724.786 6589516.591 0.000... | 0.174936 | 414165 | 110005.0 |
21 | 5 | wk_dlt | wk_dlt (5) | MULTIPOLYGON Z (((472494.982 6589516.591 0.000... | 0.289111 | 635387 | 203872.0 |
22 | 6 | wk_dlt | wk_dlt (6) | MULTIPOLYGON Z (((472564.698 6589580.943 0.000... | 0.443642 | 848266 | 418355.0 |
23 | 7 | wk_dlt | wk_dlt (7) | MULTIPOLYGON Z (((493553.819 6599807.593 0.000... | 0.590395 | 981487 | 637436.0 |
24 | 8 | wk_dlt | wk_dlt (8) | MULTIPOLYGON Z (((473671.769 6591103.149 0.000... | 0.715595 | 1101532 | 846803.0 |
25 | 9 | wk_dlt | wk_dlt (9) | MULTIPOLYGON Z (((476530.027 6596779.296 0.000... | 0.815347 | 1187831 | 972452.0 |
26 | 10 | wk_dlt | wk_dlt (10) | MULTIPOLYGON Z (((476482.114 6596873.951 0.000... | 0.879264 | 1210670 | 1098180.0 |
27 | 11 | wk_dlt | wk_dlt (11) | MULTIPOLYGON Z (((490695.520 6580882.981 0.000... | 0.953038 | 1213891 | 1181744.0 |
28 | 12 | wk_dlt | wk_dlt (12) | POLYGON Z ((490272.766 6580680.180 0.000, 4902... | 0.995431 | 1215012 | 1211916.0 |
29 | 13 | wk_dlt | wk_dlt (13) | POLYGON Z ((490153.384 6580624.282 0.000, 4901... | 0.995469 | 1215012 | 1211916.0 |
30 | 14 | wk_dlt | wk_dlt (14) | POLYGON Z ((490095.702 6580597.450 0.000, 4900... | 0.996354 | 1215012 | 1211916.0 |
31 | 15 | wk_dlt | wk_dlt (15) | POLYGON Z ((487973.997 6580326.428 0.000, 4879... | 0.997144 | 1215012 | 1215003.0 |
32 | 0 | wk_stjean | wk_stjean (0) | MULTIPOLYGON Z (((486617.235 6587097.114 0.000... | 0.004919 | 40482 | NaN |
33 | 1 | wk_stjean | wk_stjean (1) | MULTIPOLYGON Z (((485699.963 6586855.258 0.000... | 0.013364 | 62493 | NaN |
34 | 2 | wk_stjean | wk_stjean (2) | MULTIPOLYGON Z (((484782.690 6586855.258 0.000... | 0.037814 | 115059 | 6531.0 |
35 | 3 | wk_stjean | wk_stjean (3) | MULTIPOLYGON Z (((484018.370 6586855.258 0.000... | 0.083721 | 240034 | 42880.0 |
36 | 4 | wk_stjean | wk_stjean (4) | MULTIPOLYGON Z (((483215.757 6586673.827 0.000... | 0.166324 | 407187 | 103639.0 |
37 | 5 | wk_stjean | wk_stjean (5) | MULTIPOLYGON Z (((472494.982 6589516.591 0.000... | 0.265980 | 591916 | 213934.0 |
38 | 6 | wk_stjean | wk_stjean (6) | MULTIPOLYGON Z (((472564.698 6589580.943 0.000... | 0.391895 | 803034 | 404368.0 |
39 | 7 | wk_stjean | wk_stjean (7) | MULTIPOLYGON Z (((493553.819 6599807.593 0.000... | 0.517370 | 961334 | 572539.0 |
40 | 8 | wk_stjean | wk_stjean (8) | MULTIPOLYGON Z (((473671.769 6591103.149 0.000... | 0.647418 | 1072816 | 805849.0 |
41 | 9 | wk_stjean | wk_stjean (9) | MULTIPOLYGON Z (((498063.662 6589546.327 0.000... | 0.759945 | 1175234 | 956737.0 |
42 | 10 | wk_stjean | wk_stjean (10) | MULTIPOLYGON Z (((476482.114 6596873.951 0.000... | 0.834904 | 1206451 | 1058014.0 |
43 | 11 | wk_stjean | wk_stjean (11) | MULTIPOLYGON Z (((497631.154 6585740.184 0.000... | 0.899680 | 1213891 | 1167395.0 |
44 | 12 | wk_stjean | wk_stjean (12) | POLYGON Z ((494471.048 6587540.720 0.000, 4944... | 0.959214 | 1215012 | 1207596.0 |
45 | 13 | wk_stjean | wk_stjean (13) | POLYGON Z ((490272.766 6580680.180 0.000, 4902... | 0.960128 | 1215012 | 1207596.0 |
46 | 14 | wk_stjean | wk_stjean (14) | POLYGON Z ((490272.766 6580680.180 0.000, 4902... | 0.991851 | 1215012 | 1211867.0 |
47 | 15 | wk_stjean | wk_stjean (15) | POLYGON Z ((496420.362 6583771.647 0.000, 4964... | 0.992651 | 1215012 | 1214954.0 |
pv_pop_max=isolines_grouped_pop.pivot(index="range", values="POP_intersects", columns= "subset").fillna(0) #["label", "traffic"])#
pv_pop_min=isolines_grouped_pop.pivot(index="range", values="POP_contains", columns= "subset").fillna(0) #["label", "traffic"])#
ax = pv_pop_max.plot(figsize=(10, 6), title=f"Population reach ratio ({dataset}{', '+option if option else ''})")
pv_pop_min.fillna(0).plot(ax=ax, alpha=0.5, color=["C0", "C1", "C2", "C3"], legend=False)
for c in range(pv_pop_min.shape[1]):
ax.fill_between(pv_pop_max.index, pv_pop_max.iloc[:,c], pv_pop_min.iloc[:,c], color=f"C{c}", alpha=0.1)
plt.savefig(f"{output_path}/pop_reach_ratio_{dataset_label}.png") #set_index("label", "range").reach_ratio.plot()
isolines_grouped_pop
range | subset | label | geometry | reach_ratio | POP_intersects | POP_contains | |
---|---|---|---|---|---|---|---|
0 | 0 | wk_xl | wk_xl (0) | MULTIPOLYGON Z (((486617.235 6587097.114 0.000... | 0.004646 | 41477 | NaN |
1 | 1 | wk_xl | wk_xl (1) | MULTIPOLYGON Z (((485699.963 6586855.258 0.000... | 0.013392 | 72069 | NaN |
2 | 2 | wk_xl | wk_xl (2) | MULTIPOLYGON Z (((484782.690 6586855.258 0.000... | 0.038319 | 134367 | 8203.0 |
3 | 3 | wk_xl | wk_xl (3) | MULTIPOLYGON Z (((473947.741 6589516.591 0.000... | 0.084320 | 251772 | 50275.0 |
4 | 4 | wk_xl | wk_xl (4) | MULTIPOLYGON Z (((472724.786 6589516.591 0.000... | 0.164261 | 398218 | 113472.0 |
5 | 5 | wk_xl | wk_xl (5) | MULTIPOLYGON Z (((472494.982 6589516.591 0.000... | 0.263241 | 587283 | 208604.0 |
6 | 6 | wk_xl | wk_xl (6) | MULTIPOLYGON Z (((472564.698 6589580.943 0.000... | 0.391179 | 796436 | 400301.0 |
7 | 7 | wk_xl | wk_xl (7) | MULTIPOLYGON Z (((493553.819 6599807.593 0.000... | 0.521778 | 947700 | 570897.0 |
8 | 8 | wk_xl | wk_xl (8) | MULTIPOLYGON Z (((473671.769 6591103.149 0.000... | 0.657004 | 1083504 | 800134.0 |
9 | 9 | wk_xl | wk_xl (9) | MULTIPOLYGON Z (((498063.662 6589546.327 0.000... | 0.770608 | 1184193 | 936152.0 |
10 | 10 | wk_xl | wk_xl (10) | MULTIPOLYGON Z (((476482.114 6596873.951 0.000... | 0.847635 | 1210670 | 1077733.0 |
11 | 11 | wk_xl | wk_xl (11) | MULTIPOLYGON Z (((490695.520 6580882.981 0.000... | 0.901917 | 1213891 | 1175911.0 |
12 | 12 | wk_xl | wk_xl (12) | POLYGON Z ((494776.876 6587218.115 0.000, 4947... | 0.959866 | 1215012 | 1207596.0 |
13 | 13 | wk_xl | wk_xl (13) | POLYGON Z ((498075.863 6584378.605 0.000, 4980... | 0.960386 | 1215012 | 1207596.0 |
14 | 14 | wk_xl | wk_xl (14) | POLYGON Z ((497467.835 6584120.443 0.000, 4974... | 0.991851 | 1215012 | 1211867.0 |
15 | 15 | wk_xl | wk_xl (15) | POLYGON Z ((496917.071 6583885.914 0.000, 4969... | 0.992651 | 1215012 | 1214954.0 |
16 | 0 | wk_dlt | wk_dlt (0) | MULTIPOLYGON Z (((486617.235 6587097.114 0.000... | 0.005697 | 39278 | NaN |
17 | 1 | wk_dlt | wk_dlt (1) | MULTIPOLYGON Z (((485699.963 6586855.258 0.000... | 0.015875 | 72488 | NaN |
18 | 2 | wk_dlt | wk_dlt (2) | MULTIPOLYGON Z (((474253.536 6589032.744 0.000... | 0.043314 | 127170 | 7569.0 |
19 | 3 | wk_dlt | wk_dlt (3) | MULTIPOLYGON Z (((473947.741 6589516.591 0.000... | 0.088788 | 239388 | 45786.0 |
20 | 4 | wk_dlt | wk_dlt (4) | MULTIPOLYGON Z (((472724.786 6589516.591 0.000... | 0.174936 | 414165 | 110005.0 |
21 | 5 | wk_dlt | wk_dlt (5) | MULTIPOLYGON Z (((472494.982 6589516.591 0.000... | 0.289111 | 635387 | 203872.0 |
22 | 6 | wk_dlt | wk_dlt (6) | MULTIPOLYGON Z (((472564.698 6589580.943 0.000... | 0.443642 | 848266 | 418355.0 |
23 | 7 | wk_dlt | wk_dlt (7) | MULTIPOLYGON Z (((493553.819 6599807.593 0.000... | 0.590395 | 981487 | 637436.0 |
24 | 8 | wk_dlt | wk_dlt (8) | MULTIPOLYGON Z (((473671.769 6591103.149 0.000... | 0.715595 | 1101532 | 846803.0 |
25 | 9 | wk_dlt | wk_dlt (9) | MULTIPOLYGON Z (((476530.027 6596779.296 0.000... | 0.815347 | 1187831 | 972452.0 |
26 | 10 | wk_dlt | wk_dlt (10) | MULTIPOLYGON Z (((476482.114 6596873.951 0.000... | 0.879264 | 1210670 | 1098180.0 |
27 | 11 | wk_dlt | wk_dlt (11) | MULTIPOLYGON Z (((490695.520 6580882.981 0.000... | 0.953038 | 1213891 | 1181744.0 |
28 | 12 | wk_dlt | wk_dlt (12) | POLYGON Z ((490272.766 6580680.180 0.000, 4902... | 0.995431 | 1215012 | 1211916.0 |
29 | 13 | wk_dlt | wk_dlt (13) | POLYGON Z ((490153.384 6580624.282 0.000, 4901... | 0.995469 | 1215012 | 1211916.0 |
30 | 14 | wk_dlt | wk_dlt (14) | POLYGON Z ((490095.702 6580597.450 0.000, 4900... | 0.996354 | 1215012 | 1211916.0 |
31 | 15 | wk_dlt | wk_dlt (15) | POLYGON Z ((487973.997 6580326.428 0.000, 4879... | 0.997144 | 1215012 | 1215003.0 |
32 | 0 | wk_stjean | wk_stjean (0) | MULTIPOLYGON Z (((486617.235 6587097.114 0.000... | 0.004919 | 40482 | NaN |
33 | 1 | wk_stjean | wk_stjean (1) | MULTIPOLYGON Z (((485699.963 6586855.258 0.000... | 0.013364 | 62493 | NaN |
34 | 2 | wk_stjean | wk_stjean (2) | MULTIPOLYGON Z (((484782.690 6586855.258 0.000... | 0.037814 | 115059 | 6531.0 |
35 | 3 | wk_stjean | wk_stjean (3) | MULTIPOLYGON Z (((484018.370 6586855.258 0.000... | 0.083721 | 240034 | 42880.0 |
36 | 4 | wk_stjean | wk_stjean (4) | MULTIPOLYGON Z (((483215.757 6586673.827 0.000... | 0.166324 | 407187 | 103639.0 |
37 | 5 | wk_stjean | wk_stjean (5) | MULTIPOLYGON Z (((472494.982 6589516.591 0.000... | 0.265980 | 591916 | 213934.0 |
38 | 6 | wk_stjean | wk_stjean (6) | MULTIPOLYGON Z (((472564.698 6589580.943 0.000... | 0.391895 | 803034 | 404368.0 |
39 | 7 | wk_stjean | wk_stjean (7) | MULTIPOLYGON Z (((493553.819 6599807.593 0.000... | 0.517370 | 961334 | 572539.0 |
40 | 8 | wk_stjean | wk_stjean (8) | MULTIPOLYGON Z (((473671.769 6591103.149 0.000... | 0.647418 | 1072816 | 805849.0 |
41 | 9 | wk_stjean | wk_stjean (9) | MULTIPOLYGON Z (((498063.662 6589546.327 0.000... | 0.759945 | 1175234 | 956737.0 |
42 | 10 | wk_stjean | wk_stjean (10) | MULTIPOLYGON Z (((476482.114 6596873.951 0.000... | 0.834904 | 1206451 | 1058014.0 |
43 | 11 | wk_stjean | wk_stjean (11) | MULTIPOLYGON Z (((497631.154 6585740.184 0.000... | 0.899680 | 1213891 | 1167395.0 |
44 | 12 | wk_stjean | wk_stjean (12) | POLYGON Z ((494471.048 6587540.720 0.000, 4944... | 0.959214 | 1215012 | 1207596.0 |
45 | 13 | wk_stjean | wk_stjean (13) | POLYGON Z ((490272.766 6580680.180 0.000, 4902... | 0.960128 | 1215012 | 1207596.0 |
46 | 14 | wk_stjean | wk_stjean (14) | POLYGON Z ((490272.766 6580680.180 0.000, 4902... | 0.991851 | 1215012 | 1211867.0 |
47 | 15 | wk_stjean | wk_stjean (15) | POLYGON Z ((496420.362 6583771.647 0.000, 4964... | 0.992651 | 1215012 | 1214954.0 |
nb_levels=5
geo_df_levels = None
def get_crossed_sectors(geo_df, shape):
"""
Compute which polygons of geo_df are crossed by shape
Parameters
----------
geo_df: geopandas.GeoDataFrame
a GeoDataFrame with a string column "CD_SECTOR"
shape: shapely.geometry.linestring.LineString
a line corresponding typically to a geofg.loc[i].geometry.boundary
Returns
-------
pd.Series of boolean, with geo_df.shape[0] rows
"""
# Simple version. Might be very long, especially for "maternity" use case
# return pop_per_sector_geo.geometry.crosses(shape)
# Hierarchical version. Should be much more efficient
# First tries to see which high level boundaries such as provinces are crossed by "shape".
# Then, only within those "crossed" provinces, see which districts are crossed by "shape", ..., down to statistical sectors
global geo_df_levels
if geo_df_levels is None:
geo_df_levels = multilevel_split(geo_df)
ppsg = geo_df_levels[0]
for i in tqdm(range(nb_levels), leave=False):
crosses = ppsg.geometry.crosses(shape)
ppsg = ppsg[crosses]
if i < nb_levels-1:
ppsg = gpd.GeoDataFrame(ppsg[[f"CD_{i}"]].merge(geo_df_levels[i+1]))
return geo_df.CD_SECTOR.isin(ppsg[f"CD_{nb_levels-1}"])
# Multilevels
# Group sectors in several levels of larger sections
def multilevel_split(geo_df):
geo_df = geo_df.assign(CD_0 = geo_df.CD_SECTOR.str[0:2],
CD_1 = geo_df.CD_SECTOR.str[0:4],
CD_2 = geo_df.CD_SECTOR.str[0:6],
CD_3 = geo_df.CD_SECTOR.str[0:7],
CD_4 = geo_df.CD_SECTOR)
geo_df_levels= {}
for i in tqdm(range(nb_levels), leave=False, desc="Getting multilevel split..."):
df = geo_df.dissolve(f"CD_{i}")[["geometry"]].reset_index()
if i > 0:
df = df.merge(geo_df[[f"CD_{i-1}", f"CD_{i}"]].drop_duplicates())
geo_df_levels[i] = df
return geo_df_levels
if from_files:
isolines_grouped_pop = pd.read_pickle(f"data/isolines_reach_pop_{dataset_label}.pkl.gz")
else:
for i in tqdm(isolines_grouped_pop.index):
isoline_sel=isolines_grouped_pop.loc[i].geometry
# crosses = pop_per_sector_geo.simplify(1000).geometry.crosses(isoline_sel.simplify(1000).boundary)
crosses = get_crossed_sectors(pop_per_sector_geo, isoline_sel.buffer(security_buffer).boundary)
isolines_grouped_pop.loc[i, "POP_crossed"] = pop_per_sector_geo[crosses].POPULATION.sum()
common_part = pop_per_sector_geo[crosses].intersection(isoline_sel)#.plot()
covered_ratio = common_part.area / pop_per_sector_geo[crosses].area
covered_pop = int((pop_per_sector_geo[crosses].POPULATION * covered_ratio).sum())
isolines_grouped_pop.loc[i, "POP_mean"] = isolines_grouped_pop.POP_contains.fillna(0).loc[i] + covered_pop
isolines_grouped_pop.to_pickle(f"data/isolines_reach_pop_{dataset_label}.pkl.gz")
0%| | 0/48 [00:00<?, ?it/s]
Getting multilevel split...: 0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
# Solid line: assumes population is uniformly spread accross all sectors
# Top fill: assumes we reach all the population of a sector as soon as the isolines cross the sector (best case assumption)
# Bottom fill: assumes we reach no one in a sector before the full sector is included in the isoline (worst case assumption)
pv_pop_max= isolines_grouped_pop.pivot(index="range", values="POP_intersects", columns= "subset").fillna(0) / pop_per_sector_geo.POPULATION.sum()
pv_pop_min= isolines_grouped_pop.pivot(index="range", values="POP_contains", columns= "subset").fillna(0) / pop_per_sector_geo.POPULATION.sum()
pv_pop_mean=isolines_grouped_pop.pivot(index="range", values="POP_mean", columns= "subset").fillna(0) / pop_per_sector_geo.POPULATION.sum()
# With Matplotlib
ax = pv_pop_mean.plot(figsize=(10, 6),
title=f"Population reach ratio ({dataset}{', '+option if option else ''})",
colormap= ListedColormap(["C1", "C2", "C3"])) #set_index("label", "range").reach_ratio.plot()
for c in range(pv_pop_mean.shape[1]):
ax.fill_between(pv_pop_max.index, pv_pop_max.iloc[:,c], pv_pop_min.iloc[:,c], color=f"C{c+1}", alpha=0.1)
plt.savefig(f"{output_path}/pop_reach_mean_ratio_{dataset_label}.png") #set_index("label", "range").reach_ratio.plot()
# With Plotly express
fig = go.Figure()
for subs, color in zip(pv_pop_mean.columns[::-1], ["44, 160, 44", "255, 127, 14", "31, 119, 180", "255,0,0"]):
fig.add_trace(
go.Scatter(
name=f'{subs} (mean)',
x=pv_pop_mean.index,
y=pv_pop_mean[subs],
mode='lines',
line=dict(color=f'rgb({color})'),
))
fig.add_trace( go.Scatter(
name=f'{subs} (UB)',
x=pv_pop_max.index,
y=pv_pop_max[subs],
mode='lines',
marker=dict(color=f'rgb({color})'),
line=dict(width=0),
showlegend=False
))#,
fig.add_trace( go.Scatter(
name=f'{subs} (LB)',
x=pv_pop_min.index,
y=pv_pop_min[subs],
marker=dict(color=f'rgb({color})'),
line=dict(width=0),
mode='lines',
fillcolor=f'rgba({color}, 0.2)',
fill='tonexty',
showlegend=False
))
# ])
fig.update_layout(
yaxis_title='Reached population (ratio)',
xaxis_title='Range (minutes)',
title=f"Population reach ratio ({dataset}{', '+option if option else ''})",
hovermode="x",
yaxis_tickformat= ',.2%',
)
fig.update_layout(legend=dict(
yanchor="top",
y=0.99,
xanchor="left",
x=0.01
))
fig.show()
fig.write_html(f"{output_path}/pop_reach_mean_ratio_{dataset_label}.html",include_plotlyjs="cdn" ) #set_index("label", "range").reach_ratio.plot()
# Works for "Vaccin", no traffic
if dataset == "vaccin" and option is None:
aud_wb = statistical_sectors[statistical_sectors.T_MUN_FR.isin(["Auderghem", "Watermael-Boitsfort"])].copy()
il = isolines[(isolines.name == "Watermael") & (isolines.range==3)]
aud_wb["color"] = "red"
aud_wb.loc[aud_wb.crosses(il.boundary.iloc[0]), "color"] ="orange"
aud_wb.loc[aud_wb.within(il.iloc[0].geometry), "color"] ="green"
ax = aud_wb.plot(figsize=(10,10), edgecolor="black", color=aud_wb.color, alpha=0.2)
il.boundary.plot(ax=ax, color="blue")
plt.axis('off')
plt.savefig(f"{output_path}/crosses.png")
def get_osm_polyg(name):
url = f"http://nominatim.openstreetmap.org/search.php?q={name}&format=json&limit=1&polygon_geojson=1"
with urllib.request.urlopen(url) as response:
return shape(json.loads(response.read())[0]["geojson"])
foret = get_osm_polyg("Foret+de+Soignes")
bxl = get_osm_polyg("Bruxelles-Capitale")
inters = bxl.intersection(foret)
df = gpd.GeoDataFrame([{"name":"ForĂȘt de Soignes", "geometry": foret},
{"name":"Bruxelles", "geometry": bxl},
{"name":"Intersection", "geometry": inters}])
df = df.set_geometry("geometry")
df.crs = osm_crs
df.plot("name", legend=True)
print(f"Proportion of 'ForĂȘt de Soignes' in Brussels: {inters.area/bxl.area:.02%}")