Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Voronoi regions of schools in East Germany. An example using the geovoronoi package (https://pypi.org/project/geovoronoi/).
"""
Voronoi regions of schools in East Germany.
An example using the geovoronoi package (https://pypi.org/project/geovoronoi/).
Feb. 2021
Markus Konrad <markus.konrad@wzb.eu>
"""
import os
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
pd.set_option('display.max_columns', 100)
pd.set_option('display.width', 180)
import fiona
import pyproj
from shapely.geometry import shape
from geovoronoi import coords_to_points, voronoi_regions_from_coords, calculate_polygon_areas
from geovoronoi.plotting import subplot_for_map, plot_polygon_collection_with_color, plot_points,\
plot_voronoi_polys_with_points_in_area
#%%
SCHOOLTYPES = {
'gs': 'primary schools',
'ohne_gym_os': 'schools w/o senior High School classes',
'mit_gym_os': 'schools w/ senior High School classes',
}
COLOR_PUBL = '#0380B5'
COLOR_PRIV = '#9E3173'
INPUT_DATA = 'data/schools.csv'
INPUT_SHAPE = 'data/eastgermany.geojson'
INPUT_SHAPE_UNION = 'data/eastgermany_union.geojson'
FIGSIZE = (8, 12)
#%%
def plot_locations(geom, df):
fig, ax = subplot_for_map(figsize=FIGSIZE)
schools_publ = df.loc[df.traeger == 'oeff', :]
schools_priv = df.loc[df.traeger == 'priv', :]
plot_polygon_collection_with_color(ax, geom, edgecolor='black', facecolor='white')
plot_points(ax, coords_to_points(schools_publ[['x', 'y']].to_numpy()), label='public',
color=COLOR_PUBL, alpha=0.5, markersize=16)
plot_points(ax, coords_to_points(schools_priv[['x', 'y']].to_numpy()), label='private',
color=COLOR_PRIV, alpha=0.5, markersize=16)
return fig, ax
def setup_voronoi_fig():
fig = plt.figure(figsize=FIGSIZE, constrained_layout=False)
figgrid = fig.add_gridspec(nrows=6, ncols=1)
ax_vormap = fig.add_subplot(figgrid[:-1, :])
ax_vormap.set_aspect('equal')
ax_vormap.get_xaxis().set_visible(False)
ax_vormap.get_yaxis().set_visible(False)
for sp in ax_vormap.spines.values():
sp.set_visible(False)
ax_area_ts = fig.add_subplot(figgrid[-1:, :])
return fig, ax_vormap, ax_area_ts
def plot_voronoi(ax, geom, geom_surrounding, region_polys, yearpts, region_pts, region_colors, pts_colors):
plot_voronoi_polys_with_points_in_area(ax, None, region_polys, yearpts, region_pts,
voronoi_color=region_colors,
voronoi_edgecolor='gray',
points_color=pts_colors.tolist(),
plot_voronoi_opts={'alpha': 0.5},
plot_points_opts={'alpha': 0.5}
)
plot_polygon_collection_with_color(ax, geom, edgecolor='black', facecolor=(1, 1, 1, 0), linewidth=0.5)
plot_polygon_collection_with_color(ax, geom_surrounding, edgecolor='black', facecolor=(1, 1, 1, 0))
return ax
def plot_voronoi_area_ts(ax, area_df):
ax.bar(area_df.year, area_df.priv_sum, color=COLOR_PRIV, label='private')
ax.bar(area_df.year, area_df.publ_sum, color=COLOR_PUBL, label='public', bottom=area_df.priv_sum)
return ax
def finish_figure(fig, ax_main, ax_area_ts, year, schooltype, schooltype_lbl, plotfile):
fs = 24 if schooltype == 'gs' else 15
scnd_row = 0.92 if schooltype == 'gs' else 0.935
xpos = [0.025, 0.16, 0.25, 0.41] if schooltype == 'gs' else [0.025, 0.11, 0.162, 0.26]
fig.text(xpos[0], 0.955, 'Public', fontsize=fs, color=COLOR_PUBL)
fig.text(xpos[1], 0.955, 'and', fontsize=fs, color='k')
fig.text(xpos[2], 0.955, 'private', fontsize=fs, color=COLOR_PRIV)
fig.text(xpos[3], 0.955, schooltype_lbl, fontsize=fs, color='k')
fig.text(xpos[0], scnd_row, 'in East Germany (except Berlin)', fontsize=fs, color='k')
fig.text(0.95, 0.05, 'Data: Helbig/Konrad/Nikolai 2018. https://schulenkarte.wzb.eu', ha='right', fontsize=9)
ax_main.set_title(str(year), loc='right', fontdict=dict(fontsize=fs, fontweight='bold'))
if ax_area_ts is not None:
ax_area_ts.set_title('Total covered area in km²')
print(f'>> saving plot to {plotfile}')
#fig.tight_layout()
fig.savefig(plotfile)
plt.close(fig)
#%%
print(f'loading {INPUT_DATA}')
schools = pd.read_csv(INPUT_DATA)
schools = schools[['art_reduziert', 'jahr', 'traeger', 'lng', 'lat']]\
.sort_values(['art_reduziert', 'jahr', 'traeger', 'lng', 'lat'])
n_duplicated = schools.duplicated(['art_reduziert', 'jahr', 'lng', 'lat']).sum()
print(f'dropping {n_duplicated} duplicates')
schools.drop_duplicates(['art_reduziert', 'jahr', 'lng', 'lat'], keep='last', inplace=True)
coords = schools[['lng', 'lat']].to_numpy()
schools.drop(columns=['lng', 'lat'], inplace=True)
#%%
source_crs = pyproj.CRS.from_epsg(4326) # school coord. CRS: standard WGS 84
print(f'loading {INPUT_SHAPE}')
with fiona.open(INPUT_SHAPE) as f:
target_crs = pyproj.CRS.from_dict(f.meta['crs'])
eastgermany = [shape(feat['geometry']) for feat in f]
assert all(shape.is_valid for shape in eastgermany)
print(f'loading {INPUT_SHAPE_UNION}')
with fiona.open(INPUT_SHAPE_UNION) as f:
assert pyproj.CRS.from_dict(f.meta['crs']) == target_crs
assert len(f) == 1
eastgermany_union = shape(f[0]['geometry'])
assert eastgermany_union.is_valid
# transform to shape file CRS (EPSG 3035 - LAEA Europe equal area projection)
print('transforming coordinates to shape file CRS (EPSG 3035 - LAEA Europe equal area projection)')
crs_transformer = pyproj.Transformer.from_crs(source_crs, target_crs)
schools['x'], schools['y'] = crs_transformer.transform(coords[:, 1], coords[:, 0])
del coords
#%%
for schooltype, schooltype_lbl in SCHOOLTYPES.items():
print(f'generating plots for {schooltype_lbl}')
area_df = pd.DataFrame(columns=['year', 'publ_median', 'priv_median', 'publ_sum', 'priv_sum'])
all_years = sorted(schools.jahr.unique())
for year in all_years:
print(f'> year: {year}')
yeardf = schools.loc[(schools.jahr == year) & (schools.art_reduziert == schooltype), :]
plotfile = f'output/voronoi/{schooltype}_{year}.png'
if os.path.exists(plotfile):
print(f'>> file exists: {plotfile} - skipping voronoi plot')
else:
print('>> checking points')
yearpts = coords_to_points(yeardf[['x', 'y']].to_numpy())
yearpts_in_geom = np.array([eastgermany_union.contains(p) for p in yearpts])
yearpts = np.asarray(yearpts, dtype=object)[yearpts_in_geom].tolist()
yeardf = yeardf.loc[yearpts_in_geom, :]
print('>> generating voronoi regions')
region_polys, region_pts = voronoi_regions_from_coords(yearpts, eastgermany_union, per_geom=False)
region_publ = {i_reg: (yeardf.iloc[reg_pts[0]].traeger == 'oeff') if len(reg_pts) == 1 else None
for i_reg, reg_pts in region_pts.items()}
assert None not in region_publ.values()
region_colors = {i_reg: COLOR_PUBL if is_publ else COLOR_PRIV
for i_reg, is_publ in region_publ.items()}
pts_colors = np.repeat(COLOR_PUBL, len(yeardf))
pts_colors[yeardf.traeger == 'priv'] = COLOR_PRIV
region_area = calculate_polygon_areas(region_polys, m2_to_km2=True)
publ_area = np.array([a for i_reg, a in region_area.items() if region_publ[i_reg]])
priv_area = np.array([a for i_reg, a in region_area.items() if not region_publ[i_reg]])
publ_med_area = np.median(publ_area)
priv_med_area = np.median(priv_area)
area_df = pd.concat(
(area_df, pd.DataFrame([[year, np.median(publ_area), np.median(priv_area),
np.sum(publ_area), np.sum(priv_area)]],
columns=area_df.columns)), ignore_index=True
)
print('>> plotting voronoi regions')
fig, ax_vormap, ax_area_ts = setup_voronoi_fig()
plot_voronoi(ax_vormap, eastgermany, eastgermany_union, region_polys, yearpts, region_pts, region_colors,
pts_colors)
print('>> plotting area statistics')
plot_voronoi_area_ts(ax_area_ts, area_df)
ax_area_ts.set_xlim(min(all_years), max(all_years))
finish_figure(fig, ax_vormap, ax_area_ts, year, schooltype, schooltype_lbl, plotfile)
plotfile = f'output/locations/{schooltype}_{year}.png'
if os.path.exists(plotfile):
print(f'>> file exists: {plotfile} - skipping locations plot')
else:
print('>> plotting locations')
fig, ax = plot_locations(eastgermany, yeardf)
finish_figure(fig, ax, None, year, schooltype, schooltype_lbl, plotfile)
#%%
print('done.')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment