<jupyter_start><jupyter_text>Nearest neighbors[Index](0-index.ipynb)Here we compute matrices the matrix of distances between communities, $[d_{ab}]$.  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 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 show_image, geo_dist
resdir = Path('../results/')
if not resdir.is_dir():
    raise ValueError('No results directory!')
resfile = resdir / 'safegraph_analysis.hdf5'
complevel=7
complib='zlib'
with pd.HDFStore(resfile, complevel=complevel, complib=complib) as store:
    print(f"File {resfile.stem} has {len(store.keys())} entries.")
tfmt = '%Y-%m-%d'
figdir = Path('../figures') / '4-distances'
if not figdir.is_dir():
    figdir.mkdir(parents=True, exist_ok=True)
exts = ['.png', '.svg']<jupyter_output><empty_output><jupyter_text>FunctionsWe use the [formula](https://en.wikipedia.org/wiki/Geographical_distance) of the Spherical Earth projected to a plane. for the 'Spherical Earth projected to a plane':\begin{equation}d(M_1,M_2) = R \sqrt{\Delta \phi^2 + \cos{\bar{\phi}} \Delta \lambda^2},\end{equation}where $\phi$ is the latitude, $\lambda$ is the longitude, and $\bar{\phi} = (\phi_1 + \phi_2)/2$. $R = 6, 371.009 \, \mathrm{km}$ is the radius of the Earth.Here `X` is the longitude and `Y` is the latitude (see this [notebook](06-cross_compare_cssegi.ipynb)).<jupyter_code>R_earth = 6.371009e6<jupyter_output><empty_output><jupyter_text>Clusters  Load clusters<jupyter_code>key = "/clustering/clusters"
with pd.HDFStore(resfile, complevel=complevel, complib=complib) as store:
    clusters = store[key]
# clusters = pd.read_hdf(resfile, key)
N = len(clusters)
print(f"N = {N}")
clusters<jupyter_output>N = 1024<jupyter_text>Compute distance matrix using Spherical Earth projected to a plane formula<jupyter_code>XY = clusters.loc[:, ['X','Y']].astype('float64').to_numpy()

Long, Lat = XY.T * np.pi/180.

LLat1, LLat2 = np.meshgrid(Lat, Lat, indexing='ij')
dLLat = LLat1 - LLat2
LLat_bar = 0.5*(LLat1 + LLat2)

LLong1, LLong2 = np.meshgrid(Long, Long, indexing='ij')
dLLong = LLong1 - LLong2

D = R_earth * np.sqrt(dLLat**2 + np.cos(LLat_bar)*dLLong**2)
df_dist = pd.DataFrame(data=D, index=clusters.index, columns=clusters.index)

key = "/clustering/distances"
with pd.HDFStore(resfile, complevel=complevel, complib=complib) as store:
    store[key] = df_dist<jupyter_output><empty_output><jupyter_text>Compare with [equirectangular projection](https://en.wikipedia.org/wiki/Equirectangular_projection) distances<jupyter_code>P = clusters.loc[:, ['Y', 'X']].to_numpy()
PP = np.einsum('ik,jk->ij',P,P)
U = np.ones((N,1))
Delta = PP.diagonal().reshape(-1,1)
Dequi = (U*Delta.T) + (Delta*U.T) - 2*PP
Dequi = R_earth * np.pi/180.* np.sqrt(Dequi)
err = np.linalg.norm(D-Dequi)/np.linalg.norm(D)
print("Equirectangular projection error = {:.1f} %".format(err*100))<jupyter_output>Equirectangular projection error = 15.1 %<jupyter_text>Some statistics on the actual separation<jupyter_code>key = "/clustering/distances"
with pd.HDFStore(resfile, complevel=complevel, complib=complib) as store:
    df_dist = store[key]

df_dist
Ns = df_dist.index.to_numpy()
dists = df_dist.to_numpy()
dists_sorted = np.sort(dists, axis=-1)

# median distance of 1st, 2nd, 3rd...etc neighbor
dist_med = np.array([np.median(dists_sorted[:,n]) for n in Ns[1:]])
# parameters
figsize = (6,4.5)
dpi = 300
ms=2
lw=1

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

ax.plot(Ns[1:], dist_med/1.0e3, '-', ms=ms, color='darkblue')

# ax.set_xscale('log')
# ax.set_yscale('log')
ax.set_xlim(Ns[1],None)
ax.set_xlabel("kth nearest neighbor", fontsize="medium")
ax.set_ylabel("median distance (km)", fontsize="medium")
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 = 'median_distance_nearest_neighbors'
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>
{'fraction_non_alphanumeric': 0.15627668659265584, 'fraction_numerical': 0.019427839453458582, 'mean_word_length': 5.982116244411326, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 24, 'https://': 2, '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': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '413592', 'license': 'MIT', 'n_tokens_mistral': 1851, 'n_tokens_neox': 1777, 'n_words': 380, 'path': 'notebooks/4-distances.ipynb', 'repo_name': 'gletreut/epidemiology_flux_model', 'size': '4684'}