# Import modules
import numpy as np
import pandas as pd
import geopandas as gpd
import seaborn as sns
import networkx as nx
import osmnx as ox
import re
from shapely.geometry import Point, LineString, Polygon
import matplotlib.pyplot as plt
import os
import time
import warnings
import IPython
import requests
from IPython.display import display, clear_output
from shapely.ops import nearest_points   #for pharmacy_setting function
from scipy.stats import kruskal
from matplotlib import colors
from matplotlib.lines import Line2D
from matplotlib.animation import FuncAnimation
from matplotlib.colors import Normalize
from matplotlib.cm import ScalarMappable
import matplotlib.patheffects as path_effects
from tabulate import tabulate
from matplotlib.colors import ListedColormap
from mpl_toolkits.axes_grid1.anchored_artists import AnchoredSizeBar
from matplotlib.font_manager import FontProperties
from matplotlib.patches import Patch
import imageio

warnings.filterwarnings("ignore")
print('\n'.join(f'{m.__name__}=={m.__version__}' for m in globals().values() if getattr(m, '__version__', None)))

numpy==1.22.0
pandas==1.3.5
geopandas==0.10.2
seaborn==0.11.2
networkx==2.6.3
osmnx==1.1.2
re==2.2.1
IPython==8.3.0
requests==2.27.1
imageio==2.13.5

# Check working directory
os.getcwd()

# Use to set work directory properly
if os.path.basename(os.getcwd()) == 'code':
    os.chdir('../../')
os.getcwd()

# Read in pharmacy data
pharmacies = gpd.read_file('./data/raw/public/pharmacy/pharmacies_df.gpkg')
pharmacies.head()

pharmacies.explore()

# staffing_interpolated.csv is confidential. 
# with permission from researchers, download from OSF at https://osf.io/download/b956m/ 
#if not os.path.exists("./data/raw/private/staffing_interpolated.csv"):
#    print("Loading pharmacy staffing levels", flush=True)
#    url = 'https://osf.io/download/b956m/'
#    r = requests.get(url, allow_redirects=True)
#    open('./data/raw/private/staffing_interpolated.csv', 'wb').write(r.content)
    
# Load pharmacy staffing data 
pharm_staffing = pd.read_csv('./data/raw/private/staffing_interpolated.csv')
pharm_staffing.head()

# Join private staffing data to pharmacies dataset (locations, hours of operations, etc.) 
pharmacies_df = pd.merge(pharmacies, pharm_staffing, on='pharmid', how = 'left')
pharmacies_df

# Read in population data by town
population_df = gpd.read_file('./data/raw/public/population/tidycensus_population.gpkg')

# Read in metropolitan / micropolitan classifications (NECTAS)
nectas_df =  gpd.read_file('./data/raw/public/population/nectas.csv')

# Join NECTAS classifications to population data with subdivision and county FIPS
pop_df = pd.merge(population_df, nectas_df[['fips_subdivision', 'necta', 'fips_county', 'fips_state']], on=['fips_subdivision','fips_county', 'fips_state'], how = 'left')

# Save as geopackage into public/derived
pop_df.to_file('./data//derived/public/pop_data.gpkg', driver='GPKG')

pop_df.head()

# All geometries for which we have obtained population data
pop_df.plot()

<AxesSubplot:>

%%time
# To create a new graph from OpenStreetMap, delete or rename the graph file (if it exists)
# AND set OSM to True
# This is more likely to work on a local computer than CyberGISX 
# This network was created on Windows with Python 3.12.1 and osmnx 1.8.1 on 2024-01-24
# using the 00_road_network.ipynb script

OSM = False

# Define the place name for Vermont
place_name_vermont = 'Vermont, USA'

roads_path = "./data/raw/private/osm_roads.graphml"

# if buffered street network is not saved, and OSM is preferred, generate a new graph from OpenStreetMap and save it
if not os.path.exists(roads_path) and OSM:
    print("Loading buffered Vermont road network from OpenStreetMap. Please wait... runtime may exceed 9min...", flush=True)
    G = ox.graph_from_place(place_name_vermont, network_type='drive', buffer_dist=64373.8) 
    print("Saving road network to", roads_path, " Please wait...", flush=True)
    ox.save_graphml(G, roads_path)
    print("Data saved.")
    
# otherwise, if buffered street network is not saved, download graph from the OSF project
elif not os.path.exists(roads_path):
    print("Downloading buffered Vermont road network from OSF...", flush=True)
    url = 'https://osf.io/download/bwrmu/'  
    r = requests.get(url, allow_redirects=True)
    print("Saving road network to", roads_path, " Please wait...", flush=True)
    open(roads_path, 'wb').write(r.content)
    
# load the saved network graph
if os.path.exists(roads_path):
    print("Loading road network from", roads_path, "Please wait...", flush=True)
    G = ox.load_graphml(roads_path) 
    print("Data loaded.") 
else:
    print("Error: could not load the road network from file.")

Loading road network from ./data/raw/private/osm_roads.graphml Please wait...
Data loaded.
CPU times: user 1min 21s, sys: 5.51 s, total: 1min 27s
Wall time: 1min 41s

%%time
ox.plot_graph(G, node_size = 1, bgcolor = 'white', node_color = 'black', edge_color = "#333333", node_alpha = 0.2, edge_linewidth = 0.5)

CPU times: user 1min 8s, sys: 2.67 s, total: 1min 11s
Wall time: 1min 17s

(<Figure size 576x576 with 1 Axes>, <AxesSubplot:>)

%%time
# Turn network edges into a geodataframe
edges = ox.graph_to_gdfs(G, nodes=False, edges=True)

# Count frequency of each speed value
speed_values = edges['maxspeed'].value_counts()

# Ouput number of edges and frequences of speed values
print(str(len(edges)) + " edges in graph")
print(speed_values.to_string())

# view all highway types
print(edges['highway'].value_counts())

%%time
ox.speed.add_edge_speeds(G)

%%time
ox.speed.add_edge_travel_times(G)

%%time
# Create point geometries for each node in the graph, to make constructing catchment area polygons easier
for node, data in G.nodes(data=True):
    data['geometry']=Point(data['x'], data['y'])

%%time
# Turn network edges into a geodataframe
nodes,edges = ox.graph_to_gdfs(G, nodes=True, edges=True)

# Count frequency of each speed value
speed_values = edges['speed_kph'].value_counts()

# Ouput number of edges and frequences of speed values
print(str(len(edges)) + " edges in graph")
print(speed_values.to_string())

# Drop four dead-end nodes for parking lot entrances
nodes.drop([205007938, 5976921845, 204567135, 61755372], inplace = True)

# calculate a color scheme for edges based on speed
ec = ox.plot.get_edge_colors_by_attr(G, attr="speed_kph", cmap='YlOrRd')

# plot edge speeds
fig, ax = ox.plot_graph(G, node_size=0, edge_color=ec, bgcolor="white")

nodes_osmid = gpd.GeoDataFrame(nodes[["geometry"]]).reset_index()
# print(nodes_osmid.head())

pharmacies_osm = gpd.sjoin_nearest(pharmacies_df, nodes_osmid, distance_col="distances")

#rename column from osmid to nearest_osm, so that it works with other code
pharmacies_osm = pharmacies_osm.rename(columns={"osmid": "nearest_osm"})

distances = [600, 1200, 1800] # Distances in travel time (seconds!) (10 min, 20 min, 30 min)
weights = [1.0, 0.68, 0.22] # Relative weights of 10 minute catchment, 20 minute catchment, and 30 minute catchment

# initialize catchment list, composed of 3 empty geodataframes for 3 distance bands
catchments = []
for distance in distances:
    catchments.append(gpd.GeoDataFrame())
    
# initialize results geodataframe to store final catchment areas
results = gpd.GeoDataFrame(columns = ["geometry","pharmid","weight"], crs = "EPSG:4326", geometry = "geometry")

# calculate catchment areas
print("Working on pharmacy:")
for ind in pharmacies_osm.index:
    print(str(ind + 1) + "/" + str(len(pharmacies_osm)), end = '\r')
    # print(pharmacies_osm['pharmid'][ind], end=' ')
    ## CREATE DICTIONARIES ##
    # create dictionary of nearest nodes
    nearest_nodes_30 = nx.single_source_dijkstra_path_length(G, pharmacies_osm['nearest_osm'][ind], distances[2], "travel_time") # creating the largest graph from which 10 and 20 minute drive times can be extracted from

    # extract values within 20 and 10 (respectively) minutes drive times
    nearest_nodes_20 = dict()
    nearest_nodes_10 = dict()
    for key, value in nearest_nodes_30.items():
        if value <= distances[1]:
            nearest_nodes_20[key] = value
        if value <= distances[0]:
            nearest_nodes_10[key] = value

    ## CREATE POLYGONS FOR 3 DISTANCE CATEGORIES (10 min, 20 min, 30 min) ##

    # 30 MIN
    # If the graph already has a geometry attribute with point data,
    # this line will create a GeoPandas GeoDataFrame from the nearest_nodes_30 dictionary
    points_30 = gpd.GeoDataFrame(gpd.GeoSeries(nx.get_node_attributes(G.subgraph(nearest_nodes_30), 'geometry')))

    # This line converts the nearest_nodes_30 dictionary into a Pandas data frame and joins it to points
    # left_index=True and right_index=True are options for merge() to join on the index values
    points_30 = points_30.merge(pd.Series(nearest_nodes_30).to_frame(), left_index=True, right_index=True)

    # Re-name the columns and set the geodataframe geometry to the geometry column
    points_30 = points_30.rename(columns={'0_x':'geometry','0_y':'z'}).set_geometry('geometry')

    # Create a convex hull polygon from the points
    polygon_30 = gpd.GeoDataFrame(gpd.GeoSeries(points_30.unary_union.convex_hull))
    polygon_30 = polygon_30.rename(columns={0:'geometry'}).set_geometry('geometry')
    polygon_30["weight"] = weights[2]

    # 20 MIN # 1200 seconds!
    # Select nodes less than or equal to 20
    points_20 = points_30.query("z <= " + str(distances[1]))

    # Create a convex hull polygon from the points
    polygon_20 = gpd.GeoDataFrame(gpd.GeoSeries(points_20.unary_union.convex_hull))
    polygon_20 = polygon_20.rename(columns={0:'geometry'}).set_geometry('geometry')
    polygon_20["weight"] = weights[1]

    # 10 MIN # 600 seconds!
    # Select nodes less than or equal to 10
    points_10 = points_30.query("z <= " + str(distances[0]))

    # Create a convex hull polygon from the points
    polygon_10 = gpd.GeoDataFrame(gpd.GeoSeries(points_10.unary_union.convex_hull))
    polygon_10 = polygon_10.rename(columns={0:'geometry'}).set_geometry('geometry')
    polygon_10["weight"] = weights[0]
    
    # Clip the overlapping distance ploygons (create two donuts + hole)
    polygon_30_hole = gpd.overlay(polygon_30, polygon_20, how ="difference")
    polygon_20_hole = gpd.overlay(polygon_20, polygon_10, how ="difference")
    
    # Create dataframe combining polygon_10, polygon_20, polygon_30
    polygons = pd.concat([polygon_10, polygon_20_hole, polygon_30_hole])
    polygons.set_crs(crs="EPSG:4326", inplace = True)
    # polygons = pd.concat([polygon_10, polygon_20, polygon_30])
    # polygons.set_crs(crs="EPSG:4326", inplace = True)
    
    polygons_gdf = polygons
    polygons_gdf["pharmid"] = pharmacies_osm['pharmid'][ind]
    #polygons_gdf["weight"] = weights

    results = pd.concat([results, polygons_gdf])

print(len(results) / 192) # Value of 3.0 confirms all pharmacies, since there are 192 pharmacies in our study area and 3 polygons per pharmacy

3.0

results.explore("weight", categorical=True, cmap='GnBu', style_kwds={'fillOpacity': 0.3})
# try results.explore("weight")

# Change CRS to match prior to analysis, should be mapped in VT State Plane projection
results.to_crs("EPSG:6589", inplace = True)
pop_df.to_crs("EPSG:6589", inplace = True)

# Calculate town areas
pop_df['town_area'] = pop_df.geometry.area
#display(pop_df)

# Run the overlay to find intersection of fragments (fragments are where each individual catchment area overlaps with a different town geometry)
fragments = gpd.overlay(pop_df, results, how = 'intersection')

# Calculate fragment areas
fragments['frag_area'] = fragments.geometry.area

# Calculate area ratios to see how much of a town is covered by the given fragment --> used to estimate the population in that fragment
fragments['area_ratio']= fragments['frag_area'] / fragments['town_area']

# Calculate fragment value by multiplying area_ratio by distance weight
fragments['value'] = fragments['weight'] * fragments['area_ratio']

# Calculate population served in each fragment
fragments['pop_value'] = fragments['total_pop']*fragments['value']

fragments.head()

# Sum population served per pharmacy 
sum_values = fragments[["pharmid","pop_value"]]
sum_values = sum_values.groupby(by = ['pharmid']).sum('pop_value')
sum_values.rename(columns={'pop_value': 'total_value'}, inplace = True) 
sum_values.head()

# Join summed values by pharmacy to pharmacies_df
pop_join = sum_values.merge(pharmacies_df, on='pharmid', how = 'left')
pd.set_option('display.max_columns', None)

# Only display this data for testing. It reveals confidential information at the pharmacy level.
# pop_join.head()

# Calculate number of staff that work at each pharmacy
# Weekday
pop_join['week_staff'] = pop_join['week_pharm'] + 0.5 * pop_join['week_tech']

#Saturday
pop_join['sat_staff'] = pop_join['sat_pharm'] + 0.5 * pop_join['sat_tech']

# Sunday
pop_join['sun_staff'] = pop_join['sun_pharm'] + 0.5 * pop_join['sun_tech']

# Calculate service to pop ratio for each pharmacy for Weekdays, Saturdays, and Sundays between 7am-11pm
def calculate_sp_ratios(pop_join, day_type):
    pop_join[f'sp_ratio_{day_type}'] = pop_join[f'{day_type}_staff'] / pop_join['total_value']
    
    for hour in range(7, 24):
        pop_join[f'sp_ratio_{hour}{day_type}'] = np.where((pop_join[f'{day_type}_open'] <= hour) & (pop_join[f'{day_type}_close'] >= hour), 
                                                       pop_join[f'{day_type}_staff'] / pop_join['total_value'], 
                                                       np.nan)
# Weekday
calculate_sp_ratios(pop_join, 'week')

# Saturday
calculate_sp_ratios(pop_join, 'sat')

# Sunday
calculate_sp_ratios(pop_join, 'sun')

# Only display this data for testing. It reveals confidential information at the pharmacy level.
# pop_join.head()

# Join to the fragments df
frag_join = fragments.merge(pop_join, on = 'pharmid', how = 'left')

# Only display this data for testing. It reveals confidential information at the pharmacy level.
# frag_join.head()

# Calculate weighted service to population ratio for each slice of time (weight correspond to distance from the pharmacy still)
frag_join['sp_weighted_w'] =  frag_join['sp_ratio_week'] * frag_join['value'] 
frag_join['sp_weighted_s'] =  frag_join['sp_ratio_sat'] * frag_join['value']  
frag_join['sp_weighted_su'] =  frag_join['sp_ratio_sun'] * frag_join['value'] 

def calculate_weighted_sp_ratios(frag_join, day_type):
    for hour in range(7, 24):
        frag_join[f'sp_weighted_{hour}{day_type}'] = frag_join[f'sp_ratio_{hour}{day_type}'] * frag_join['value']

# Weekday
calculate_weighted_sp_ratios(frag_join, 'week')

# Saturday
calculate_weighted_sp_ratios(frag_join, 'sat')

# Sunday
calculate_weighted_sp_ratios(frag_join, 'sun')

# Make all NaN values 0, since these values represent 0 access
frag_join.fillna(0, inplace = True)

# Sum weighted service to population ratio by town
accessibility_towns = frag_join.groupby('GEOID').sum()

accessibility_towns.head()

#Normalize (multiply by large constant)
#Full Days
accessibility_towns['access_w'] = accessibility_towns['sp_weighted_w']*10000
accessibility_towns['access_s'] = accessibility_towns['sp_weighted_s']*10000
accessibility_towns['access_su'] = accessibility_towns['sp_weighted_su']*10000

def calculate_accessibility_times(accessibility_towns, day_type):
    for hour in range(7, 24):
        accessibility_towns[f'access_{hour}{day_type}'] = accessibility_towns[f'sp_weighted_{hour}{day_type}'] * 10000

# Weekday times
calculate_accessibility_times(accessibility_towns, 'week')

# Saturday times
calculate_accessibility_times(accessibility_towns, 'sat')

# Sunday times
calculate_accessibility_times(accessibility_towns, 'sun')

#Select just access columns 
accessibility_towns = accessibility_towns[[col for col in accessibility_towns.columns if col.startswith('access')]]
accessibility_towns.head()

# Join back to geometries in pop_df for mapping. Pop_df contains town geometries
mapping_df = pd.merge(accessibility_towns, pop_df[['GEOID', 'geometry', 'fips_state', 'necta', 'pct_elderly']], on='GEOID', how='left')
mapping_df = mapping_df[mapping_df['fips_state'].isin(['50'])]
mapping_df = gpd.GeoDataFrame(mapping_df, crs="EPSG:6589", geometry='geometry')
mapping_df.explore("access_w", tooltip=["access_w", "access_s", "access_su"], cmap="Greens", scheme='Quantiles', k=6, legend_kwds={'scale':False}) # View geometries of interest with weekday, sat, sun summaries

# Prepare pharmacies for mapping
pharmacies_df.to_crs("EPSG:6589", inplace = True) # Match CRS for Vermont
pharmacies_map = pharmacies_df[pharmacies_df['state'] == 'VT']

# Replace NaN values in 'necta' column with 'Rural'. All towns that are not metropolitan or micropolitan are rural. 
mapping_df['necta'].fillna('Rural', inplace=True)

# Replace 'Metropolitan NECTA' and 'Micropolitan NECTA' with 'Metropolitan' and 'Micropolitan', respectively
mapping_df['necta'] = mapping_df['necta'].replace({'Metropolitan NECTA': 'Metropolitan', 'Micropolitan NECTA': 'Micropolitan'})

necta_summary = mapping_df[['necta', 'GEOID']].groupby("necta").count()
necta_summary.rename(columns={"GEOID": "N"}, inplace=True)
print(necta_summary)

                N
necta            
Metropolitan   35
Micropolitan   48
Rural         172

# Calculate population denisty by town
popdensity_df = pop_df.copy()
popdensity_df = popdensity_df[popdensity_df['fips_state'].isin(['50'])]
popdensity_df.to_crs("EPSG:6589", inplace = True)
popdensity_df['town_area_km2'] = popdensity_df['town_area'] / 10**6
popdensity_df['pop_density'] = popdensity_df['total_pop']/popdensity_df['town_area_km2']
popdensity_df['pop_density'] = popdensity_df['pop_density'].replace([np.inf, -np.inf], 0) #make inf values zero since they arise from total_pop = 0

# Custom color map for pop density
base_cmap = plt.get_cmap('Blues')
num_colors = 256

custom_colors = []
for i in range(num_colors):
    color = colors.rgb_to_hsv(base_cmap(i)[:3])  # Convert RGB to HSV
    color[1] *= 0.7  # Reduce saturation by one third
    custom_colors.append(colors.hsv_to_rgb(color))  # Convert back to RGB

custom_cmap = colors.LinearSegmentedColormap.from_list("custom_blues", custom_colors, N=num_colors)

# Map Population Density
fig1, ax = plt.subplots(figsize=(12.3, 10.8))
popdensity_df.plot(column='pop_density', scheme = 'jenkscaspallforced',k=5, cmap=custom_cmap, legend = False, ax=ax)

metropolitan_df = mapping_df[mapping_df['necta'] == 'Metropolitan']
micropolitan_df = mapping_df[mapping_df['necta'] == 'Micropolitan']

metropolitan_boundary = metropolitan_df.dissolve(by='necta')['geometry'].boundary
micropolitan_boundary = micropolitan_df.dissolve(by='necta')['geometry'].boundary

# Plot only the merged exterior boundaries of the 'Metropolitan NECTA' group
metropolitan_boundary.plot(ax=ax, color='black', linewidth=1.5)

# Plot only the merged exterior boundaries of the 'Micropolitan NECTA' group
micropolitan_boundary.plot(ax=ax, color='black', linestyle = 'dashed', linewidth=1.5)

vt_boundary = popdensity_df.dissolve().boundary
vt_boundary.plot(ax=ax, color='black', linewidth=.75)

# Legend
legend_elements = [
    Line2D([0], [0], color='black', lw=2, label='Metropolitan'),
    Line2D([0], [0], color='black', lw=2, linestyle= 'dashed', label='Micropolitan'),
    Line2D([0], [0], marker='o', color='w', markerfacecolor='white', markeredgecolor='black', markersize=9, label='Pharmacy')
]

pd_colorbar = plt.cm.ScalarMappable(cmap=custom_cmap, norm=plt.Normalize(vmin=0, vmax=2046.6))
cbar = plt.colorbar(pd_colorbar, shrink = .38, pad=.001, label = 'Population Density (person/km^2)', location="bottom")
ax.legend(handles=legend_elements, loc='lower right', fontsize=10, bbox_to_anchor=(.95,.2))

# Plot pharmacy locations
pharmacies_map.plot(ax=ax, marker='o', facecolor='white', edgecolor='black', markersize=50, label='VT Pharmacies', legend = True)
plt.text(0.22, 0.71, 'Burlington', transform=ax.transAxes, size=11, weight="bold", color='white',
         path_effects=[path_effects.Stroke(linewidth=2, foreground='black'), path_effects.Normal()])

#fig.patch.set_edgecolor('black') #Border
#fig.patch.set_linewidth(2)

plt.axis('off')
plt.show()

#Export figure to results/figures
fig1.savefig("./results/figures/figure1.jpg", dpi=300)

# Weekday accessibility by classification table
means_by_metro = mapping_df.groupby('necta').mean()
weekdaymean_by_metro = means_by_metro[['access_w']]

weekdaymean_by_metro.columns = ['Mean Access']
table_1 = tabulate(weekdaymean_by_metro, headers='keys', tablefmt='simple_grid')
print(table_1)

necta           Mean Access
------------  -------------
Metropolitan        4.19849
Micropolitan        4.2507
Rural               3.01231

# set global maximum for mapping accessibility scores
maxacc = 8

# Map of weekday pharmacy accessibility
fig2, ax = plt.subplots(figsize=(12, 12), facecolor = 'white')
mapping_df.plot(column='access_w', cmap='Greens', legend = False, ax=ax, vmax=maxacc)
mapping_df.dissolve().boundary.plot(ax=ax, color='black', linewidth=1)

# Dissolve the geometries for each group to merge adjacent polygons 
metropolitan_df = mapping_df[mapping_df['necta'] == 'Metropolitan']
micropolitan_df = mapping_df[mapping_df['necta'] == 'Micropolitan']

metropolitan_boundary = metropolitan_df.dissolve(by='necta')['geometry'].boundary
micropolitan_boundary = micropolitan_df.dissolve(by='necta')['geometry'].boundary

# Plot exterior boundaries of the 'Metropolitan NECTA' group
metropolitan_boundary.plot(ax=ax, color='black', linewidth=1.2)

# Plot exterior boundaries of the 'Micropolitan NECTA' group
micropolitan_boundary.plot(ax=ax, color='black',linestyle = "dashed", linewidth=1.5)

legend_elements = [
    Line2D([0], [0], color='black', lw=2, label='Metropolitan'),
    Line2D([0], [0], color='black', linestyle = 'dashed',lw=2, label='Micropolitan')
]
ax.legend(handles=legend_elements, loc='lower right', fontsize=12, bbox_to_anchor=(1,.2))
colorbar = plt.cm.ScalarMappable(cmap='Greens', norm=plt.Normalize(vmin=0, vmax=maxacc))

#fig.patch.set_edgecolor('black') #Border
#fig.patch.set_linewidth(2)  
#plt.subplots_adjust(right = .9, top = .9)

# Colorbar settings
cbar = plt.colorbar(colorbar, shrink = .5, pad=.02, label = 'Accessibility', location="bottom")
cbar.set_ticks([0, 2, 4, 6, 8])
    
plt.axis('off')
plt.show()

# Save figure
fig2.savefig("./results/figures/figure2.jpg", dpi=300)

# Check for normal distribution in weekday metro, micro, rural
access_w_metro = mapping_df[mapping_df['necta'] == 'Metropolitan']['access_w']
access_w_micro = mapping_df[mapping_df['necta'] == 'Micropolitan']['access_w']
access_w_rural = mapping_df[mapping_df['necta'] == 'Rural']['access_w']

plt.hist(access_w_rural, bins=10)
plt.hist(access_w_metro, bins=10)
plt.hist(access_w_micro, bins=10)
print("Not normal distribution. Cannot use ANOVA test. Use Kruskal-Wallis instead.")

Not normal distribution. Cannot use ANOVA test. Use Kruskal-Wallis instead.

# Kruskal Wallis Test for significant difference of means between necta classification during conventional business hours
h_statistic_1, p_value_1 = kruskal(access_w_metro, access_w_micro, access_w_rural)
                 
print("Kruskal-Wallis H Statistic:", h_statistic_1)
print("P-value:", p_value_1)

alpha = 0.05
if p_value_1 < alpha:
    print("Reject the null hypothesis. There is a significant difference in mean access during conventional weekday \nbusiness hours between metropolitan, micropolitan, and rural towns.")
else:
    print("Fail to reject the null hypothesis. There is no significant difference in mean access between groups.")

Kruskal-Wallis H Statistic: 27.690600217921656
P-value: 9.706497843353757e-07
Reject the null hypothesis. There is a significant difference in mean access during conventional weekday 
business hours between metropolitan, micropolitan, and rural towns.

# Mean accessibility by day table
mean_access_day = mapping_df[['access_w', 'access_s', 'access_su']].mean()
mean_access_day_df = mean_access_day.to_frame().rename(columns={0: 'Mean Access'})
mean_access_day_df.index = ['Weekday', 'Saturday', 'Sunday']
table_2 = tabulate(mean_access_day_df, headers='keys', tablefmt='simple_grid')
print(table_2)

            Mean Access
--------  -------------
Weekday         3.40823
Saturday        2.02426
Sunday          1.43059

# Map accessibility by day of the week
mapping_df1 = mapping_df
#print(mapping_df1[['access_w', 'access_s', 'access_su']].min().min())  # Min accessibility value
#print(mapping_df1[['access_w', 'access_s', 'access_su']].max().max())  # Max accessibility value

fig3, axs = plt.subplots(1, 3, figsize=(22.5, 10), facecolor = 'white')
plt.subplots_adjust(wspace=-.4)

for i, column in enumerate(['access_w', 'access_s', 'access_su']):
    ax = axs[i]
    mapping_df1.plot(column=column, cmap='Greens', linewidth=0.2, ax=ax, edgecolor='0.8', legend=False,
                vmin=0, vmax=maxacc)
    mapping_df1.dissolve().boundary.plot(ax=ax, color='black', linewidth=1)
    
    # Plot only the merged exterior boundaries of the 'Metropolitan NECTA' group
    metropolitan_boundary.plot(ax=ax, color='black', linewidth=.9)

    # Plot only the merged exterior boundaries of the 'Micropolitan NECTA' group
    micropolitan_boundary.plot(ax=ax, color='black', linestyle = 'dashed', linewidth=.9)
    
    axs[0].set_title(label='Weekday', fontsize=14)
    axs[1].set_title('Saturday', fontsize=14)
    axs[2].set_title('Sunday', fontsize=14)
    ax.axis('off')

cbar = fig.colorbar(plt.cm.ScalarMappable(norm=plt.Normalize(vmin=0, vmax=maxacc), cmap='Greens'), ax=axs, #Max was set to 20 for visualization purposes
                    orientation='horizontal', pad=.07, shrink=.5)
cbar.set_label('Accessibility')

#fig.patch.set_edgecolor('black') # Figure Border
#fig.patch.set_linewidth(2)       # Figure Border
plt.subplots_adjust(right=1)

cbar.ax.set_position([0.26, 0.05, .5, 0.05])

ax.legend(handles=legend_elements, loc='lower right', fontsize=12, bbox_to_anchor=(.9,-.129))

plt.show()

#Save Figure
fig3.savefig("./results/figures/figure3.jpg", dpi=300)

linegraph_df = mapping_df.copy()  

# Weekday
access_columns_w = ['access_7week','access_8week','access_9week','access_10week','access_11week','access_12week', 'access_13week','access_14week','access_15week','access_16week','access_17week','access_18week','access_19week', 'access_20week', 'access_21week', 'access_22week']
values_w = linegraph_df[access_columns_w].mean()  # Calculate the mean value for each access column

# Saturday
access_columns_s = ['access_7sat','access_8sat','access_9sat','access_10sat','access_11sat','access_12sat','access_13sat','access_14sat','access_15sat','access_16sat','access_17sat','access_18sat','access_19sat', 'access_20sat', 'access_21sat', 'access_22sat']
values_s = linegraph_df[access_columns_s].mean()  # Calculate the mean value for each access column

# Sunday
access_columns_su = ['access_7sun','access_8sun','access_9sun','access_10sun', 'access_11sun','access_12sun','access_13sun','access_14sun', 'access_15sun', 'access_16sun', 'access_17sun', 'access_18sun', 'access_19sun', 'access_20sun', 'access_21sun', 'access_22sun']
values_su = linegraph_df[access_columns_su].mean()  # Calculate the mean value for each access column

# X-axis values for all three days
common_x_values = ['7am', '8am', '9am', '10am', '11am','12pm','1pm', '2pm', '3pm', '4pm', '5pm', '6pm', '7pm', '8pm', '9pm', '10pm']

# Plot together
fig, ax = plt.subplots(figsize=(10, 6))
plt.plot(common_x_values, values_w, label='Weekday', color='black', linestyle='-')
plt.plot(common_x_values, values_s, label='Saturday', color='black', linestyle='-.')
plt.plot(common_x_values, values_su, label='Sunday', color='black', linestyle=':')
#plt.title('Access by Time of Day')
#plt.xlabel('Time of Day')
plt.ylabel('Mean Accessibility')
plt.xticks() 
plt.yticks()
plt.legend()
plt.show()

# save figure
fig.savefig("./results/figures/figure4.jpg", dpi=300)

fig, ax = plt.subplots(figsize=(10, 10), facecolor = "white")
mapping_df.plot(column='access_23week', cmap='Greens', ax=ax, linewidth=0.2, edgecolor= '.8', vmax=.5)
mapping_df.dissolve().boundary.plot(ax=ax, color='black', linewidth=1)

cb = plt.cm.ScalarMappable(cmap='Greens', norm=plt.Normalize(vmin=0, vmax=.5))
#cbar = plt.colorbar(cb, shrink = .5, orientation = 'horizontal', pad=0, label = 'Accessibility')

metropolitan_boundary.plot(ax=ax, color='black', linewidth=.8)
micropolitan_boundary.plot(ax=ax, color='black',linestyle = "dashed", linewidth=.8)

#fig.patch.set_edgecolor('black') #Border
#fig.patch.set_linewidth(2)  
plt.subplots_adjust(right = .5, top = .9)

# Colorbar and Legend
cbar = plt.colorbar(cb, shrink = .5, orientation = 'horizontal', pad=0, label = 'Accessibility')
ax.legend(handles=legend_elements, loc='lower right', fontsize=12, bbox_to_anchor=(.92,.135))

plt.axis('off')
plt.tight_layout()
plt.show()

# Save Figure
fig.savefig("./results/figures/figure5.jpg", dpi=300)

access_w_test = pd.Series(mapping_df['access_w'])
access_s_test = pd.Series(mapping_df['access_s'])
access_su_test = pd.Series(mapping_df['access_su'])

# Check for distribution between days of week. 
plt.hist(access_w_test, bins=10) 
plt.hist(access_s_test , bins=10) 
plt.hist(access_su_test, bins=10) 
print("Distribution is not normal. Cannot use ANOVA test. Use Kruskal-Wallis instead.")

Distribution is not normal. Cannot use ANOVA test. Use Kruskal-Wallis instead.

# Run Kruskal-Wallis Test
h_statistic_2, p_value_2 = kruskal(access_w_test, access_s_test, access_su_test)
                 
print("Kruskal-Wallis H Statistic:", h_statistic_2)
print("P-value:", p_value_2)

alpha = 0.05
if p_value_2 < alpha:
    print("Reject the null hypothesis. There is a significant difference in mean access between weekdays, Saturdays, and Sundays.")
else:
    print("Fail to reject the null hypothesis. There is no significant difference in mean access between days.")

Kruskal-Wallis H Statistic: 202.67482244348162
P-value: 9.76610740072804e-45
Reject the null hypothesis. There is a significant difference in mean access between weekdays, Saturdays, and Sundays.

# Accessibility by Day and metropolitan/micropolitan Table
means_by_metro = mapping_df.groupby('necta').mean()
means_by_metro = necta_summary.join(means_by_metro[['access_w','access_s','access_su']])

n_rural = means_by_metro.iat[0,0] + means_by_metro.iat[1,0]

means_by_metro.loc['Non-rural'] = [n_rural, 
                                   (means_by_metro.iat[0,0] * means_by_metro.iat[0,1] + means_by_metro.iat[1,0] * means_by_metro.iat[1,1]) / n_rural,
                                   (means_by_metro.iat[0,0] * means_by_metro.iat[0,2] + means_by_metro.iat[1,0] * means_by_metro.iat[1,2]) / n_rural,
                                   (means_by_metro.iat[0,0] * means_by_metro.iat[0,3] + means_by_metro.iat[1,0] * means_by_metro.iat[1,3]) / n_rural]

means_by_metro = means_by_metro.iloc[[3,2]]

means_by_metro.loc['Percent Change'] = ['-',
                                            (means_by_metro.iat[1,1] - means_by_metro.iat[0,1]) / means_by_metro.iat[0,1] * 100,
                                            (means_by_metro.iat[1,2] - means_by_metro.iat[0,2]) / means_by_metro.iat[0,2] * 100,
                                            (means_by_metro.iat[1,3] - means_by_metro.iat[0,3]) / means_by_metro.iat[0,3] * 100]

print(tabulate(means_by_metro, tablefmt = 'fancy_grid', headers=["","N","Weekday Mean Access", "Saturday Mean Access", "Sunday Mean Access"]))
#means_by_metro

╒════════════════╤═══════╤═══════════════════════╤════════════════════════╤══════════════════════╕
│                │ N     │   Weekday Mean Access │   Saturday Mean Access │   Sunday Mean Access │
╞════════════════╪═══════╪═══════════════════════╪════════════════════════╪══════════════════════╡
│ Non-rural      │ 83.0  │               4.22868 │                2.56468 │              1.95207 │
├────────────────┼───────┼───────────────────────┼────────────────────────┼──────────────────────┤
│ Rural          │ 172.0 │               3.01231 │                1.76347 │              1.17895 │
├────────────────┼───────┼───────────────────────┼────────────────────────┼──────────────────────┤
│ Percent Change │ -     │             -28.7648  │              -31.2401  │            -39.6052  │
╘════════════════╧═══════╧═══════════════════════╧════════════════════════╧══════════════════════╛

access_columns = ['access_7week','access_8week','access_9week','access_10week','access_11week','access_12week',
                  'access_13week','access_14week','access_15week','access_16week','access_17week','access_18week',
                  'access_19week', 'access_20week', 'access_21week', 'access_22week']
linegraph_df2 = mapping_df.groupby('necta')[access_columns].mean()
linegraph_df2 = linegraph_df2.T

access_columns2 = ['access_7sat','access_8sat','access_9sat','access_10sat','access_11sat','access_12sat',
                   'access_13sat','access_14sat','access_15sat','access_16sat','access_17sat','access_18sat',
                   'access_19sat', 'access_20sat', 'access_21sat', 'access_22sat']
linegraph_df3 = mapping_df.groupby('necta')[access_columns2].mean()
linegraph_df3 = linegraph_df3.T

access_columns3 = ['access_7sun','access_8sun','access_9sun','access_10sun', 'access_11sun','access_12sun',
                   'access_13sun','access_14sun', 'access_15sun', 'access_16sun', 'access_17sun', 'access_18sun', 
                   'access_19sun', 'access_20sun', 'access_21sun', 'access_22sun']
linegraph_df4 = mapping_df.groupby('necta')[access_columns3].mean()
linegraph_df4 = linegraph_df4.T

fig, axs = plt.subplots(3, 1, figsize=(10, 12), sharey = True)  

custom_ticks = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]
time_labels = ['7am','8am', '9am','10am','11am','12pm','1pm','2pm', '3pm','4pm','5pm','6pm', '7pm','8pm','9pm','10pm']


# Weekday 
linegraph_df2.plot(time_labels, kind='line', color= ['#8AB5DB','#6787A3','#3B4E5E'] , ax=axs[0], legend = None)  
axs[0].set_title('Weekday')
axs[0].set_xlabel('')
axs[0].set_ylabel('Mean Access Value')
axs[0].set_xticks(custom_ticks)  
axs[0].set_xticklabels(time_labels, rotation= 45) 
axs[0].grid(axis='y', linestyle='-', linewidth=0.08, color='black')

# Saturday
linegraph_df3.plot(time_labels, kind='line', color = ['#8AB5DB','#6787A3','#3B4E5E'] , ax=axs[1], legend = None) 
axs[1].set_title('Saturday')
axs[1].set_xlabel('')
axs[1].set_ylabel('Mean Access Value')
axs[1].set_xticks(custom_ticks) 
axs[1].set_xticklabels(time_labels, rotation = 45) 
axs[1].grid(axis='y', linestyle='-', linewidth=0.08, color='black')

# Sunday
linegraph_df4.plot(time_labels, kind='line', color = ['#8AB5DB','#6787A3','#2A3742'], ax=axs[2], legend= None)  
axs[2].set_title('Sunday')
axs[2].set_xlabel('')
axs[2].set_ylabel('Mean Access Value')
axs[2].set_xticks(custom_ticks)  
axs[2].set_xticklabels(time_labels, rotation = 45) 
axs[2].grid(axis='y', linestyle='-', linewidth=0.08, color='black')

handles, labels = axs[0].get_legend_handles_labels()
fig.legend(handles, labels, loc='lower center', ncol=3, fontsize=12, title='', title_fontsize='large', labels = ['Metropolitan','Micropolitan','Rural'], bbox_to_anchor=(0.5, -0.04))

#fig.patch.set_edgecolor('black') #Border
#fig.patch.set_linewidth(2)  

#plt.suptitle('Mean Access by Time of Day', x=0.51, y=1.0, fontsize=16, ha='center')
plt.tight_layout()
plt.show()

# Save Figure
fig.savefig("./results/figures/figure7.jpg", dpi=300)

# WEEKDAY
mapping_df['necta'].fillna('Rural', inplace=True)

access_columns = ['access_7week','access_8week','access_9week','access_10week','access_11week',
                  'access_12week', 'access_13week','access_14week','access_15week','access_16week',
                  'access_17week','access_18week','access_19week', 'access_20week', 'access_21week', 
                  'access_22week']

hourly_diff_df = mapping_df.groupby('necta')[access_columns].mean() # Mean access by hour of the day after grouping by NECTA
hourly_diff_df = hourly_diff_df.T
hourly_diff_df['Non-rural'] = (hourly_diff_df['Metropolitan'] * 35 + hourly_diff_df['Micropolitan'] * 48) / (35 + 48) #weighted average to 
hourly_diff_df['%_diff'] = ((hourly_diff_df['Rural'] - hourly_diff_df['Non-rural']) / hourly_diff_df['Non-rural']) * 100 #Percent diff between rural and non-rural
hourly_diff_df['%_diff']['access_22week']=0 #Make this value 0. Both accessibility values too small to take meaningful % difference. Both values virtually 0 access
hourly_diff_df['%_diff'] = hourly_diff_df['%_diff'].abs() #Take absolute value for plotting
hourly_diff_df

# SATURDAY
access_columns2 = ['access_7sat','access_8sat','access_9sat','access_10sat','access_11sat','access_12sat',
                   'access_13sat','access_14sat','access_15sat','access_16sat','access_17sat','access_18sat',
                   'access_19sat', 'access_20sat', 'access_21sat', 'access_22sat']

hourly_diff_df2 = mapping_df.groupby('necta')[access_columns2].mean() # Mean access by hour of the day after grouping by NECTA
hourly_diff_df2 = hourly_diff_df2.T
hourly_diff_df2['Non-rural'] = (hourly_diff_df2['Metropolitan'] * 35 + hourly_diff_df2['Micropolitan'] * 48) / (35 + 48) #weighted average to 
hourly_diff_df2['%_diff'] = ((hourly_diff_df2['Rural'] - hourly_diff_df2['Non-rural']) / hourly_diff_df2['Non-rural']) * 100 #Percent diff between rural and non-rural
hourly_diff_df2['%_diff']['access_22sat']=0 #Make this value 0. Both accessibility values too small to take meaningful % difference. Both values virtually 0 access
hourly_diff_df2['%_diff'] = hourly_diff_df2['%_diff'].abs() #Take absolute value for plotting
hourly_diff_df2

# SUN
access_columns3 = ['access_7sun','access_8sun','access_9sun','access_10sun', 'access_11sun','access_12sun',
                   'access_13sun','access_14sun', 'access_15sun', 'access_16sun', 'access_17sun', 'access_18sun', 
                   'access_19sun', 'access_20sun', 'access_21sun', 'access_22sun']

hourly_diff_df3 = mapping_df.groupby('necta')[access_columns3].mean() # Mean access by hour of the day after grouping by NECTA
hourly_diff_df3 = hourly_diff_df3.T
hourly_diff_df3['Non-rural'] = (hourly_diff_df3['Metropolitan'] * 35 + hourly_diff_df3['Micropolitan'] * 48) / (35 + 48) #weighted average to 
hourly_diff_df3['%_diff'] = ((hourly_diff_df3['Rural'] - hourly_diff_df3['Non-rural']) / hourly_diff_df3['Non-rural']) * 100 #Percent diff between rural and non-rural
hourly_diff_df3['%_diff']['access_22sun']=0 #Make this value 0. Both accessibility values too small to take meaningful % difference. Both values virtually 0 access
hourly_diff_df3['%_diff'] = hourly_diff_df3['%_diff'].abs() #Take absolute value for plotting
hourly_diff_df3

common_x_values = ['7am', '8am', '9am', '10am', '11am','12pm','1pm', '2pm', '3pm', '4pm', '5pm', '6pm', '7pm', '8pm', '9pm', '10pm']

# PLOT TOGETHER
figure = plt.figure(figsize=(10, 6))
plt.plot(common_x_values, hourly_diff_df['%_diff'], label='Weekday', color = 'black',linestyle='-')
plt.plot(common_x_values, hourly_diff_df2['%_diff'], label='Saturday', color = 'black', linestyle='-.')
plt.plot(common_x_values, hourly_diff_df3['%_diff'], label='Sunday', color = 'black', linestyle=':')

plt.ylabel('Urban-Rural Percent Change in Access')
plt.grid(axis='y', linestyle='-', linewidth=0.08, color='black')
plt.xticks(rotation=45)
plt.legend()
plt.tight_layout()
plt.show()

# Save Figure
figure.savefig("./results/figures/figure8.jpg", dpi=300)

def plot_accessibility_map(mapping_df, metropolitan_boundary, micropolitan_boundary, time_column, ax):
    mapping_df.plot(column=time_column, cmap='Greens', legend=False, ax=ax, vmax=maxacc)
    metropolitan_boundary.plot(ax=ax, color='black', linewidth=1.2)
    micropolitan_boundary.plot(ax=ax, color='black', linestyle = "dashed", linewidth=1.2)
    mapping_df.dissolve().boundary.plot(ax=ax, color='black', linewidth=1)
    ax.set_axis_off()
    
# WEEKDAY GIF
if not os.path.exists('results/GIFs'):
    os.makedirs('results/GIFs')

time_columns = ['access_7week','access_8week','access_9week','access_10week','access_11week','access_12week', 'access_13week','access_14week','access_15week','access_16week','access_17week','access_18week','access_19week', 'access_20week', 'access_21week', 'access_22week']
time_labels = ['7 am', '8 am', '9 am','10 am', '11 am', '12pm', '1 pm','2 pm', '3 pm', '4 pm', '5 pm','6 pm', '7 pm', '8 pm','9 pm', '10 pm']

for time_column, time_label in zip(time_columns, time_labels):
    fig, ax = plt.subplots(figsize=(8, 8))
    plot_accessibility_map(mapping_df, metropolitan_boundary, micropolitan_boundary, time_column, ax)
    ax.set_title(time_label, fontsize=14)
    plt.axis('off')
    plt.savefig(f'results/GIFs/{time_column}.png', bbox_inches='tight')
    plt.close()

# Create GIF from images
images = []
for time_column in time_columns:
    images.append(imageio.imread(f'results/GIFs/{time_column}.png'))

imageio.mimsave('results/GIFs/Week.gif', images, duration=1)

# Remove temp image files
for time_column in time_columns:
    os.remove(f'results/GIFs/{time_column}.png')

# SATURDAY GIF
time_columns_saturday = ['access_7sat', 'access_8sat', 'access_9sat', 'access_10sat', 'access_11sat', 'access_12sat', 'access_13sat', 'access_14sat', 'access_15sat', 'access_16sat', 'access_17sat', 'access_18sat', 'access_19sat', 'access_20sat', 'access_21sat', 'access_22sat']
time_labels_saturday = ['7 am', '8 am', '9 am', '10 am', '11 am', '12 pm', '1 pm', '2 pm', '3 pm', '4 pm', '5 pm', '6 pm', '7 pm', '8 pm', '9 pm', '10 pm']

for time_column, time_label in zip(time_columns_saturday, time_labels_saturday):
    fig, ax = plt.subplots(figsize=(8, 8))
    plot_accessibility_map(mapping_df, metropolitan_boundary, micropolitan_boundary, time_column, ax)
    ax.set_title(time_label, fontsize=14)
    plt.axis('off')
    plt.savefig(f'results/GIFs/{time_column}.png', bbox_inches='tight')
    plt.close()

images_saturday = []
for time_column in time_columns_saturday:
    images_saturday.append(imageio.imread(f'results/GIFs/{time_column}.png'))

imageio.mimsave('results/GIFs/Saturday.gif', images_saturday, duration=1)

for time_column in time_columns_saturday:
    os.remove(f'results/GIFs/{time_column}.png')

# SUNDAY GIF
time_columns_sunday = ['access_7sun', 'access_8sun', 'access_9sun', 'access_10sun', 'access_11sun', 'access_12sun', 'access_13sun', 'access_14sun', 'access_15sun', 'access_16sun', 'access_17sun', 'access_18sun', 'access_19sun', 'access_20sun', 'access_21sun', 'access_22sun']
time_labels_sunday = ['7 am', '8 am', '9 am', '10 am', '11 am', '12 pm', '1 pm', '2 pm', '3 pm', '4 pm', '5 pm', '6 pm', '7 pm', '8 pm', '9 pm', '10 pm']

for time_column, time_label in zip(time_columns_sunday, time_labels_sunday):
    fig, ax = plt.subplots(figsize=(8, 8))
    plot_accessibility_map(mapping_df, metropolitan_boundary, micropolitan_boundary, time_column, ax)
    ax.set_title(time_label, fontsize=14)
    plt.axis('off')
    plt.savefig(f'results/GIFs/{time_column}.png', bbox_inches='tight')
    plt.close()

images_sunday = []
for time_column in time_columns_sunday:
    images_sunday.append(imageio.imread(f'results/GIFs/{time_column}.png'))

imageio.mimsave('results/GIFs/Sunday.gif', images_sunday, duration=1)

for time_column in time_columns_sunday:
    os.remove(f'results/GIFs/{time_column}.png')

Label	Definition	Type
pharmid	unique pharmacy identifier	Text
pharmacy_name	business name	Text
type	chain, independent, grocery, department, hospital	Text
address	pharmacy street address	Text
state	state in which pharmacy is located	Text
lat	GPS coordinates. Lat	Decimal
lon	GPS coordinates. Lon	Decimal
week_open	opening time (24-hr) on weekdays	Integer
week_close	closing time (24-hr) on weekdays	Integer
sat_open	opening time (24-hr) on Saturdays	Integer
sat_close	closing time (24-hr) operations on Saturdays	Integer
sun_hours	opening time (24-hr) on Sundays	Integer
sun_hours	closing time (24-hr) on Sundays	Integer
week_pharm	Typical # pharmacists on weekday	Integer
week_tech	Typical # pharm. techs on weekdays	Integer
sat_pharm	Typical # pharmacists on Saturday	Integer
sat_tech	Typical # pharm. techs on Saturdays	Integer
sun_pharm	Typical # pharmacists on Sundays	Integer
sun_tech	Typical # pharm. techs on Sundays	Integer

	pharmid	pharmacy_name	type	address	state	lat	lon	week_open	week_close	sat_open	sat_close	sun_open	sun_close	geometry
0	VT99	Walgreen Eastern Co., Inc. d/b/a Walgreens #18325	chain	30 Shelburne Shopping Park, Shelburne, VT	VT	44.379262	-73.224989	9.0	19.0	-1.0	-1	-1.0	-1	POINT (-73.22499 44.37926)
1	VT98	WALGREENS #19795	chain	133 N Main St Ste 23, Saint Albans, VT	VT	44.814168	-73.082604	9.0	19.0	9.0	18	9.0	17	POINT (-73.08260 44.81417)
2	VT97	WALGREENS #17747	chain	201 Route 7 S, Milton, VT	VT	44.682189	-73.112042	9.0	19.0	9.0	18	-1.0	-1	POINT (-73.11204 44.68219)
3	VT96	Walgreens #19346	chain	221 Main St, Enosburg Falls, VT	VT	44.905772	-72.806612	9.0	19.0	9.0	18	-1.0	-1	POINT (-72.80661 44.90577)
4	VT95	Walgeens #17183	chain	1184 Prim Rd Ste 2, Colchester, VT	VT	44.537734	-73.247286	9.0	20.0	9.0	18	9.0	17	POINT (-73.24729 44.53773)

Label	Definition	Type	Accuracy	Domain	Missing Data Value(s)	Missing Data Frequency
GEOID	Town unique identifier	Integer	NA	NA	NA	NA
total_pop	Total population in a town	Integer	Defined by Census	0 - ~45,000	NA	NA
subdivision_name	Name of county subdivision (town)	Text	NA	NA	NA	NA
elderly_pop	Total elderly population in a town	Integer	Defined by Census	0 - ~1,000	NA	NA
percent_elderly	Percent elderly residents in a town	Decimal	Defined by Census	0% - ~50%	NA	NA
geometry	Geo data for each town in study area	Text	NA	NA	NA	NA
necta	Classified as metropolitan or micropolitan	Text	Defined by Census	NA	NaN	Frequent because missing data means town is rural

Label	Alias	Definition	Type	Accuracy	Domain	Missing Data Value(s)	Missing Data Frequency
osmid	osmid	unique OpenStreetMap identifier	long integer	n/a	n/a	none	none
highway	highway	highway key from OpenStreetMap with values indicating highway type	text	unknown	residential, tertiary, secondary, primary, unclassified, trunk, etc. See https://wiki.openstreetmap.org/wiki/Key:highway	no value recorded	194,089 of 577,399 edges missing highway type
oneway	oneway	road is for one-way travel only	text representing boolean	unknown	{true, false}	n/a	n/a
reversed	reversed	road direction is opposite the sequence of coordinates	text representing boolean	unknown	{true, false}	n/a	n/a
length	length	length in meters	decimal	geodesic calculation with epsg 4326	unknown	n/a	none

Spatio-Temporal Accessibility of Pharmacy Care in Vermont, USA¶

Authors¶

Abstract¶

Study metadata¶

Study design¶

Background¶

Materials and procedure¶

Computational environment¶

Check Directories¶

Methods¶

Data and variables¶

(1) Pharmacy Dataset¶

Load and Visualize Pharmacy Data¶

(2) Population Dataset¶

Load and Visualize Population Data¶

(3) OpenStreetMap Roads for Vermont and Surroundings¶

Load the Road Network¶

Plot the Road Network¶

Check speed limits and highway types¶

Process the road network¶

Data transformations¶

Analysis¶

Creating catchment areas¶

Service-to-population ratio calculations¶

Results¶

Study Area¶

Figure 1: Mapping town types, population density, and pharmacy locations in VT¶

Hypothesis 1 - Spatial Dimension¶

Figure 2: Weekday accessibility¶

Statistical Significance¶

Hypothesis 2 - Temporal Dimension¶

Figure 3: Accessibility variation by day of the week¶

Figure 4: Accessibility by time of day¶

Figure 5: Accessibility at extreme hours¶

Statistical Significance¶

Hypothesis 3 - Spatio-Temporal Dynamics¶

Figure 7: Time of day and NECTA classification¶

Accessibility GIFs (accessory visuals)¶

Discussion¶

Bias and threats to validity¶

Conclusion¶

Acknowledgements¶

References¶