Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Download City of Austin permit data and use shapely to append neighborhood name and census tract attributes.
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
#!/usr/bin/env python3
'''
Download City of Austin permit data and append neighborhood name and
census tract attributes.
'''
import csv
import json
import pdb
import requests
from shapely.geometry import shape, point
def get_data(base_url, limit=1000):
'''
Get data from City of Austin Open Data Portal via Socrata API
as described here: https://dev.socrata.com/docs/
'''
done = False
offset = 0
data = []
while not done:
# request records until fewer records are returned than limit
url = base_url + ' limit {} offset {}'.format(limit, offset*limit)
res = requests.get(url)
res.raise_for_status()
records = res.json()
data = data + records
offset += 1
if len(records) < limit:
done = True
return data
def read_csv(f):
'''
Read CSV into memory as list of dicts
'''
with open(f, 'r') as fin:
reader = csv.DictReader(fin)
return [row for row in reader]
def read_json(f):
'''
Load (geo)JSON into memory
'''
with open(f, 'r') as fin:
return json.loads(fin.read())
def create_points(data, x_field, y_field, geom_key='geometry'):
'''
Create shapely geometry from list of dicts
'''
for row in data:
x = row.get(x_field)
y = row.get(y_field)
if x and y:
row[geom_key] = point.Point(float(x), float(y))
else:
row[geom_key] = None
return data
def point_in_poly(points=None, polys=None, key_map=[], geom_key='geometry'):
'''
Get properties of polygon that contains input point. Assumes polygons
do not overlap.
points: list of dicts with shapely geometries
polys: geojson polygon feature collection
key_map: list of dicts which identifies poly properties to extract. Each
dict must contain 'source' and 'dest' fieldname keys
Returns points with polygon properties specified in key_map.
Assigns property=None to points that are missing geometries or do not
intersect with polygon
'''
for pt in points:
if pt[geom_key]:
contained = False
p = point.Point(pt[geom_key])
for feature in polys['features']:
polygon = shape(feature[geom_key])
if polygon.contains(p):
properties = {key['dest'] : feature['properties'][key['source']] for key in key_map}
pt.update(properties)
contained = True
continue
if not contained:
# add empty keys to rows not contained by poly
properties = { key['dest'] : None for key in key_map }
pt.update(properties)
else:
# add empty keys to rows missing geometry
properties = { key['dest'] : None for key in key_map }
pt.update(properties)
return points
def write_csv(data, fname):
fieldnames = data[0].keys()
with open(fname, 'w', newline='\n') as fname:
writer = csv.DictWriter(fname, fieldnames=fieldnames)
writer.writeheader()
for row in data:
writer.writerow(row)
# map fields from containing polygon to input point
key_map = {
'hoods' : [
{
'source' : 'Name',
'dest' : 'neighborhood'
}
],
'tracts' : [
{
'source' : 'GEOID',
'dest' : 'GEOID'
}
]
}
# "City of Austin Issued Construction Permits
resource = 'x9yh-78fz'
# Query to select residential demolition permits for single and two-family homes
query = '''SELECT permit_number,work_class,calendar_year_issued,issue_date,description,permit_location,status_current,latitude,longitude,permit_class_mapped WHERE upper(permit_class) LIKE '%25645%25' OR upper(permit_class) LIKE '%25646%25' ORDER BY permit_number DESC'''
# API Endpoint
url = 'https://data.austintexas.gov/resource/{}.json?$query={}'.format(resource, query)
# Get permit data from city of austin open data portal
print('Get data')
permits = get_data(url)
# Create shapely geometries from x/y attributes
permits = create_points(permits, 'longitude', 'latitude')
# Read zillow neighborhood data into memory (converted to geojson with Mapshaper)
neighborhoods = read_json('austin_neighborhoods.json')
# Read austin census tracts data into memory (clipped/converted to geojson with mapshaper)
tracts = read_json('austin_census_tracts.json')
# Append neighborhood name to permit records
print('Get neighborhoood name')
permits = point_in_poly(
points=permits,
polys=neighborhoods,
key_map=key_map['hoods']
)
# Append census tract geoid to permit records
print('Get census tract')
permits = point_in_poly(
points=permits,
polys=tracts,
key_map=key_map['tracts']
)
write_csv(permits, 'demo_permits.csv')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment