Plots COVID risk graph, probability of infected person present given infection rate in population and number of persons present.

plotCovidRiskGraph.py

# Antti Lipponen, 18 December 2021
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import zoom
import seaborn as sns
from matplotlib import rcParams
from matplotlib.colors import Normalize
import matplotlib as mpl
from matplotlib.colors import ListedColormap
import matplotlib.patheffects as PathEffects


rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['Lato']
rcParams['hatch.linewidth'] = 8
rcParams['hatch.color'] = '#00000020'


def compute_prob(incidence, N):
    pHasCovid = incidence / 100_000
    pNoCovid = 1 - pHasCovid
    pNoCovidForGroup = pNoCovid**N
    pCovidForGroup = 1 - pNoCovidForGroup
    return pCovidForGroup * 100


def computeisoline(p, N):
    pCovidForGroup = p / 100
    pNoCovidForGroup = 1 - pCovidForGroup
    pNoCovid = pNoCovidForGroup**(1 / N)
    pHasCovid = 1 - pNoCovid
    incidence = pHasCovid * 100_000
    return incidence


xlim = [1, 250]
ylim = [50, 10000]

Xlocs = np.linspace(1, 250, 100).astype(int)
Ylocs = np.logspace(np.log10(40), np.log10(10000), 100).astype(int)
Xgrid, Ygrid = np.meshgrid(Xlocs, Ylocs)
Zgrid = np.nan * np.ones_like(Xgrid)
for ii in range(len(Xlocs)):
    for jj in range(len(Ylocs)):
        Zgrid[jj, ii] = compute_prob(Ylocs[jj], Xlocs[ii])
Xgrid = zoom(Xgrid, 2)
Ygrid = zoom(Ygrid, 2)
Zgrid = zoom(Zgrid, 2)

colors = sns.color_palette("rocket_r", n_colors=20)
for ii in range(4, 20, 2):
    meancolor = np.array(colors[ii:ii + 2]).mean(axis=0)
    colors[ii] = meancolor
    colors[ii + 1] = meancolor

cmap = ListedColormap(colors.as_hex())
norm = Normalize(vmin=0, vmax=100)
levels = np.linspace(0, 100, 500) + 0.01

fig = plt.figure(figsize=(3, 4), dpi=100)
ax = fig.add_axes([0.2, 0.26, 0.75, 0.71])
ax.contourf(Xgrid, Ygrid, Zgrid, levels=levels, cmap=cmap)
for p in [5, 10, 15, 20, 30, 40, 50, 60, 70, 80, 90]:
    Nxx = np.linspace(xlim[0], xlim[1], 250)
    yy = computeisoline(p, Nxx)
    ax.plot(Nxx, yy, 'w-', linewidth=0.3, zorder=99)
ax.grid(True, alpha=0.9, c='#a0a0a0', linestyle='--', linewidth=0.2)


p = 50
xx = np.linspace(6.5, 255, 260)
ax.fill_between(xx, y1=computeisoline(p, xx), y2=11000, hatch='/', zorder=99, fc='#dedede', clip_on=False, alpha=0.25)

ax.plot([6.5 + 0.25, 6.5, 255, 255, 250], [10000, 11000, 11000, computeisoline(p, xx)[-1], computeisoline(p, xx)[-1] + 5], linewidth=0.5, c='#000000', alpha=0.2, clip_on=False)
txt = ax.text(264, 80, 'More likely to have a person with COVID-19 present than none with COVID-19.', fontsize=5, c='#303030', ha='right', va='bottom', alpha=1.0, rotation=-90, zorder=100, clip_on=False)

ax.set_xlim(xlim[0], xlim[1])
ax.set_ylim(ylim[0], ylim[1])
ax.set_yscale('log')
ax.set_xlabel('Number of persons present')
ax.set_ylabel('Infected persons / 100,000 persons in the region', labelpad=-1)

ax.set_yticks([100, 1000, 10000])
ax.set_yticklabels(['100', '1,000', '10,000'], rotation=65)

ax.set_xticks([2, 50, 100, 150, 200, 250])
ax.set_xticklabels(['2', '50', '100', '150', '200', '250'], rotation=0)

axcb = fig.add_axes([0.05, 0.11, 0.9, 0.02])
cb = fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmap), cax=axcb, orientation='horizontal')
cb.set_label(label='Probability for having at least 1 infected person present (%)', fontsize=7.5, labelpad=1.5)
cb.ax.tick_params(labelsize=8)

ax.text(0.01, 0.005, 'Computed with basic probability theory. Assuming each person is equally likely to be infected.', fontsize=4.5, transform=fig.transFigure, ha='left', va='bottom', alpha=0.35)
ax.text(0.99, 0.005, '@anttilip', fontsize=4.5, transform=fig.transFigure, ha='right', va='bottom', alpha=0.25)

txt = ax.text(35, 53.5, '<5%', fontsize=8, c='#303030', ha='center', va='bottom', alpha=0.9, zorder=100)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
txt = ax.text(110, 51.5, '5-10%', fontsize=8, c='#303030', ha='center', va='bottom', alpha=0.9, zorder=100, rotation=-30)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
txt = ax.text(220, 49.5, '10-15%', fontsize=8, c='#303030', ha='center', va='bottom', alpha=0.9, zorder=100, rotation=-13.5)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
txt = ax.text(225, 69, '15-20%', fontsize=8, c='#303030', ha='center', va='bottom', alpha=0.9, zorder=100, rotation=-13.9)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
txt = ax.text(247.5, 123, '20-30%', fontsize=8, c='#303030', ha='right', va='center', alpha=0.9, zorder=100, rotation=-13.9)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
txt = ax.text(247.5, 188, '30-40%', fontsize=8, c='#303030', ha='right', va='center', alpha=0.9, zorder=100, rotation=-14)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
txt = ax.text(247.5, 260, '40-50%', fontsize=8, c='#303030', ha='right', va='center', alpha=0.9, zorder=100, rotation=-14.4)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
txt = ax.text(247.5, 350, '50-60%', fontsize=8, c='#303030', ha='right', va='center', alpha=0.9, zorder=100, rotation=-14.5)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
txt = ax.text(247.5, 460, '60-70%', fontsize=8, c='#303030', ha='right', va='center', alpha=0.9, zorder=100, rotation=-14.5)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
txt = ax.text(247.5, 600, '70-80%', fontsize=8, c='#303030', ha='right', va='center', alpha=0.9, zorder=100, rotation=-15)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
txt = ax.text(247.5, 835, '80-90%', fontsize=8, c='#303030', ha='right', va='center', alpha=0.9, zorder=100, rotation=-16)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
txt = ax.text(247.5, 2850, '>90%', fontsize=8, c='#303030', ha='right', va='center', alpha=0.9, zorder=100)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])

plt.savefig("graph.png", dpi=900)