<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Validated boundaries to government unit incident density comparison
The backing theory for this notebook is proving that we will be able to use the highest-density (fire count vs government unit area) government unit to determine a department's boundary for departments that do not have boundaries.
Step5:  DB migration/setup
Step8:  Processing
Step9:  Results
<ASSISTANT_TASK:>
Python Code:
import psycopg2
from psycopg2.extras import RealDictCursor
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
#pd.set_option('display.float_format', lambda x: '%.3f' % x)
%matplotlib inline
conn = psycopg2.connect('service=firecares')
nfirs = psycopg2.connect('service=nfirs')
Explanation: Validated boundaries to government unit incident density comparison
The backing theory for this notebook is proving that we will be able to use the highest-density (fire count vs government unit area) government unit to determine a department's boundary for departments that do not have boundaries.
End of explanation
# Create materialized view of all usgs govt units in FireCARES
q = 
create materialized view if not exists usgs_governmentunits as
(
    select id, population, county_name as name, 'countyorequivalent' as source, geom from usgs_countyorequivalent where geom is not null
    union
    select id, population, place_name as name, 'incorporatedplace' as source, geom from usgs_incorporatedplace where geom is not null
    union
    select id, population, minorcivildivision_name as name, 'minorcivildivision' as source, geom from usgs_minorcivildivision where geom is not null
    union
    select id, population, name, 'nativeamericanarea' as source, geom  from usgs_nativeamericanarea where geom is not null
    union
    select id, 0 as population, name, 'reserve' as source, geom  from usgs_reserve where geom is not null
    union
    select id, population, state_name as name, 'stateorterritoryhigh' as source, geom from usgs_stateorterritoryhigh where geom is not null
    union
    select id, population, place_name as name, 'unincorporatedplace' as source, geom from usgs_unincorporatedplace where geom is not null
);
create unique index on usgs_governmentunits (id, source);
create index on usgs_governmentunits using gist (geom);
    
with conn.cursor() as c:
    c.execute(q)
    
conn.commit()
# Link remote firecares usgs_governmentunits view to nfirs-local usgs_government units
q = 
create foreign table usgs_governmentunits (id integer, population integer, name character varying(120), source text, geom geometry)
server firecares
options (table_name 'usgs_governmentunits');
with nfirs.cursor() as c:
    c.execute(q)
    
nfirs.commit()
# Old nfirs.firestation_firedepartment foreign table columns needed to be synced
q = 
alter foreign TABLE firestation_firedepartment add column archived boolean NOT NULL;
alter foreign TABLE firestation_firedepartment add column domain_name character varying(255);
alter foreign TABLE firestation_firedepartment add column owned_tracts_geom public.geometry(MultiPolygon,4326);
alter foreign TABLE firestation_firedepartment add column display_metrics boolean NOT NULL;
alter foreign TABLE firestation_firedepartment add column boundary_verified boolean NOT NULL;
alter foreign TABLE firestation_firedepartment add column cfai_accredited boolean NOT NULL;
alter foreign TABLE firestation_firedepartment add column ems_transport boolean NOT NULL;
alter foreign TABLE firestation_firedepartment add column staffing_verified boolean NOT NULL;
alter foreign TABLE firestation_firedepartment add column stations_verified boolean NOT NULL;
alter foreign TABLE firestation_firedepartment add column census_override boolean NOT NULL;
alter foreign TABLE firestation_firedepartment add column additional_fdids character varying(255);
with nfirs.cursor() as c:
    c.execute(q)
nfirs.commit()
q = 
create foreign table if not exists firecares_core_address (id integer NOT NULL,
    address_line1 character varying(100) NOT NULL,
    address_line2 character varying(100),
    city character varying(50) NOT NULL,
    state_province character varying(40) NOT NULL,
    postal_code character varying(10) NOT NULL,
    geom public.geometry(Point,4326),
    geocode_results text,
    country_id character varying(2) NOT NULL)
server firecares
options (table_name 'firecares_core_address');
with nfirs.cursor() as c:
    c.execute(q)
nfirs.commit()
Explanation: DB migration/setup
End of explanation
q = 
select id, fdid, state, name
from firestation_firedepartment
where boundary_verified = true;
with nfirs.cursor(cursor_factory=RealDictCursor) as c:
    c.execute(q)
    fds = c.fetchall()
q = 
with fires as (select * from joint_buildingfires
    inner join joint_incidentaddress
    using (fdid, inc_no, inc_date, state, exp_no)
    where state = %(state)s and fdid = %(fdid)s
),
govt_units as (
    select gu.name, gu.source, gu.id, gu.geom, fd.id as fc_id, fd.geom as fd_geom, ST_Distance(addr.geom, ST_Centroid(gu.geom)) as distance_to_headquarters
    from firestation_firedepartment fd
    inner join firecares_core_address addr
        on addr.id = fd.headquarters_address_id
    join usgs_governmentunits gu
    on ST_Intersects(ST_Buffer(addr.geom, 0.05), gu.geom)
    where 
        fd.fdid = %(fdid)s and fd.state = %(state)s and source != 'stateorterritoryhigh'
    )
select gu.fc_id, count(fires) / ST_Area(gu.geom) as density, count(fires), ST_Area(ST_SymDifference(gu.fd_geom, gu.geom)) / ST_Area(gu.fd_geom) as percent_difference_to_verified_boundary, ST_Area(gu.geom), gu.distance_to_headquarters, gu.name, gu.id, gu.source from fires
inner join govt_units gu
on ST_Intersects(fires.geom, gu.geom)
group by gu.name, gu.id, gu.geom, gu.source, gu.distance_to_headquarters, gu.fd_geom, gu.fc_id
order by ST_Area(gu.geom) / count(fires) asc;
for fd in fds:
    with nfirs.cursor(cursor_factory=RealDictCursor) as c:
        print 'Analyzing: {} (id: {} fdid: {} {})'.format(fd['name'], fd['id'], fd['fdid'], fd['state'])
        c.execute(q, dict(fdid=fd['fdid'], state=fd['state']))
        items = c.fetchall()
        df = pd.DataFrame(items)
        df.to_csv('./boundary-analysis-{}.csv'.format(fd['id']))
Explanation: Processing
End of explanation
from glob import glob
df = None
for f in glob("boundary-analysis*.csv"):
    if df is not None:
        df = df.append(pd.read_csv(f))
    else:
        df = pd.read_csv(f)
df.rename(columns={'Unnamed: 0': 'rank'}, inplace=True)
selected_government_units = df[df['rank'] == 0].set_index('fc_id')
total_validated_department_count = len(selected_government_units)
perfect_fits = len(selected_government_units[selected_government_units['percent_difference_to_verified_boundary'] == 0])
print 'Perfect fits: {}/{} ({:.2%})'.format(perfect_fits, total_validated_department_count, float(perfect_fits) / total_validated_department_count)
print 'Machine-selected government unit area difference mean: {:.2%}'.format(df[df['rank'] == 0].percent_difference_to_verified_boundary.mean())
selected_government_units['percent_difference_to_verified_boundary'].hist(bins=50)
selected_government_units
df.set_index('fc_id')
df.to_csv('./validated-boundary-vs-government-unit-density.csv')
pd.read_csv('./validated-boundary-vs-government-unit-density.csv')
Explanation: Results
End of explanation
<END_TASK>