In [1]:
import psycopg2
import pandas as pd
import folium
from shapely import wkb
from shapely.geometry import mapping
from tabulate import tabulate
from IPython.display import display
pd.set_option("display.max_rows", 1000)
conn = psycopg2.connect('service=firecares-dev')
nfirs = psycopg2.connect('service=nfirs')
parcels = psycopg2.connect(service='parcels')
def display_geom(geom):
_map = folium.Map(location=[geom.centroid.y, geom.centroid.x],
tiles='Stamen Toner')
_map.choropleth(geo_str=mapping(geom), line_weight=0, fill_opacity=0.2, fill_color='green')
ll = geom.bounds[1::-1]
ur = geom.bounds[3:1:-1]
_map.fit_bounds([ll, ur])
return _map
In [ ]:
q = """select fd.id, fd.name, fd.state,
COALESCE(fd.population, 0) as population,
sum(rm.structure_count) as structure_count,
fd.population / sum(rm.structure_count)::float as people_per_structure
from firestation_firedepartment fd
inner join firestation_firedepartmentriskmodels rm
on rm.department_id = fd.id
where rm.level != 0
group by fd.id, COALESCE(fd.population, 0)
order by COALESCE(population / sum(rm.structure_count)::float, 0) desc"""
df = pd.read_sql_query(q, conn)
df.to_csv('/tmp/people_per_structure_by_department.csv')
filtered = df[df['population'] > 100000][:30]
res = tabulate(filtered, headers='keys', tablefmt='pipe')
open('/tmp/outf.md', 'w').write(res)
! cat /tmp/outf.md | pbcopy
display(filtered)
In [ ]:
q = """SELECT ST_Multi(ST_Union(bg.geom))
FROM nist.tract_years ty
INNER JOIN census_block_groups_2010 bg
ON ty.tr10_fid = ('14000US'::text || "substring"((bg.geoid10)::text, 0, 12))
WHERE ty.fc_dept_id = %(id)s
GROUP BY ty.fc_dept_id"""
geom = pd.read_sql_query(q, nfirs, params={'id': 96649})['st_multi'][0]
In [ ]:
display_geom(wkb.loads(geom, hex=True))
In [ ]:
q = """select count(1), risk_category
from parcel_risk_category_local p
where ST_Intersects(p.wkb_geometry, ST_SetSRID(%(owned_geom)s::geometry, 4326))
group by risk_category"""
df = pd.read_sql_query(q, nfirs, params={'owned_geom': geom})
display(df)
In [61]:
q = """select parcel_id, risk_category
from parcel_risk_category_local l
where ST_Intersects(l.wkb_geometry, ST_SetSRID(%(owned_geom)s::geometry, 4326))"""
owned_parcels = pd.read_sql_query(q, nfirs, params={'owned_geom': geom})
In [70]:
import geopandas
q = """select p.*, rc.risk_category as hazard_level from parcels p
inner join parcel_risk_category_local rc using (parcel_id)
where p.parcel_id in %(ids)s"""
res = map(lambda x: x[0], owned_parcels.values)
gdf = geopandas.read_postgis(q, nfirs, geom_col='wkb_geometry', params={'ids': tuple(res)})
gdf.drop('risk_category', 1)
gdf.crs = {'init': 'epsg:4326'}
gdf.to_file('/tmp/tamarac-parcels.shp', driver='ESRI Shapefile')