import seaborn as sns
import metapack as mp
import pandas as pd
import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display 
from shapely.geometry import Point

%matplotlib inline
sns.set_context('notebook')


pkg = mp.jupyter.open_package()
#pkg = mp.jupyter.open_source_package()
pkg

pkg.resource('collisions')

col = pkg.resource('collisions').read_csv()


col_sd = col[col.juris.isin( ['3711','3714','3797', '3725'] )]
len(col_sd)


## Create a new GeoPandas frame, converting the targetlongitude and targetlatitude
## colums to a Shapely Point and assigning it to the frame's geometry

t = col_sd.dropna(subset=['longitude', 'latitude'])

gdf = gpd.GeoDataFrame(t, geometry=
                        [Point(x,y) for x,y in zip(-t.longitude, t.latitude) ])

len(gdf),len(col_sd)

"""
Explanation: Traffic Accidents by San Diego Community
This analysis links traffic accident records from SWITRS to San Diego planning communities to see which communityies have the most traffoc accidents. Unfortunately, the data, from the years 2014 and 2015, has very poor lat/lon values -- most records have none -- so the number of records that can be linked into communities is too small for a worthwhile answer
End of explanation
"""


len(col)

comm  = pkg.reference('communities').geoframe()
ax = comm.plot()
gdf.plot(ax=ax, color='red')

"""
Explanation: There are only 106 records with valid lat/lon, out of the original 8919, which is a pretty poort sample. Let's plot those points over a map of the county to see where they are.
End of explanation
"""


comm = pkg.reference('communities').geoframe()
comm.head()

sd_comm = comm[ comm['type'] == 'sd_community']
ax = sd_comm.plot()
gdf.plot(ax=ax, color='red')

"""
Explanation: Not only are there only 106 valid points, but a lot of the points are in crazy places, outside of the city, outside of the county, and in the ocean.
End of explanation
"""


j = gpd.sjoin(gdf.dropna(subset=['geometry']), sd_comm)
sum(j.code.notnull())

"""
Explanation: Spatially join the accident points with the communities.
End of explanation
"""


j.groupby('name').count().case_id.sort_values(ascending=False)

col.head().T

# SD County
col_sdc = col[ ( (col.juris.str.startswith('37')) |  (col.juris == '9645'))] # 9645 is CHP in San Diego City
len(col_sdc)

col_sdc.juris.value_counts().sort_index()


"""
Explanation: Out of the small number with real points, only 84 of them link to a community in San Diego, we can get an answer, but it will not be a good one.
End of explanation
"""


col[col.juris == '3720']

pkg.resource('sd_juris')

juris = pkg.resource('sd_juris').read_csv()
#juris['code'] = juris.code.astype('str')

t = col_sdc.merge(juris, left_on='cnty_city_loc', right_on='code')
df = t.groupby('name').count().case_id.sort_values(ascending=False).to_frame()
df.columns = ['accidents']
df = df.join(juris.set_index('name'))
df['area'] = df.area / 1_000_000 # sq meters to sq km

# Rate by area and population 
df['area_rate'] = df.accidents / df.area
df['pop_rate'] = df.accidents / df.population * 100_000
df['dvmt_rate'] = df.accidents / df.daily_vmt
df['road_rate'] = df.accidents / df.road_miles

df.sort_values('pop_rate', ascending=False)

from tabulate import tabulate
print(tabulate(df[['accidents','area_rate','pop_rate']], headers='keys', tablefmt='pipe'))

from tabulate import tabulate
print(tabulate(df[['dvmt_rate', 'road_rate']].sort_values('road_rate', ascending=False), headers='keys', tablefmt='pipe'))

"""
Explanation: What's interesting here is that there are a lot of city codes that are in the jurisdiction documentation that aren't in this set: 3780 Poway, 3720 Lemon grove. From the look for the documentation, these values might be mapped to 3700.
End of explanation
"""


