In [ ]:

    

# This notebook assumes to be running from your local machine at this folder location using a Python 2 kernel,
# with the "nfirs" and "firecares-dev" postgres service configured in your pg_service.conf, pointed at the
# NIFRS database as well as the FireCARES database instance
# see requirements.txt in the root of the repository for python dependencies

import sys
import os
import time
import pandas as pd
import numpy as np
import psycopg2
import folium
import matplotlib as plt
from shapely import wkb
from shapely.geometry import mapping
from pretty import pprint
from IPython.display import display
#pd.set_option("display.max_rows",100)

def display_df_geom(df, col, row=0):
    return display_geom(wkb.loads(df[col][row], hex=True))

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

nfirs = psycopg2.connect("service=nfirs")
fc = psycopg2.connect("service=firecares-dev")


    

In [ ]:

    

# Pull featured departments that have a computed census tract boundary

df = pd.read_sql_query("select id, name, owned_tracts_geom from firestation_firedepartment where featured and not archived;", fc)
print 'TOTAL NON-ARCHIVED FEATURED: {}'.format(len(df))

print 'DEPARTMENTS MISSING CENSUS TRACT JURISDICTION'
display(df[df.owned_tracts_geom.isnull()][['id', 'name']])

df2 = pd.read_sql_query("select id, name, state, owned_tracts_geom from firestation_firedepartment where featured and not archived and owned_tracts_geom is not null;", fc)
records = df2.to_dict(orient='records')

print 'TOTAL NON-ARCHIVED FEATURED WITH CENSUS GEOM {}'.format(len(records))


    

In [ ]:

    

# Create holding table
create = """
create table if not exists firedepartment_parcel_unions (
    fd_id integer not null,
    parcel_union_geom geometry(MultiPolygon, 4326),
    owned_tracts_geom geometry(MultiPolygon, 4326)
);

create unique index on firedepartment_parcel_unions (fd_id);
"""
with nfirs.cursor() as c:
    c.execute(create)
    nfirs.commit()


    

In [ ]:

    

PARCEL_COVERAGE = """
select ST_Union(p.wkb_geometry)
from parcel_risk_category_local p
where ST_Intersects(ST_SetSRID(%(owned_geom)s::geometry, 4326), p.wkb_geometry)
"""

records = pd.read_sql_query("select id, name, state, owned_tracts_geom from firestation_firedepartment where featured and not archived and owned_tracts_geom is not null;", fc).to_dict(orient='records')

filtered = records # filter(lambda x: x.get('id') == 95982, records)

def already_done(fd_id):
    with nfirs.cursor() as c:
        c.execute("select fd_id from firedepartment_parcel_unions where fd_id = %(fd_id)s;", {'fd_id': fd_id})
        return c.fetchone() is not None

for r in filtered:
    if already_done(r.get('id')):
        continue
    print 'Working on {}'.format(r.get('id'))
    jurisdiction = wkb.loads(r.get('owned_tracts_geom'), hex=True)
    df = pd.read_sql_query(PARCEL_COVERAGE, nfirs, params={'owned_geom': jurisdiction.wkb_hex})
    parcels = wkb.loads(df['st_union'][0], hex=True)
    #_map = display_geom(parcels)
    #_map.choropleth(geo_str=mapping(jurisdiction), fill_opacity=0, line_weight=10, line_color='red')
    #display(_map)
    
    with nfirs.cursor() as c:
        c.execute("insert into firedepartment_parcel_unions values (%(fd_id)s, ST_SetSRID(%(parcels)s::geometry, 4326), ST_SetSRID(%(tracts)s::geometry, 4326));",
                  {
                      'fd_id': r.get('id'),
                      'parcels': parcels.wkb_hex,
                      'tracts': jurisdiction.wkb_hex
                  })
        nfirs.commit()
    
    print 'DEPT {} - parcel coverage: {} census jurisdiction: {} % coverage: {}'.format(r.get('id'),
                                                                                        parcels.area,
                                                                                        jurisdiction.area,
                                                                                        parcels.area / jurisdiction.area)


    

In [ ]:

    

%matplotlib inline
q = """
select fd_id, ST_Area(parcel_union_geom) / ST_Area(owned_tracts_geom) as parcel_percent_coverage from firedepartment_parcel_unions;
"""

df = pd.read_sql_query(q, nfirs)
df[['parcel_percent_coverage']].hist()


    

In [ ]:

    

df.sort(['parcel_percent_coverage'])


    

In [ ]:

    

q = """
select fd_id, parcel_union_geom, owned_tracts_geom from firedepartment_parcel_unions where fd_id = %(fd_id)s;
"""

df = pd.read_sql_query(q, nfirs, params={'fd_id': 80595})

d = df.to_dict(orient='records')[0]

jurisdiction = wkb.loads(d.get('owned_tracts_geom'), hex=True)
parcels = wkb.loads(d.get('parcel_union_geom'), hex=True)

#_map = display_geom(parcels)
#_map.choropleth(geo_str=mapping(jurisdiction), fill_opacity=0, line_weight=2, line_color='red')
_map = display_geom(jurisdiction)
display(_map)