In [ ]:
# This notebook assumes to be running from your FireCARES VM (eg. python manage.py shell_plus --notebook --no-browser)

import sys
import os
import time
import pandas as pd
import numpy as np
import psycopg2
import django
sys.path.insert(0, os.path.realpath('..'))
import folium
django.setup()
from IPython.display import display
from django.contrib.gis.geos import GEOSGeometry

pd.set_option("display.max_rows",100)

def display_geom(geom):
    _map = folium.Map(location=[geom.centroid.y, geom.centroid.x],
                      tiles='Stamen Toner')
    _map.choropleth(geo_str=geom.geojson, line_weight=0, fill_opacity=0.2, fill_color='green')
    ll = geom.extent[1::-1]
    ur = geom.extent[3:1:-1]
    _map.fit_bounds([ll, ur])

    return _map

nfirs = psycopg2.connect('service=nfirs')
fc = psycopg2.connect('service=firecares')
parcels = psycopg2.connect('service=parcels')

Total departments/stations in FireCARES


In [ ]:
q = """
select *
from
    (select count(1) as total_departments from firestation_firedepartment) total_departments,
    (select count(1) as total_stations from firestation_firestation) total_stations
"""

pd.read_sql_query(q, fc)

Covered population


In [ ]:
q = """
select sum(fd.population) as total_population
from firestation_firedepartment fd
"""

pd.read_sql_query(q, fc)

Covered area


In [ ]:
from firecares.firestation.models import FireDepartment

deps = FireDepartment.objects.filter(population__isnull=False).exclude(population=0)

areas = filter(lambda x: x[1] is not None, map(lambda x: (x.id, x.geom_area), deps))

In [ ]:
sum(a[1] for a in areas)

In [ ]:
len(areas)

Total incidents and building fires


In [ ]:
display(pd.read_sql_query("select count(1) as buildingfires from buildingfires", nfirs))
display(pd.read_sql_query("select count(1) as incidents from incidentaddress", nfirs))

Coverage of departments based on department-owned census geometry that have populations


In [ ]:
# Number of departments with owned census tracts

q = """
select *, (null_tracts_count + not_null_tracts_count) as total
from
    (select count(1) as null_tracts_count from firestation_firedepartment where owned_tracts_geom is null) null_owned_tracts,
    (select count(1) as not_null_tracts_count from firestation_firedepartment where owned_tracts_geom is not null) not_null_owned_tracts
    """

print "Owned tracts geom"
display(pd.read_sql_query(q, fc))


q = """
select *, (null_jurisdictions + not_null_jurisdictions) as total
from
    (select count(1) as null_jurisdictions from firestation_firedepartment where geom is null) null_jurisdiction,
    (select count(1) as not_null_jurisdictions from firestation_firedepartment where geom is not null) not_null_jurisdiction
"""
print "Jurisdictions"
display(pd.read_sql_query(q, fc))

In [ ]:
%%time

q = """
select ST_Union(ST_Union(fd.owned_tracts_geom), ST_Union(fd.geom))
from firestation_firedepartment fd
where fd.state = 'MO'
"""

df = pd.read_sql_query(q, fc)

poly = GEOSGeometry(df.values[0][0])
display(display_geom(poly.simplify()))

In [ ]:
%%time
# Total square meters in the USA (territories excluded)

us = 9147593069344.

pd.set_option('display.float_format', lambda x: '%f' % x)

q = """
select ST_Area(geography(ST_Union(ST_Union(ST_MakeValid(fd.owned_tracts_geom)), ST_Union(ST_MakeValid(fd.geom)))))
from firestation_firedepartment fd
where fd.population is not null and fd.population != 0
and fd.state = 'CT'
"""

display(pd.read_sql_query(q, fc).values[0][0] / us)

Invalid geometries


In [ ]:
%%time

q = """
select count(1), state from firestation_firedepartment fd where not ST_IsValid(fd.geom) group by fd.state
"""

display(pd.read_sql_query(q, fc))

In [ ]:
%%time

q = """
select count(1), state from firestation_firedepartment fd where not ST_IsValid(fd.owned_tracts_geom) group by fd.state
"""

display(pd.read_sql_query(q, fc))

Bad address geocodes


In [ ]:
%%time

q = """
select id, ST_SRID(a.geom), ST_X(a.geom), ST_Y(a.geom), a.state_province
from firecares_core_address a
where ST_X(a.geom) > 180 or ST_X(a.geom) < -180 or ST_Y(a.geom) > 180 or ST_Y(a.geom) < -180
"""

display(pd.read_sql_query(q, fc))

In [ ]:
# Fix bad addresses
%%time

q = """
update firecares_core_address
set geom = ST_Transform(ST_SetSRID(geom, 3857), 4326)
where ST_X(geom) > 180 or ST_X(geom) < -180 or ST_Y(geom) > 180 or ST_Y(geom) < -180
"""

c = fc.cursor()
c.execute(q)
fc.commit()

Bad addresses (missing country/etc)


In [ ]:
%%time

q = """
select count(1)
from firecares_core_address
where address_line2 = 'None'
"""

display(pd.read_sql_query(q, fc))

Verify station lat/lons


In [ ]:
from firecares.firestation.models import FireStation
map(lambda x: (x.get('geom').coords, x.get('name')), FireStation.objects.filter(department_id=97963).values('geom', 'name'))

Jurisdiction boundary statistics


In [ ]:
%%time

from firecares.firestation.models import FireDepartment
from firecares.utils import lenient_summation
import json
import os
d = []
fname = '/firecares/outf.json'

if os.path.exists(fname):
    d = json.load(open(fname))
ids = map(lambda x: x.get('fcid'), d)

for idx, fd in enumerate(FireDepartment.objects.defer('owned_tracts_geom').filter(archived=False)):
    if fd.id in ids:
        continue
    if idx % 10 == 0:
        print idx
    with open(fname, 'w') as f:
        json.dump(d, f)
    d.append({'fcid': fd.id,
              'name': fd.name,
              'state': fd.state,
              'boundary': True if fd.geom else False,
              'population': fd.population,
              'region': fd.region,
              'fires_count': lenient_summation(*map(lambda x: x.get('count'), fd.metrics.residential_structure_fire_counts['all'])),
              'casualty_count': fd.metrics.nfirs_deaths_and_injuries_sum['all'],
              'station_count': fd.firestation_set.count()})

print 'done'

In [ ]:
deduped = []
for i in d:
    if i.get('fcid') in map(lambda x: x.get('fcid'), deduped):
        continue
    deduped.append(i)

In [ ]:
pd.DataFrame(deduped).to_csv('/tmp/stations.csv')