US traffic accidents spatial autocorrelation

import pandas as pd # to create spatial data import geopandas as gpd # for basemaps import contextily as ctx # For spatial statistics import esda from esda.moran import Moran, Moran_Local import splot from splot.esda import moran_scatterplot, plot_moran, lisa_cluster,plot_moran_simulation import libpysal as lps # Graphics import matplotlib.pyplot as plt import plotly.express as px

df=pd.read_csv('/work/US_Accidents_Dec20_updated.csv')

df.head()

IDobject

Severityint64

A-2716600

A-2716601

A-2716602

A-2716603

A-2716604

This chart is empty

Chart was probably not set up properly in the notebook

df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 1516064 entries, 0 to 1516063
Data columns (total 47 columns):
 #   Column                 Non-Null Count    Dtype  
---  ------                 --------------    -----  
 0   ID                     1516064 non-null  object 
 1   Severity               1516064 non-null  int64  
 2   Start_Time             1516064 non-null  object 
 3   End_Time               1516064 non-null  object 
 4   Start_Lat              1516064 non-null  float64
 5   Start_Lng              1516064 non-null  float64
 6   End_Lat                1516064 non-null  float64
 7   End_Lng                1516064 non-null  float64
 8   Distance(mi)           1516064 non-null  float64
 9   Description            1516064 non-null  object 
 10  Number                 469969 non-null   float64
 11  Street                 1516064 non-null  object 
 12  Side                   1516064 non-null  object 
 13  City                   1515981 non-null  object 
 14  County                 1516064 non-null  object 
 15  State                  1516064 non-null  object 
 16  Zipcode                1515129 non-null  object 
 17  Country                1516064 non-null  object 
 18  Timezone               1513762 non-null  object 
 19  Airport_Code           1511816 non-null  object 
 20  Weather_Timestamp      1485800 non-null  object 
 21  Temperature(F)         1473031 non-null  float64
 22  Wind_Chill(F)          1066748 non-null  float64
 23  Humidity(%)            1470555 non-null  float64
 24  Pressure(in)           1479790 non-null  float64
 25  Visibility(mi)         1471853 non-null  float64
 26  Wind_Direction         1474206 non-null  object 
 27  Wind_Speed(mph)        1387202 non-null  float64
 28  Precipitation(in)      1005515 non-null  float64
 29  Weather_Condition      1472057 non-null  object 
 30  Amenity                1516064 non-null  bool   
 31  Bump                   1516064 non-null  bool   
 32  Crossing               1516064 non-null  bool   
 33  Give_Way               1516064 non-null  bool   
 34  Junction               1516064 non-null  bool   
 35  No_Exit                1516064 non-null  bool   
 36  Railway                1516064 non-null  bool   
 37  Roundabout             1516064 non-null  bool   
 38  Station                1516064 non-null  bool   
 39  Stop                   1516064 non-null  bool   
 40  Traffic_Calming        1516064 non-null  bool   
 41  Traffic_Signal         1516064 non-null  bool   
 42  Turning_Loop           1516064 non-null  bool   
 43  Sunrise_Sunset         1515981 non-null  object 
 44  Civil_Twilight         1515981 non-null  object 
 45  Nautical_Twilight      1515981 non-null  object 
 46  Astronomical_Twilight  1515981 non-null  object 
dtypes: bool(13), float64(13), int64(1), object(20)
memory usage: 412.1+ MB

# we will analyze these columns df=df[['ID','Start_Lat','Start_Lng','City','State','Temperature(F)','Weather_Condition','Sunrise_Sunset']]

df.head()

IDobject

Start_Latfloat64

A-2716600

40.10891

A-2716601

39.86542

A-2716602

39.10266

A-2716603

39.10148

A-2716604

41.06213

cities_by_accident = df.City.value_counts() cities_by_accident

states_by_accident = df.State.value_counts() states_by_accident

# plot top 20 city and states by accidents cities_by_accident[:20].plot(kind='bar')

states_by_accident[:20].plot(kind='bar')

US_states=gpd.read_file('/work/us-state-boundaries.geojson')

US_states

centlatfloat64

-13.9638307 - 63.4205306

intptlatfloat64

-14.2668475 - 63.3461909

18.2164715

18.217648

16.7974379

14.9367835

34.8940208

34.8955256

38.641186

38.6472854

37.1547352

37.1551773

43.9396009

43.9717125

41.5940253

41.5974187

47.4148989

47.4073238

30.902586

30.8634368

46.3164845

46.3158148

US_states.drop(US_states.index[[0,1,18,20,26,36,49]], inplace=True) # drop unnecessary rows

US_states.to_file('/work/US_states.geojson', driver="GeoJSON") # save changed file

US_states=gpd.read_file('/work/US_states.geojson') # read it

US_states.info()

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 49 entries, 0 to 48
Data columns (total 21 columns):
 #   Column     Non-Null Count  Dtype   
---  ------     --------------  -----   
 0   centlat    49 non-null     float64 
 1   intptlat   49 non-null     float64 
 2   intptlon   49 non-null     float64 
 3   gid        49 non-null     int64   
 4   division   49 non-null     int64   
 5   name       49 non-null     object  
 6   region     49 non-null     int64   
 7   oid        49 non-null     object  
 8   funcstat   49 non-null     object  
 9   geoid      49 non-null     object  
 10  state      49 non-null     object  
 11  basename   49 non-null     object  
 12  arealand   49 non-null     object  
 13  stusab     49 non-null     object  
 14  objectid   49 non-null     int64   
 15  mtfcc      49 non-null     object  
 16  areawater  49 non-null     object  
 17  centlon    49 non-null     float64 
 18  statens    49 non-null     object  
 19  lsadc      49 non-null     object  
 20  geometry   49 non-null     geometry
dtypes: float64(4), geometry(1), int64(4), object(12)
memory usage: 8.2+ KB

US_states.crs

US_states=US_states.to_crs(epsg=3857)

#plot it fig, ax=plt.subplots(figsize=(15,15)) US_states.plot(ax=ax, color='black', edgecolor='white', lw=0.5, alpha=0.4 ) ax.axis('off') #no axis # add a basemap ctx.add_basemap(ax,source=ctx.providers.CartoDB.Positron)

# convert accidents .csv into geopandas dataframe from shapely import geometry gdf=gpd.GeoDataFrame(df,crs='EPSG:4326', geometry=gpd.points_from_xy(df.Start_Lng, df.Start_Lat))

gdf=gdf.to_crs(epsg=3857)

gdf.crs

minx, miny, maxx, maxy=gdf.geometry.total_bounds print(minx) print(miny) print(maxx) print(maxy)

-13859005.76344146
2823047.6948689437
-7471003.575963932
6274959.808444523

#sub plot for multilayer maps fig,ax=plt.subplots(1,1, figsize=(15,15)) US_states.plot(ax=ax, color='black', edgecolor='white', alpha=0.5 ) gdf.plot(ax=ax, color='red', markersize=1, alpha=0.2) # use the bounding box coordinates to set the x and y limits ax.set_xlim(minx-1000,maxx+1000) # added/substracted value is to give some margin around total bounds ax.set_ylim(miny-1000,maxy+1000) ax.axis('off') # add a basemap # ctx.add_basemap(ax, source=ctx.providers.CartoDB.Positron)

# Do the spatial join join=gpd.sjoin(gdf,US_states, how='left') join.head()

Unnamed: 0int64

IDobject

A-2716600

A-2716601

A-2716602

A-2716603

A-2716604

accidents_by_states=join.State.value_counts().rename_axis('State').reset_index(name='accidents_count')

accidents_by_states.head()

Stateobject

accidents_countint64

448833

153007

87484

75142

60976

#bar chart of 20 geographies accidents_by_states[:20].plot.bar(figsize=(20,4), x='State', y='accidents_count')

#rename stusab column name to State US_states.rename(columns = {'stusab':'State'}, inplace = True)

US_states.head()

centlatfloat64

intptlatfloat64

34.8940208

34.8955256

38.641186

38.6472854

37.1547352

37.1551773

43.9396009

43.9717125

41.5940253

41.5974187

## join value counts to the US_states US_states=US_states.merge(accidents_by_states, on='State')

US_states.head()

centlatfloat64

intptlatfloat64

34.8940208

34.8955256

38.641186

38.6472854

37.1547352

37.1551773

43.9396009

43.9717125

41.5940253

41.5974187

US_states.to_file('/work/US_states.geojson', driver="GeoJSON")

US_states=gpd.read_file('/work/US_states.geojson')

US_states.head()

centlatfloat64

intptlatfloat64

34.8940208

34.8955256

38.641186

38.6472854

37.1547352

37.1551773

43.9396009

43.9717125

41.5940253

41.5974187

Choropleth map of accidents by states

fig, ax=plt.subplots(figsize=(15,15)) US_states.plot(ax=ax, column='accidents_count', legend=True, alpha=0.8, cmap='RdYlGn_r', scheme='quantiles') ax.axis('off') ax.set_title('US accidents per state', fontsize=22) ctx.add_basemap(ax,source=ctx.providers.CartoDB.Positron)