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


---------------------------------------------------------------------------
OperationalError                          Traceback (most recent call last)
<ipython-input-1-29ee4a320f40> in <module>
      9 pd.set_option("display.max_rows", 1000)
     10 
---> 11 conn = psycopg2.connect('service=firecares-dev')
     12 nfirs = psycopg2.connect('service=nfirs')
     13 parcels = psycopg2.connect(service='parcels')

/opt/conda/lib/python3.7/site-packages/psycopg2/__init__.py in connect(dsn, connection_factory, cursor_factory, **kwargs)
    124 
    125     dsn = _ext.make_dsn(dsn, **kwargs)
--> 126     conn = _connect(dsn, connection_factory=connection_factory, **kwasync)
    127     if cursor_factory is not None:
    128         conn.cursor_factory = cursor_factory

OperationalError: definition of service "firecares-dev" not found

People per structure used for sanity check on parcel counts by department


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)

Get owned census-tracts for department and count parcels


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))

Structure count by hazard category by owned census tract geometries


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)

Export of parcels by owned census tracts


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')