In [ ]:

    

import psycopg2
import pandas as pd
import geopandas as gpd
from shapely import wkb
from shapely.geometry import mapping as to_geojson
import folium

pd.options.display.max_columns = None
pd.options.display.max_rows = None


    

In [ ]:

    

# Assumes that you have a 'firecares' service setup in your ~/.pg_service.conf w/ SELECT access to 'firestation_firedepartment'
conn = psycopg2.connect('service=firecares')


    

In [ ]:

    

q = """
select id, name, state, fdid, population, geom, owned_tracts_geom
from firestation_firedepartment
where boundary_verified = true
    and owned_tracts_geom is not null
"""
df = gpd.read_postgis(q, conn, crs=4326)
df['owned_tracts_geom'] = df['owned_tracts_geom'].apply(lambda x: wkb.loads(x, hex=True))


    

In [ ]:

    

q = """
select
    id,
    name,
    state,
    fdid,
    ST_Area(geom) as barea,
    ST_Area(owned_tracts_geom) as oarea,
    ST_Area(ST_SymDifference(geom, owned_tracts_geom)) as diff,
    ST_Area(ST_SymDifference(geom, owned_tracts_geom)) / ST_Area(geom) * 100 as perdiff_vs_boundaries
from firestation_firedepartment
where boundary_verified = true
    and owned_tracts_geom is not null
order by perdiff_vs_boundaries desc;
"""

diffdf = pd.read_sql(q, conn)
display(diffdf)
print 'Mean % difference in validated boundary vs owned tracts {}%'.format(diffdf.perdiff_vs_boundaries.mean())


    

In [ ]:

    

def show_department_polys(_id):
    m = folium.Map(tiles='Stamen Toner')
    boundary = df[df.id == _id]['geom']
    folium.GeoJson(
        to_geojson(boundary),
        name='Validated boundary',
        style_function=lambda x: {'fillColor': '#0000ff', 'weight': 2, 'color': '#0000ff'}
    ).add_to(m)
    
    owned_tracts = df[df.id == _id]['owned_tracts_geom']
    
    folium.GeoJson(
        to_geojson(gpd.GeoSeries(owned_tracts)),
        name='Owned census tracts',
        style_function=lambda x: {'fillColor': '#00cc00', 'weight': 0.4, 'color': '#00cc00'}
    ).add_to(m)
    
    folium.LayerControl().add_to(m)
    
    extents_geom = owned_tracts.values[0] if owned_tracts.values[0].area > boundary.values[0].area else boundary.values[0]
    
    x1, y1, x2, y2 = extents_geom.bounds
    m.fit_bounds([(y1, x1), (y2, x2)])
    
    display(m)


    

In [ ]:

    

show_department_polys(91743)