get_nhs_miles_by_county_and_type.py

import geopandas as gpd
import pandas as pd
import fiona

###
#read in NHS data
###
#download data: #http://www.dot.ca.gov/hq/tsip/gis/datalibrary/zip/highway/NHS_Map21.zip
nhs = gpd.read_file("data/NHS_Map21/NHS_Map21.shp")

###
#read tomtom county table
###
#download data: https://mtcdrive.box.com/s/wo86uw3w9m7c189sxk3ofumgw4013pw0
cnty = gpd.read_file('C:/projects/merge_highways/TomTom_CA.gdb/TomTom_CA.gdb', driver="OpenFileGDB", layer="mn_a8")

cnty = cnty.to_crs({'init': 'epsg:4269'})

nhs_cnty = gpd.sjoin(nhs, cnty, how="inner", op='intersects')

field_select = ['CT_DISTRIC', 'FULLNAME', 'HWY_NUM', 'NHS', 'RAMP', 'RAM_ROUTE1',
       'RAM_ROUTE_', 'SHIELD', 'Type', 'geometry', 'index_right',
       'ID', 'NAME', 'NAMELC', 'ORDER00','ORDER08',
       'ORDER01', 'POP', 'POPCLASS']

nhs_cnty = nhs_cnty.loc[:,field_select]

# nhs_cnty.to_file('nhs_cnty.shp')

county_select = ["CA041", "CA055", "CA095", "CA097", "CA081", "CA001", "CA085", "CA013", "CA075"]

nhs_cnty_ba = nhs_cnty.loc[nhs_cnty.ORDER08.isin(county_select),:]

nhs_cnty_ba = nhs_cnty_ba.to_crs({'init': 'epsg:26910'})

nhs_cnty_ba["length"] = nhs_cnty_ba["geometry"].length

meter_in_miles = 0.000621

nhs_cnty_ba["length_miles"] = nhs_cnty_ba["length"]*meter_in_miles

cnty_grpd = nhs_cnty_ba.groupby(["NAME","Type"])["length_miles"].sum()

cnty_grpd.unstack(level=-1).to_excel("nhs_lengths_by_county_and_type.xls")