<jupyter_start><jupyter_text>Construct a node-edge network of all the intersections within our bounding box. This allows us to compute distance calulations and aggregates with point data on this network<jupyter_code>net = osm.pdna_network_from_bbox(45.1453, -123.1095, 45.7859, -122.2575, network_type='drive')
# precompute a distance of 3000 meters, we won't be doing calulations further than that
net.precompute(3000)<jupyter_output><empty_output><jupyter_text>Preprocess point of interest data and add to the networkWe can 'set' points of interest to the network (e.g. grocery stores) to later perform aggregations on. We do this by simply passing vectors of lats and longs, a label for the points and the max distance we care about in future calulations. Some of the data is in ESRI Shapefile format so they need to be loaded and have their coordinates extracted.<jupyter_code>def process_pois(rootdir, network, max_map=250):
    # iterate through point shapefiles
    for data in glob(join('./data/pois', '*')):
        print('Processing {}...'.format(data))
        df = gpd.read_file(data)
        # confirm the projection is WGS84
        if df.crs != {'init': 'epsg:4326'}:
            df = df.to_crs({'init': 'epsg:4326'})
        xs = df.apply(lambda row: row['geometry'].x, axis=1)
        ys = df.apply(lambda row: row['geometry'].y, axis=1)
        
        network.set(network.get_node_ids(xs, ys, max_map), variable=None, name=data.split('/')[3])
    print('Done!')
process_pois('./data/pois', net)<jupyter_output>Processing ./data/pois/grocery_store...
Processing ./data/pois/lrt_stop...
Processing ./data/pois/tran_cen...
Done!<jupyter_text>The OSM data already has lat longs but needs to be separated into groups (i.e. unconcatenated)<jupyter_code># OSM already has lat longs so process these can essentially be added to the network directly
osm = pd.read_csv('./data/osm.csv')
osm.head(2)
osm.alias.unique()<jupyter_output><empty_output><jupyter_text>Separate out different types of amenities<jupyter_code>fd = osm[osm['alias'] == 'food_drink'][['id', 'lon', 'lat']]
fd.index.name = 'fd'
cafe = osm[osm['alias'] == 'cafe'][['id', 'lon', 'lat']]
cafe.index.name = 'cafe'
college = osm[osm['alias'] == 'college'][['id', 'lon', 'lat']]
college.index.name = 'college'
ent = osm[osm['alias'] == 'entertainment'][['id', 'lon', 'lat']]
ent.index.name = 'ent'
# set points to the net
for df in [fd, cafe, college, ent]:
    net.set(net.get_node_ids(df['lon'], df['lat'], 250), variable=None, name=df.index.name)<jupyter_output><empty_output><jupyter_text>Load the data we actually care about (observations aka properties) and assign each to a node<jupyter_code>props = gpd.read_file('./data/properties')
props.crs
xs = props.apply(lambda row: row['geometry'].x, axis=1)
ys = props.apply(lambda row: row['geometry'].y, axis=1)
props['node_id'] = net.get_node_ids(xs, ys)
props.head()<jupyter_output><empty_output><jupyter_text>Perform aggregatesWe now have the following ammeneties set on the net:* grocery_store* lrt_stop* tran_cen* fd* cafe* college* entThis is nice because we can quickly and easily change the type of aggregate and the radius of the search<jupyter_code>grocery = net.aggregate(2000, type='count', decay='flat', name='grocery_store')
grocery.name = 'grocery_2k'
light_rail = net.aggregate(2000, type='count', decay='flat', name='lrt_stop')
light_rail.name = 'lightrail_2k'
transit_centers = net.aggregate(2000, type='count', decay='flat', name='tran_cen')
transit_centers.name = 'transit_cen_2k'
food_drink = net.aggregate(2000, type='count', decay='flat', name='fd')
food_drink.name = 'food_drink_2k'
caf = net.aggregate(2000, type='count', decay='flat', name='cafe')
caf.name = 'cafe_2k'
coll = net.aggregate(2000, type='count', decay='flat', name='college')
coll.name = 'college_2k'
entertain = net.aggregate(2000, type='count', decay='flat', name='ent')
entertain.name = 'entertainment_2k'<jupyter_output><empty_output><jupyter_text>Join the aggregates to our observations<jupyter_code>props.shape
sub = props
for agg in [grocery, light_rail, transit_centers, food_drink, caf, coll, entertain]:
    sub = sub.join(agg, on='node_id', how='inner')
sub.head()
# confirm join didn't exclude anyone
sub.shape
# we can drop the geometry
sub.drop('geometry', axis=1, inplace=True)<jupyter_output><empty_output><jupyter_text>ParksParks require a different approach: being polygons they can't accurately be represented with one point, such as a centroid, for our purposes. How many park centroids an observation has within distance x is unecessarily restrictive. In practice what is more important is how many parks could be reached by traveling distance x? Put another way, how many park perimeters are within distance x? So here it seems traditional polygon to polygon intersection is the way to go. It is worth noting however that this will be a 'as the crow flies' measurement instead of a network traversal distance.  Read in the park layer<jupyter_code>parks = gpd.read_file('./data/orca_sites')
parks.head(2)
parks = parks.rename(columns=lambda x: x.lower())
parks['type'].unique()<jupyter_output><empty_output><jupyter_text>Filter down to just parks<jupyter_code>parks = parks[parks['type'] == 'Park and/or Natural Area']<jupyter_output><empty_output><jupyter_text>Check projections<jupyter_code>props.crs
parks.crs<jupyter_output><empty_output><jupyter_text>It generally seems best to use a local projection for these operations<jupyter_code>props = props.to_crs({'init': 'epsg:2913'})
parks = parks.to_crs({'init': 'epsg:2913'})<jupyter_output><empty_output><jupyter_text>Create a 2k buffer for each observation<jupyter_code># the units of this projection are feet: 6561.68 feet == 2k meters
props['geometry'] = props.apply(lambda row: row['geometry'].buffer(6561.68), axis=1)
props.head(2)
parks.head(2)<jupyter_output><empty_output><jupyter_text>Do the spatial join<jupyter_code>parks_sjoin = gpd.sjoin(props, parks[['sitename', 'geometry']], how='inner', op='intersects').drop('index_right', axis=1)
parks_sjoin.head(2)
parks_sjoin.shape<jupyter_output><empty_output><jupyter_text>Group our observations to get back to the correct number of rows<jupyter_code>parks_grp = parks_sjoin[['propertyid', 'sitename']].groupby(by='propertyid').count()
parks_grp.head(2)
# not all properties are within 2k meters of a park so the number of rows here is less than 409
parks_grp.shape
# rename the column
parks_grp.columns = ['parks_2k']
parks_grp.head(2)
sub.head(2)<jupyter_output><empty_output><jupyter_text>Merge parks_grp and sub to create the final dataframe<jupyter_code># the left join will leave some NaNs in the park column so they get filled as zeros
final = pd.merge(sub, parks_grp, how='left', left_on='propertyid', right_index=True).fillna(value=0)
final.head()<jupyter_output><empty_output><jupyter_text>Add tract and blockgroup fips codes so observations can be joined to census data later Read observations, tracts and blockgroups and confirm projections<jupyter_code>props.head(2)
# reread the property file to get the original geometries
props = gpd.read_file('./data/properties')
props.crs
props = props.to_crs({'init': 'epsg:2913'})<jupyter_output><empty_output><jupyter_text>Tracts<jupyter_code># OR
tracts41 = gpd.read_file('./data/tracts_41')
tracts41.crs
tracts41 = tracts41.to_crs({'init': 'epsg:2913'})
# WA
tracts53 = gpd.read_file('./data/tracts_53')
tracts53.crs
tracts53 = tracts53.to_crs({'init': 'epsg:2913'})
tracts = pd.concat([tracts41, tracts53], ignore_index=True)
# pd operations do not always preserve gpd dataframes
type(tracts)
tracts.head(2)
tracts = tracts.rename(columns=lambda x: x.lower())<jupyter_output><empty_output><jupyter_text>Block groups<jupyter_code># OR
bgs41 = gpd.read_file('./data/bgs_41')
bgs41.crs
bgs41 = bgs41.to_crs({'init': 'epsg:2913'})
# WA
bgs53 = gpd.read_file('./data/bgs_53')
bgs53.crs
bgs53 = bgs53.to_crs({'init': 'epsg:2913'})
bgs = pd.concat([bgs41, bgs53], ignore_index=True)
type(bgs)
bgs = bgs.rename(columns=lambda x: x.lower())<jupyter_output><empty_output><jupyter_text>Perform spatial joins to assign each observation to a geography<jupyter_code>tracts.head(2)
bgs.head(2)
census_geos = gpd.sjoin(props, tracts[['geoid', 'geometry']], how='inner', op='within')
census_geos.shape
census_geos.head(2)
census_geos.drop('index_right', axis=1, inplace=True)
census_geos = census_geos.rename(columns={'geoid': 'tract_fips'})
census_geos.head(2)
census_geos = gpd.sjoin(census_geos, bgs[['geoid', 'geometry']], how='inner', op='within')
census_geos.drop('index_right', axis=1, inplace=True)
census_geos = census_geos.rename(columns={'geoid': 'bg_fips'})
census_geos.head(2)<jupyter_output><empty_output><jupyter_text>One final merge of all the data<jupyter_code>master = pd.merge(final, census_geos, how='inner', on='propertyid')
master.shape
master.head(2)
# geometry is no longer needed
master.drop('geometry', axis=1, inplace=True)
master.to_csv('./output.csv', index=False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1373876373321496, 'fraction_numerical': 0.023082898679391853, 'mean_word_length': 6.585858585858586, 'pattern_counts': {'":': 0, '<': 80, '<?xml version=': 0, '>': 80, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '666585', 'license': 'MIT', 'n_tokens_mistral': 3317, 'n_tokens_neox': 3149, 'n_words': 802, 'path': 'spatial_work.ipynb', 'repo_name': 'tadhgfendt/spatial_econometrics', 'size': '9011'}