Script to show the N values from https://pavelmayer.de/covid/risks/ as a Map. Also attempts to scale the N value with the area of the districts. Also inspired by https://rstats-tips.net/2020/08/09/visualization-of-corona-incidence-in-germany-per-county/

covid_pavel_N_map_plot.py

# -*- coding: utf-8 -*-

import pandas as pd
import geopandas as gpd
import geoplot as gplt

import numpy as np

f = gpd.read_file('./gadm36_DEU_shp/gadm36_DEU_2.shp')


from datetime import datetime

columns = [
        'IdBundesland',
        'Bundesland',
        'IdLandkreis',
        'Landkreis',
        'Meldedatum',
        'AnzahlFall',
        'Bevoelkerung',
        'Refdatum',
        'Datenstand',
    ]

utctimestamparser = lambda x: datetime.utcfromtimestamp(int(x)/1000)

df = pd.read_csv('full-data.csv', usecols = columns, parse_dates=['Meldedatum', 'Refdatum'], date_parser=utctimestamparser)

df["Datenstand"] = pd.to_datetime(df['Datenstand'],format="%d.%m.%Y, %H:%S Uhr")


#Combine all data for Berlin
df = df.assign(
        IdLandkreis = np.where(df['IdBundesland'] == 11, 11000, df['IdLandkreis']),
        Bevoelkerung = np.where(df['IdBundesland'] == 11,
                                df[df['IdBundesland'] == 11].groupby('IdLandkreis').first()['Bevoelkerung'].sum(),
                                df['Bevoelkerung']
                        )
    )

dff = df[['IdBundesland','IdLandkreis','AnzahlFall','Meldedatum']]
dff = dff[dff['AnzahlFall'] > 0]

level = 'IdLandkreis'

dff = dff.set_index([level, 'Meldedatum'])['AnzahlFall']

dff = dff.groupby(level=level).resample('1d',level='Meldedatum').sum()

dff = dff.reindex(pd.MultiIndex.from_product([dff.index.levels[0],dff.index.levels[1]],names=[level,'Meldedatum']), fill_value=0)

sum_7day = dff.groupby(level=level).apply(lambda x: x.rolling(7).sum())
sum_14day = dff.groupby(level=level).apply(lambda x: x.rolling(14).sum())

rwk = sum_7day.groupby(level=level).apply(lambda y: y.rolling(8).apply(
            lambda x: (x[-1] + 5) / (x[0] + 5), raw=True
        )).rename('RwK')

bevoelkerung = df[['IdBundesland','IdLandkreis','Bevoelkerung']].groupby(['IdBundesland','IdLandkreis']).first()['Bevoelkerung']

N = bevoelkerung.groupby(level).sum() / 6.25 / (rwk * sum_14day)
N = N.rename('N')

f['IdLandkreis'] = pd.to_numeric(f['CC_2'])
f = f[~np.isnan(f['IdLandkreis'])]
f['IdLandkreis'] = f['IdLandkreis'].astype(int)
f.set_index('IdLandkreis',inplace=True)

#Fix Göttingen and Osterode am Harz:
f.loc[3159] = f.loc[3152]
f.loc[3159,'geometry'] = f.loc[3156]['geometry'].union(f.loc[3152]['geometry'])
f.drop(3156,inplace=True)
f.drop(3152,inplace=True)


f = f.assign(area = f.area)
normalized_area = f.area / f.area.mean()

N_area = (N * normalized_area).clip(lower=1,upper=99999)
N_area = N_area.rename('N')
N = N.clip(lower=1,upper=99999)

from matplotlib import pyplot as plt
import matplotlib.colors
import matplotlib.cm
import geoplot.crs as gcrs


norm = matplotlib.colors.LogNorm(vmin=10, vmax=99999)


plt.ioff()

shape_bundeslaender = gpd.read_file('./gadm36_DEU_shp/gadm36_DEU_1.shp')

projection = gcrs.LambertAzimuthalEqualArea()

for d in pd.date_range(min(df['Meldedatum']), max(df['Meldedatum']))[13:]:
    heute = N.loc[:,d]

    fig = plt.figure(1, figsize=(7,3.5))
    fig.clf()
    
    fig.patch.set_facecolor('gray')
    
    
    fig.suptitle(d.strftime('%Y-%m-%d'),x=0,ha='left')
    ax1 = plt.subplot(121, projection = projection)
    ax1.set_title("N")
    ax1 = gplt.choropleth(f.join(heute, on='IdLandkreis'),hue='N',norm=norm, cmap='viridis_r',ax=ax1)
    ax1 = gplt.polyplot(shape_bundeslaender,ax=ax1,edgecolor='white',linewidth=0.25,zorder=2)
    
    
    
    ax2 = plt.subplot(122, projection = projection)
    ax2.set_title("N * area")
    ax2 = gplt.choropleth(f.join(N_area.loc[:,d], on='IdLandkreis'),hue='N',norm=norm, cmap='viridis_r',ax=ax2)
    ax2 = gplt.polyplot(shape_bundeslaender,ax=ax2,edgecolor='white',linewidth=0.25,zorder=2)

    fig.subplots_adjust(left=0.01, right=0.99, top=0.99, bottom=0.01,wspace=0,hspace=0)
    
    plt.colorbar(matplotlib.cm.ScalarMappable(norm=norm, cmap='viridis_r'),ax = [ax1,ax2], shrink=0.75)  

    fig.savefig('output/' + d.strftime('%Y-%m-%d') + '.png',dpi=200)
    
    
    print("Done",d)