<jupyter_start><jupyter_text>Importing CSSEGI data into communities[Index](0-index.ipynb)In this notebook, we assign CSSEGI data to the communities constructed in [1-clustering](1-clustering.ipynb). Data reported in the CSSEGI dataset is assigned to a community based on the latitude and longitude. Data can be accessed [here](https://github.com/CSSEGISandData/COVID-19).  Imports and global variables<jupyter_code>%load_ext autoreload
%autoreload 2
%matplotlib inline
from pathlib import Path
import sys
import numpy as np
import pandas as pd
import geopandas as gpd
import datetime
import json

import matplotlib.pyplot as plt
import matplotlib.colors as mco
import matplotlib.gridspec as mgs
import matplotlib.ticker as ticker
import matplotlib.cm as cm
from matplotlib import animation
plt.rcParams['svg.fonttype'] = 'none'

from IPython.display import HTML
from IPython.display import Image
sys.path.append('../code')
from functions import plot_omega_profile, plot_omega_map
datadir = Path('../data')
if not datadir.is_dir():
    raise ValueError("Data dir doesn'nt exist!")
resdir = Path('../results/')
if not resdir.is_dir():
    raise ValueError('No results directory!')
complevel=7
complib='zlib'

resfile = resdir / 'safegraph_analysis.hdf5'
# resfile = resdir / 'safegraph_analysis_monthly.hdf5'
with pd.HDFStore(resfile, complevel=complevel, complib=complib) as store:
    print(f"File {resfile.stem} has {len(store.keys())} entries.")

cssefile = resdir / 'cssegi_sanddata_analysis.hdf5'
with pd.HDFStore(cssefile, complevel=complevel, complib=complib) as store:
    print(f"File {cssefile.stem} has {len(store.keys())} entries.")

tfmt = '%Y-%m-%d'
exts = ['.png', '.svg']<jupyter_output><empty_output><jupyter_text>Load clustered CBGs<jupyter_code>key = "/clustering/clusters"
with pd.HDFStore(resfile, complevel=complevel, complib=complib) as store:
    clusters = store[key]
clusters
clusters_pos = clusters.loc[:, ['X','Y']].rename(columns={'X': 'longitude', 'Y': 'latitude'})
clusters_pos<jupyter_output><empty_output><jupyter_text>Load data from geometry<jupyter_code>geofile = datadir / 'safegraph_open_census_data'  / 'geometry' / 'cbg.geojson'
if not geofile.is_file():
    raise ValueError("Geo file doesn't exist!")
geo = gpd.read_file(geofile).astype({'CensusBlockGroup': 'int64'})
geo.set_index('CensusBlockGroup', inplace=True)
XY = geo.representative_point().apply(lambda X: [X.x, X.y]).tolist()
df_xy = pd.DataFrame(data=np.array(XY), index=geo.index, columns=['longitude','latitude'])
df_xy<jupyter_output><empty_output><jupyter_text>Load the confirmed cases from CSSEGI data<jupyter_code>dfile = datadir / 'time_series_covid19_confirmed_US.csv'
df_confirmed_us = pd.read_csv(dfile)
df_confirmed_us = df_confirmed_us.loc[df_confirmed_us['iso2'] == 'US'] # only keep the US
#.astype({'FIPS': 'int64'})
df_confirmed_us<jupyter_output><empty_output><jupyter_text>Drop the locations with no latitude or longitude information<jupyter_code>idx = ( (df_confirmed_us['Lat'] == 0.) | (df_confirmed_us['Long_'] == 0.) )
df_confirmed_us.drop(index=df_confirmed_us.index[idx], inplace=True)
df_confirmed_us<jupyter_output><empty_output><jupyter_text>Mapping clusters to CSSEGI FIPS based on (Longitude, Latitude)<jupyter_code>clusters_pos.sort_index(inplace=True)
XY_clusters = clusters_pos.loc[:, ['longitude', 'latitude']].to_numpy()
XY_csse = df_confirmed_us.loc[:, ['Long_', 'Lat']].to_numpy()
# iterate over CSSEGI data
cluster_ids = []
for i in range(len(df_confirmed_us.index)):
    xy = XY_csse[i]
    imin = np.argmin(np.linalg.norm(XY_clusters - xy, axis=1))
    cluster_ids.append(clusters_pos.index[imin])
col = 'cluster_id'
if col in df_confirmed_us.columns:
    del df_confirmed_us[col]

df_confirmed_us.insert(1, col, cluster_ids)
df_confirmed_us<jupyter_output><empty_output><jupyter_text>Retain only the dates and cluster_ID<jupyter_code>dates = []
for c in df_confirmed_us.columns:
    try:
        t = datetime.datetime.strptime(c, '%m/%d/%y')
        dates.append(c)
    except ValueError:
        continue
columns = ['cluster_id'] + dates

df_clusters = df_confirmed_us.loc[:, columns].groupby('cluster_id').apply(sum)
del df_clusters['cluster_id']
df_clusters
clusters_csse = pd.DataFrame(data=np.zeros((len(clusters.index),len(df_clusters.columns)), np.uint), index=clusters.index, columns=df_clusters.columns)
clusters_csse.loc[df_clusters.index] = df_clusters
clusters_csse = clusters_csse.T
clusters_csse.index = [datetime.datetime.strptime(s, '%m/%d/%y') for s in clusters_csse.index]
clusters_csse<jupyter_output><empty_output><jupyter_text>I currently have one empty community.<jupyter_code>clusters.loc[1017].to_frame().T<jupyter_output><empty_output><jupyter_text>But the reported number of cases there is zero.<jupyter_code>clusters_csse.loc[clusters_csse.index[-1], 1017]<jupyter_output><empty_output><jupyter_text>There is one community where the reported number of cases is larger than the population...<jupyter_code>idx = (clusters_csse.loc[clusters_csse.index[-1]] > clusters['population']).to_numpy()
np.sum(idx)
clusters.loc[clusters.index[idx]]
clusters_csse.loc[clusters_csse.index[-1], clusters.index[idx]].to_frame().T<jupyter_output><empty_output><jupyter_text>So I correct this:<jupyter_code>for t in clusters_csse.index:
    x = np.min(np.array([clusters_csse.loc[t].to_numpy(), clusters['population'].to_numpy()], dtype=np.int_), axis=0)
    clusters_csse.loc[t] = x

clusters_csse
idx = (clusters_csse.loc[clusters_csse.index[-1]] > clusters['population']).to_numpy()
np.sum(idx)
total_csse = clusters_csse.sum(axis=1).to_frame()
total_csse.rename(columns={0: 'omega'}, inplace=True)
total_csse<jupyter_output><empty_output><jupyter_text>Write<jupyter_code>key = '/clustering/cssegi'
with pd.HDFStore(resfile, complevel=complevel, complib=complib) as store:
    store[str(key)] = clusters_csse<jupyter_output><empty_output><jupyter_text>Plots<jupyter_code>figdir = Path('../figures') / '3-import_cssegi'
if not figdir.is_dir():
    figdir.mkdir(parents=True, exist_ok=True)<jupyter_output><empty_output><jupyter_text>Plot the total evolution<jupyter_code># parameters
figsize = (6,4.5)
dpi = 300
ms=2
lw=1


fig = plt.figure(facecolor='w', figsize=figsize)
ax = fig.gca()

ax.plot(total_csse.index.to_numpy(), total_csse.sum(axis=1).to_numpy(), '-', ms=ms, color='darkblue')

ax.set_yscale('log')
ax.set_xlim(total_csse.index[0],None)
ax.set_ylabel("$T$", fontsize="medium")
plt.xticks(rotation=45)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.tick_params(left=True, labelleft=True, bottom=True, labelbottom=True)
ax.tick_params(axis='both', length=4)
fig.tight_layout()

fname = 'omega_tot_cssegi'
for ext in exts:
    filepath = figdir / (fname + ext)
    fig.savefig(filepath, bbox_inches='tight', pad_inches=0, dpi=dpi)
    print("Written file: {:s}".format(str(filepath)))
fig.clf()
plt.close('all')
filepath = figdir / (fname + '.png')
Image(filename=filepath, width=4./3*360)<jupyter_output><empty_output><jupyter_text>Plot the spatial evolution as a profile<jupyter_code>clusters_csse.iloc[-1].max()/1.0e6
# parameters
dpi=150
fps=10
figsize=(6, 4.5)
lw=0.5
ms=4
ylabel="$T_a$"

fileout = figdir / 'dynamic_profile.mp4'
mydir = figdir / 'profiles'
if not mydir.is_dir():
    mydir.mkdir(parents=True, exist_ok=True)
plot_omega_profile(np.array([clusters_csse.to_numpy()]), clusters_csse.index.to_list(), colors=['red'], \
                   fileout=fileout, tpdir=mydir, dpi=dpi, fps=fps, figsize=figsize, ylabel=ylabel, \
                   lw=lw, ms=ms, styles=['o'], deletetp=False, exts=['.png','.svg'], ymin=1., ymax=1.3e6)
HTML("""
<video height="480" controls>
  <source src="{:s}" type="video/mp4">
</video>
""".format(str(fileout)))<jupyter_output><empty_output><jupyter_text>Plot the spatial evolution on a map<jupyter_code>clusters_csse.iloc[-1].max()/1.0e6
fileout = figdir / 'dynamic_profile.mp4'
mydir = figdir / 'maps'
if not mydir.is_dir():
    mydir.mkdir(parents=True, exist_ok=True)
fileout = figdir / 'dynamic_map.mp4'
# parameters
dpi=150
fps=10
figsize=(6, 4.5)
lw=0.5
ms=4
idump=1


fileout = figdir / 'map_T.mp4'
plot_omega_map(clusters_csse.to_numpy(), clusters_csse.index.to_list(), XY=clusters.loc[:, ['X', 'Y']].to_numpy().T, \
fileout=fileout, tpdir=mydir, dpi=dpi, fps=fps, figsize=figsize, idump=idump, \
               clabel="$T$", vmin=1., vmax=1.3e6, deletetp=False, exts=['.png','.svg'])
HTML("""
<video height="480" controls>
  <source src="{:s}" type="video/mp4">
</video>
""".format(str(fileout)))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15257202846144874, 'fraction_numerical': 0.012481045141724018, 'mean_word_length': 6.703504043126685, 'pattern_counts': {'":': 0, '<': 71, '<?xml version=': 0, '>': 73, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '506912', 'license': 'MIT', 'n_tokens_mistral': 3379, 'n_tokens_neox': 3281, 'n_words': 576, 'path': 'notebooks/3-import_cssegi.ipynb', 'repo_name': 'gletreut/epidemiology_flux_model', 'size': '8573'}