Rates of Covid19 in Madrid
Initial exploration of infection data
- Covid19 Weekly Data from the 'Comunidad de Madrid'
- Covid 19-TIA Zonas Básicas de Salud (CSV)
- Histograms of number of zones in each category in the last 6 weeks
- Maps for the most recent week
- Treating the rate as a surface and plotting contours
- Confusion matrix for rate > 1000 and restrictions
- Weak relationship between population density and rate
- Comparison of individual zones to the best and worst
!pip install geopandas
import pandas as pd
import numpy as np
import re
import geopandas
import plotnine
from plotnine import *
import matplotlib.pyplot as plt
import matplotlib.tri as tri
from mpl_toolkits.axes_grid1 import make_axes_locatable
from sklearn.metrics import confusion_matrix
%matplotlib inline
Covid19 Weekly Data from the 'Comunidad de Madrid'
https://datos.comunidad.madrid/catalogo/dataset/covid19_tia_zonas_basicas_salud
Covid 19-TIA Zonas Básicas de Salud (CSV)
On a Tuesday morning, weekly data are provided for the 286 basic health areas:
- confirmed cases and rate of accumulated incidence for basic health areas
- the source of data is the Madrid region's network of epidemiological vigilence
- the data are continually updated; you should always show the date when the data were accessed
- the polygon data includes the population figures from the region's statistical institute from the population register of 1 January 2019.
"Información epidemiológica Covid-19: Casos confirmados y tasa de incidencia acumulada (TIA) por Zonas Básicas de Salud (Decreto 52/2010, de 29 de julio, BOCM 9 de agosto de 2010). Informe diario o semanal con datos de casos confirmados, casos confirmados con infección activa y tasas de incidencia acumulada de los últimos 14 días y desde el inicio de la epidemia (25 de febrero de 2020). La fuente de los datos es la Red de Vigilancia Epidemiológica de la Comunidad de Madrid. Cuando se utilicen los datos se debe indicar en qué fecha se ha accedido a los mismos, dada su actualización continua. Hasta el 1 de julio de 2020 los registros de datos correspondientes a cada fecha de informe se han añadido diariamente. A partir del 2 de julio de 2020 la actualización pasa a ser semanal. A fecha 1 de junio, se ha procedido a actualizar el conjunto de datos de Zonas Básicas de Salud, con los datos actualizados de geometrías proporcionados por el Instituto de Estadística de la Comunidad de Madrid. Las nuevas geometrías incluyen los datos poblacionales actualizados por el Instituto de Estadística de la Comunidad de Madrid, relativos al Padrón Continuo a 1 de enero de 2019."
Here we will read the CSV data dated 29 September 2020 and explore:
zones_df=pd.read_csv('covid19_tia_zonas_basicas_salud_s.csv', delimiter=';', encoding='latin')
print ('Columns:')
[print(col) for col in list(zones_df.columns)]
print()
print('Shape: ',zones_df.shape)
print()
print('Dates: ',list(zones_df.fecha_informe.unique()))
Tidy-up the dataframe by
- replacing decimal
,with decimal. - shortening key column names
- dropping time from the date column
for col in ['tasa_incidencia_acumulada_activos_ultimos_14dias', 'tasa_incidencia_acumulada_ultimos_14dias', 'tasa_incidencia_acumulada_total']:
zones_df[col] = pd.to_numeric(zones_df[col].apply(lambda x: re.sub(',', '.', x)))
zones_df.rename(columns={'zona_basica_salud':'zona_basic','tasa_incidencia_acumulada_ultimos_14dias':'tasa', 'casos_confirmados_ultimos_14dias':'casos'}, inplace=True)
zones_df.fecha_informe=zones_df.fecha_informe.apply(lambda x: x[5:10])
Data for the centre of Madrid (where the main office of the regional authority are located).
zones_df[zones_df.zona_basic=='Cortes'].T
fig, ax = plt.subplots(1,1,figsize=(10,5))
plt.plot(zones_df[zones_df.zona_basic=='Cortes'].fecha_informe[::-1], zones_df[zones_df.zona_basic=='Cortes'].casos[::-1])
plt.title('Cases confirmed in the last 14 days - Cortes');
fig, ax = plt.subplots(1,1,figsize=(10,5))
plt.plot(zones_df[zones_df.zona_basic=='Cortes'].fecha_informe[::-1], zones_df[zones_df.zona_basic=='Cortes'].tasa[::-1])
plt.ylim(0,1000)
plt.title('Rate of accumulated incidence for 100 000 people in the last 14 days - Cortes');
print(len(zones_df.zona_basic.unique()))
print(zones_df.zona_basic.unique()[:10])
Allocate to bins based on rate
zones_df.tasa.fillna(0, inplace=True)
cut_labels=['<200','200-400','400-600','600-800','800-1000','>1000']
cut_bins = [-1., 200., 400., 600., 800., 1000., max(zones_df.tasa)]
zones_df['tasa_bin'] = pd.cut(zones_df.tasa, bins=cut_bins, labels=cut_labels)
zones_df.tasa_bin.value_counts()
Helper function to check for errors in the presentation of the names, e.g. spaces, 'de',...
def check_zones(zonas_restricted):
lst=[]
for zona_basic in zones_df.zona_basic[:286]:
if zona_basic in zonas_restricted: lst.append(zona_basic)
print(len(lst), set(zonas_restricted)-set(lst))
Identify the initial 37 areas restricted from 21 September 2020:
zonas_restricted_1=["Puerta Bonita","Vista Alegre","Guayaba", "Almendrales", "Las Calesas", "Zofío", "Orcasur", "San Fermín", "San Andrés", "San Cristóbal",
"El Espinillo", "Los Rosales", "Villa Vallecas", "Entrevías", "Martínez de la Riva", "San Diego",
"Numancia", "Peña Prieta", "Pozo de Tío Raimundo", "Ángela Uriarte", "Alcalá de Guadaira", "Federica Montseny", "Doctor Cirajas",
"Gandhi", "Daroca", "La Elipa", "Alicante", "Cuzco", "Francia", "Humanes de Madrid", "San Blas", "Isabel II",
"Las Margaritas", "Sánchez Morate", "Reyes Católicos", "Alcobendas - Chopera", "Miraflores"]
check_zones(zonas_restricted_1)
Identify the additional 8 areas restricted from 28 September 2020:
zonas_restricted_2=["Panaderas", "Doctor Trueta", "Miguel Servet", "Campo de la Paloma", "Rafael Alberti", "García Noblejas", "Vicálvaro - Artilleros","Orcasitas"]
check_zones(zonas_restricted_2)
Label zones with restricted group
zones_df['restricted']='0'
for i in range(len(zonas_restricted_1)): zones_df.loc[zones_df.zona_basic==zonas_restricted_1[i], 'restricted']='1'
for i in range(len(zonas_restricted_2)): zones_df.loc[zones_df.zona_basic==zonas_restricted_2[i], 'restricted']='2'
zones_df.restricted.value_counts()/len(zones_df.fecha_informe.unique())
Helper function to reorder the presentation of the dates
def col_func(fecha): return fecha[3]+fecha[4]+fecha[2]+fecha[0]+fecha[1]
plotnine.options.figure_size = (14, 8)
ggplot(zones_df[zones_df.fecha_informe>'08/24'], aes(x='tasa_bin', fill='restricted')) \
+ geom_histogram(binwidth=1, alpha=0.6, position='stack') \
+ facet_wrap('fecha_informe', labeller=labeller(cols=col_func)) \
+ theme_minimal() \
+ labs(title="Zonas básicas de salud",
x='Tasa incidencia acumulada ultimos 14 días',
y="Número de zonas") \
+ scale_fill_brewer(palette = 'Reds')
week_df=zones_df[zones_df.fecha_informe==zones_df.fecha_informe.unique().max()][['zona_basic','tasa','casos','tasa_bin','restricted']]
df=(geopandas.read_file('./maps/zonas_basicas_salud.shp')).merge(week_df)
df['geometry_pt']=df.geometry.centroid
df.head().T
Maps of rate of accumulated instance and restricted zones
fig, (ax1, ax2) = plt.subplots(ncols=2, sharex=True, sharey=True, figsize=(12,6))
df.plot(ax=ax1, column=df.tasa_bin, cmap='Reds', legend=True)
df.plot(ax=ax1, color='white', edgecolor='grey', alpha=0.1)
ax1.set_title('Rate of accumulated incidence')
ax1.axis('off')
df.plot(ax=ax2, column=df.restricted, cmap='Reds', legend=True)
df.plot(ax=ax2, color='white', edgecolor='grey', alpha=0.1)
ax2.set_title('Restricted zones')
ax2.axis('off');
df.tasa_bin.value_counts()
df.restricted.value_counts()
def getXY(pt): return (pt.x, pt.y)
def data_for_week(fecha):
week_df=zones_df[zones_df.fecha_informe==fecha][['zona_basic','tasa','casos','tasa_bin','restricted']]
df=geopandas.read_file('./maps/zonas_basicas_salud.shp').merge(week_df)
x, y = [list(t) for t in zip(*map(getXY, df.geometry.centroid))]
z = df.tasa
return x,y,z
def multi_plot(i):
fecha=weeks[i]
x,y,z = data_for_week(fecha)
fig, ax = plt.subplots(1,1,figsize=(10,8))
ax.set(xlim=(360000, 500000), ylim=(4420000, 4550000))
cntr = ax.tricontourf(x, y, z, [0,200,400,800,1000,2000], cmap="Reds") # plot contour fill
df.plot(color='white', edgecolor='grey', alpha=0.1, ax=ax)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.0)
fig.colorbar(cntr, drawedges=False, cax=cax)
ax.set_title(f'Tasa incidencia acumulada en los últimos 14 días\n fecha {fecha[3]+fecha[4]+fecha[2]+fecha[0]+fecha[1]}, maximum {str(int(max(z)))}/100000 habitantes')
ax.axis('off')
weeks = zones_df.fecha_informe.unique()[:8][::-1]
for i in range(len(weeks)): multi_plot(i)
Plots for 4 weeks
def plot_4_weeks(plot_weeks):
fig, axs = plt.subplots(1, 4, figsize=(16, 4))
for i in range(4):
fecha=plot_weeks[i]
x,y,z = data_for_week(plot_weeks[i])
df.plot(color='white', edgecolor='grey', alpha=0.1, ax=axs[i])
axs[i].set_title(f'\nfecha {fecha[3]+fecha[4]+fecha[2]+fecha[0]+fecha[1]}\nmax. {str(int(max(z)))}/100000 hab.', fontsize=8)
axs[i].axis('off')
cntr=axs[i].tricontourf(x,y,z,[0,200,400,800,1000,2000], cmap="Reds")
divider = make_axes_locatable(axs[-1])
cax = divider.append_axes("right", size="3%", pad=0.0)
fig.colorbar(cntr, drawedges=False, cax=cax)
fig.suptitle('Tasa incidencia acumulada en los últimos 14 días', fontsize=12, va='baseline')
plt.savefig('four_weeks.png', bbox_inches='tight');
plot_4_weeks(weeks[:4])
plot_4_weeks(weeks[4:])
week_df=zones_df[zones_df.fecha_informe==zones_df.fecha_informe.max()][['zona_basic','tasa','casos','tasa_bin','restricted']]
df=(geopandas.read_file('./maps/zonas_basicas_salud.shp')).merge(week_df)
print(confusion_matrix((df.restricted != '0')*1, (df.tasa>1000)*1))
- 228 zones with rate less than 1000 and no restrictions.
- 36 zones with rate over 1000 and restrictions.
13 zones with rate over 1000 but no restrictions
df[(df.restricted == '0') & (df.tasa>1000)][['zona_basic','tasa']].sort_values('tasa', ascending=False)
9 zones with rate under 1000 but restricted
df[(df.restricted != '0') & (df.tasa<1000)][['zona_basic','tasa']].sort_values('tasa')
df['area_sq_km'] = df.geometry.area/1_000_000
df['pob_densidad_sq_km'] = df.pob_pad19/df.area_sq_km
df[['tasa','pob_densidad_sq_km']].corr()
ggplot(df, aes(x='pob_densidad_sq_km',y='tasa',color='restricted')) \
+ geom_point() \
+ geom_smooth(aes(x='pob_densidad_sq_km',y='tasa'), color='black',method = "lm")
df_week=zones_df[zones_df.fecha_informe==zones_df.fecha_informe.max()]
most = df_week[df_week.tasa==df_week.tasa.max()].zona_basic.to_list()[0]
least = df_week[df_week.tasa==df_week.tasa.min()].zona_basic.to_list()[0]
most, least
zona = 'Cortes'
df_tmp = zones_df.pivot_table(values='tasa', index='fecha_informe', columns='zona_basic').reset_index()
zbs = [most, zona, least]
title= zbs[1] + ' compared to least (' + zbs[2] + ') and most (' + zbs[0] + ')'
ggplot(df_tmp, aes(x='fecha_informe',y=zbs[1],group=1)) \
+ geom_point(colour='blue', alpha = 0.5) \
+ geom_line(aes(y=zbs[1]), colour='blue') \
+ geom_smooth(aes(y=zbs[0]), colour='red' ) \
+ geom_smooth(aes(y=zbs[2]), colour='green') \
+ labs(title=title,
x='Fecha informe',
y="Tasa") \
+ geom_text(data=df_tmp[df_tmp.fecha_informe=='09/22'],label=zbs[1], color='blue', hjust= 'right', vjust='bottom') \
+ theme_minimal()
