# 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
gdf=gpd.read_file('/work/census tract.geojson')
gdf.head() # show first 5 rows
gdf.info()
gdf.crs
# trim the data to the bare minimum columns
gdf=gdf[['geoid','B01003001','geometry']]
# rename the columns
gdf.columns=['FIPS','TotalPopu','geometry']
gdf.tail()
gdf=gdf.drop(625)
# delete last row which is for the entire city of detroit
gdf['FIPS']=gdf['FIPS'].str.replace('15000US','')
gdf.tail()
gdf.sort_values(by='TotalPopu').head(20)
# sort by total pop
# delete less than 100 population geographies
gdf = gdf[gdf['TotalPopu']>100]
# get the layers into a web mercator projection
# reproject to web mercator
gdf=gdf.to_crs(epsg=3857)
gdf.crs
# plot it
fig, ax=plt.subplots(figsize=(15,15))
gdf.plot(ax=ax,
color='black',
edgecolor='white',
lw=0.5,
alpha=0.4
)
ax.axis('off') #no axis
#add basemap
ctx.add_basemap(ax, source=ctx.providers.CartoDB.Positron)
GI = gpd.read_file('/work/Green_Stormwater_Infrastructure_Locations_OpenData.geojson')
GI.head() # show first 5 rows
GI.info() # show columns and data types
# trim the data to the bare minimum columns
GI = GI[['OBJECTID','PRACTICE_TYPE','GALLONS_MANAGED_MG','PRACTICE_AREA','DISTRICT','geometry']]
# last 5 rows
GI.tail()
GI.crs
GI=GI.to_crs(epsg=3857)
minx1, miny2, maxx1, maxy2=gdf.geometry.total_bounds
print(minx1)
print(miny2)
print(maxx1)
print(maxy2)
GI.crs
#plot it
fig, ax=plt.subplots(figsize=(12,12))
GI.plot(ax=ax,
color='red',
markersize=5
)
ax.axis('off') #no axis
# add a basemap
ctx.add_basemap(ax,source=ctx.providers.CartoDB.Positron)
#create a two layer map
# get the bounding box for GI
minx, miny, maxx, maxy=GI.geometry.total_bounds
print(minx)
print(miny)
print(maxx)
print(maxy)
#sub plot for multilayer maps
fig, ax=plt.subplots(1,1,figsize=(12,12))
gdf.plot(ax=ax,
color='gray',
edgecolor='black',
alpha=0.5
)
GI.plot(ax=ax,
color='black',
markersize=10
)
# 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') #no axis
# add a basemap
ctx.add_basemap(ax,source=ctx.providers.CartoDB.Positron)
# do the spatial join
join=gpd.sjoin(GI, gdf, how='left')
join.head()
GI_by_gdf=join.FIPS.value_counts().rename_axis('FIPS').reset_index(name='GI_count')
GI_by_gdf.head()
# bar chart of 20 geographies
GI_by_gdf[:20].plot.bar(figsize=(20,4),
x='FIPS',
y='GI_count')
# join GI value counts back to the gdf
gdf=gdf.merge(GI_by_gdf, on='FIPS')
gdf.head()
# save this geodataframe
gdf.to_file('/work/gdf.geojson', driver="GeoJSON")
#read saved file
gdf=gpd.read_file('/work/gdf.geojson')
gdf.head()
gdf.info()
fig,ax=plt.subplots(figsize=(15,15))
gdf.plot(ax=ax,
column='GI_count',
legend=True,
alpha=0.8,
edgecolor='black',
cmap='RdYlGn_r', # a diverging color scheme
scheme='quantiles')
ax.axis('off')
ax.set_title('GI count per census tract', fontsize=22)
ctx.add_basemap(ax,source=ctx.providers.CartoDB.Positron)