Voronoi regions of schools in East Germany. An example using the geovoronoi package (https://pypi.org/project/geovoronoi/).

voronoize.py

"""
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.')