<jupyter_start><jupyter_text>Choropleth map of COVID-19 cases in Spain using GeopandasThis tutorial shows how to display official data from Spanish government on map of Spain using Geopandas package. To do that we need to get COVID-19 data, map of regions of Spain, combine them and display. This approach may be used for visual representation of modelling results (e.g. SIR, SEIR models as mentioned [here](https://github.com/DataForScience/Epidemiology101?fbclid=IwAR2mhfvS6XvPm1JUgdUnDKMmgkIJE65_IfG4pMoquh_9kMwavFh59FHJYCQ)) Import libraries we will be using<jupyter_code>import matplotlib.pyplot as plt
import geopandas as gpd
import pandas as pd
import numpy as np
from pathlib import Path
import datetime
plt.rcParams['figure.figsize'] = [10, 5]
%matplotlib inline<jupyter_output><empty_output><jupyter_text>1. Shape files dataShape files define geographical regions and would be used for visualization.To download shape files clik on the [link](https://www.diva-gis.org/gdata) and:1) Select from drop-down list `Country`: `Spain`2) Select from drop-down list `Subject`: `Administrative areas`3) Click button `Ok`4) At new page click `Download` link and then save `ESP_adm.zip` archive (7.8 MB) to folder with current notebook.5) Unzip `ESP_adm.zip` to `ESP_adm` folder in the current folder (may be done by uncommenting and running the following line of code)A `.shp` file is actually just one of the 3 files that are necessary for your shape data to render — you may have noticed that a `.shp` file comes in a zip file when you download. This is because you also need a `.dbf` and a `.shx` file. While the `.shp` files gives figures their geometry, `.dbf` stores attribute data and object ids, and `.shx` provides an index that is helpful with certain software and packages [(ref)](https://medium.com/thoughtful-data-science/geopandas-an-introduction-c544a352c662).<jupyter_code>#!unzip ESP_adm.zip -d ESP_adm<jupyter_output><empty_output><jupyter_text>Dictionary of standard region names shortcuts. We will use it later during merging `.shp` and tabular COVID-19 data.<jupyter_code>codes_dictionary = {
    'Andalucía': 'AN',
    'Aragón': 'AR',
    'Cantabria': 'CB',
    'Castilla-La Mancha': 'CM',
    'Castilla y León': 'CL',
    'Cataluña': 'CT',
    'Ceuta y Melilla': 'CE_ML',
    'Comunidad de Madrid': 'MD',
    'Comunidad Foral de Navarra': 'NC',
    'Comunidad Valenciana': 'VC',
    'Extremadura': 'EX',
    'Galicia': 'GA',
    'Islas Baleares': 'IB',
    'Islas Canarias': 'CN',
    'La Rioja': 'RI',
    'País Vasco': 'PV',
    'Principado de Asturias': 'AS',
    'Región de Murcia': 'MC'
}
path = Path('ESP_adm/ESP_adm1.shp')
map_df = gpd.read_file(path, encoding='latin9')
map_df.head()
type(map_df)<jupyter_output><empty_output><jupyter_text>We see above that `map_df` is not usual `Pandas.DataFrame`, but `geopandas.geodataframe.GeoDataFrame` object.<jupyter_code>map_df['NAME_1'] = codes_dictionary.keys()
len(map_df['NAME_1'].values), map_df.NAME_1.values<jupyter_output><empty_output><jupyter_text>So 18 regions of Spain are represented at our `.shp` file (actually there are 17 regions + 2 autonomous cities in Spain).<jupyter_code>map_df.columns<jupyter_output><empty_output><jupyter_text>Geopandas supports easy to use `plot` function:<jupyter_code>map_df.plot();<jupyter_output><empty_output><jupyter_text>Plot separately Spain without Islas Canarias:<jupyter_code>map_df[map_df['NAME_1'] != 'Islas Canarias'].plot();<jupyter_output><empty_output><jupyter_text>Islas Canarias:<jupyter_code>map_df[map_df['NAME_1'] == 'Islas Canarias'].plot();<jupyter_output><empty_output><jupyter_text>2. COVID-19 dataDownload official COVID-19 data provided by Spanish government by region/date by running following cell:<jupyter_code>import requests
covid_url = 'https://covid19.isciii.es/resources/serie_historica_acumulados.csv'
response = requests.get(covid_url)
with open("spain.csv", 'wb') as file:
    file.write(response.content)<jupyter_output><empty_output><jupyter_text>Read COVID-19 data 'Magic number' `-8` reflects that so far 8 last strings are natural language notes, not a data. If you run code and receive value starting from 'NOTE1:' at `CCAA` column - please change `-8` to `-9`, `-10` until you would not read those notes into DataFrame.`NaN` values replaced by zeros, because they they testify lack of known values at the moment.<jupyter_code>df_covid = pd.read_csv('spain.csv', encoding='latin1')[:-8].fillna(0)
ccaa_num = len(np.unique(df_covid['CCAA']))
df_covid.tail(ccaa_num)<jupyter_output><empty_output><jupyter_text>We are interested in number of cases per region (`CASOS` column), but this column has `NaN` values. Explanation may be found in [data description (in Spanish)](https://rubenfcasal.github.io/COVID-19/):>On 2020-04-25 two new PCR + and TestAC + variables are included with the number of positive results in PCR tests and in antibody tests (it is assumed that with symptoms) and NA is reported in some values of the variable Cases from 2020 -04-15. Although in previous versions it was assumed that only positives were considered confirmed by PCR, apparently this is not the case, so these values were replaced by the sum of PCR + and TestAC + (since 2020-04-28 Galicia already reports these two variables).According to description above, `CASOS` = `PCR+` + `TestAc+` (however it may be the place of mistake - what if some individuals were taking both PCR and TestAc tests?):<jupyter_code>df_covid['CASOS'] = df_covid['PCR+'] +  df_covid['TestAc+']
df_covid.tail(ccaa_num)
ccaa_num<jupyter_output><empty_output><jupyter_text>Well, we have covid-19 data for 17+2 regions of Spain, but only 18 regions of Spain are represented in our .shp file. Via eyeballing we found out that `CE` (Ceuta) and `ML` (Melilla) regions should be merged.<jupyter_code>df_tmp = df_covid[df_covid['CCAA'].isin(['CE', 'ML'])]
df_tmp.index = range(len(df_tmp))
df_tmp
assert len(df_tmp) / 2 == len(df_covid) / ccaa_num<jupyter_output><empty_output><jupyter_text>Create dataframe for `CE_ML` ('Ceuta y Melilla') rows<jupyter_code>df_covid_ce_ml = pd.DataFrame(columns=df_covid.columns)
for idx in range(len(df_tmp) // 2):
    series_ce = df_tmp.iloc[2 * idx]
    series_ml = df_tmp.iloc[2 * idx + 1]
    series_ce_ml = pd.Series({'CCAA': 'CE_ML',
                              'FECHA': series_ce.FECHA,
                              'CASOS': series_ce.CASOS + series_ml.CASOS,
                              'PCR+': series_ce['PCR+'] + series_ml['PCR+'],
                              'TestAc+': series_ce['TestAc+'] + series_ml['TestAc+'],
                              'Hospitalizados': series_ce.Hospitalizados + series_ml.Hospitalizados,
                              'UCI': series_ce.UCI + series_ml.UCI,
                              'Fallecidos': series_ce.Fallecidos + series_ml.Fallecidos,
                              'Recuperados': series_ce.Recuperados + series_ml.Recuperados})
    df_covid_ce_ml = df_covid_ce_ml.append(series_ce_ml, ignore_index=True)
df_covid_ce_ml<jupyter_output><empty_output><jupyter_text>Drop `CE` and `ML` rows from `df_covid` dataframe.<jupyter_code>df_covid_dropped = df_covid[df_covid['CCAA'].isin(
    set(np.unique(df_covid['CCAA'])).difference(set(['CE', 'ML'])))]
df_covid_dropped<jupyter_output><empty_output><jupyter_text>Add instead combined `CE_ML` rows<jupyter_code>df_combined = df_covid_dropped.append(df_covid_ce_ml, ignore_index=True)
df_combined['FECHA'] = pd.to_datetime(df_combined['FECHA'], dayfirst=True)
df_combined = df_combined.sort_values(by=['FECHA'], ignore_index=True)
df_combined.tail(18)<jupyter_output><empty_output><jupyter_text>3. Comparison of active cases numbers by regionLet's visualize number of active cases by region. At firs we define function that extracts data by required region from all of our COVID-19 data.<jupyter_code>def get_region_df(df, region_code):
    """
    Returns region dataframe from whole country dataframe,
    computes 'CASOS_ACTIVOS' column.
    """
    df = df_combined[df_combined['CCAA'] == region_code]
    df = df[['FECHA','CASOS', 'Fallecidos', 'Recuperados']].set_index('FECHA')
    df['CASOS_ACTIVOS'] = df['CASOS'] - df['Fallecidos'] - df['Recuperados']
    return df<jupyter_output><empty_output><jupyter_text>Let us get data from Madrid and show 5 last days, for example:<jupyter_code>df_md = get_region_df(df_combined, 'MD')
df_md.tail()<jupyter_output><empty_output><jupyter_text>Following is graph of number of active COVID-19 cases in Madrid, Cataluña and Castilla-La Mancha regions.<jupyter_code>df_md = get_region_df(df_combined, 'MD')
df_ct = get_region_df(df_combined, 'CT')
df_cm = get_region_df(df_combined, 'CM')
plt.plot(df_md['CASOS_ACTIVOS'], label='Madrid');
plt.plot(df_ct['CASOS_ACTIVOS'], label='Cataluña');
plt.plot(df_cm['CASOS_ACTIVOS'], label='Castilla-La Mancha');
plt.title('Covid-19 active cases')
plt.legend();<jupyter_output><empty_output><jupyter_text>4. Visualization of the last day of observations with reference to the mapIn our case - May 17, 2020.<jupyter_code>def long2short(long_code):
    return codes_dictionary[long_code]<jupyter_output><empty_output><jupyter_text>Grab the last 18 rows of the COVID-19 dataframe. They represent information over regions of Spain at the last observation day.<jupyter_code>df_covid_day = df_combined.iloc[-18:]
df_covid_day
type(map_df)<jupyter_output><empty_output><jupyter_text>Now add to GeoDataFrame `map_df` column `cases` - values that we would visualize.<jupyter_code>map_df['cases'] = pd.Series([df_covid_day[df_covid_day['CCAA'] == long2short(code)]['CASOS'].iloc[0]
                                      for code in map_df['NAME_1']])
map_df.head()<jupyter_output><empty_output><jupyter_text>Following commented cell introduces other labels that also may be visualized via choropleth map. If interested in those values - uncomment following cell, and during visualization change value of `var` to corresponding value, e. g. `hospitalized`, `dead` or `recovered`.<jupyter_code># map_df['hospitalized'] = pd.Series([df_covid_day[df_covid_day['CCAA'] == long2short(code)]['Hospitalizados']
#                                       for code in map_df['NAME_1']])
# map_df['dead'] = pd.Series([df_covid_day[df_covid_day['CCAA'] == long2short(code)]['Fallecidos'].iloc[0]
#                                       for code in map_df['NAME_1']])
# map_df['recovered'] = pd.Series([df_covid_day[df_covid_day['CCAA'] == long2short(code)]['Recuperados'].iloc[0]
#                                       for code in map_df['NAME_1']])
df_mainland = map_df[map_df['NAME_1'] != 'Islas Canarias']
df_canarias = map_df[map_df['NAME_1'] == 'Islas Canarias']
# variable to display
var = 'cases'
# set the range for choropleth
vmin, vmax = df_mainland['cases'].min(), df_mainland['cases'].max()
# create figure and axes for Matplotlib
fig, ax = plt.subplots(1, figsize=(20, 12))
# create map
# df_mainland.plot(column=var, cmap='Blues', linewidth=0.8, ax=ax, edgecolor='0.8')
df_mainland.plot(column=var, cmap='Blues', ax=ax)
ax.set_title(f"{var}", fontdict={'fontsize': '25', 'fontweight': '3'});
# create colorbar as legend
sm = plt.cm.ScalarMappable(cmap='Blues' ,norm=plt.Normalize(vmin=vmin, vmax=vmax))

# empty array for the data range
sm._A = []

# add the colorbar to the figure
cbar = fig.colorbar(sm)

# save to file
# plt.savefig('cases_spain_may_17.png')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12210729312762973, 'fraction_numerical': 0.01393758765778401, 'mean_word_length': 5.107601713062098, 'pattern_counts': {'":': 0, '<': 89, '<?xml version=': 0, '>': 90, 'https://': 5, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '19891', 'license': 'Apache-2.0', 'n_tokens_mistral': 4113, 'n_tokens_neox': 3896, 'n_words': 1028, 'path': '_notebooks/2020-05-20-Spain-covid19-map.ipynb', 'repo_name': 'CalmScout/blog', 'size': '11408'}