Skip to content

Instantly share code, notes, and snippets.

@jaklinger
Created January 5, 2018 15:55
Show Gist options
  • Save jaklinger/6bf635e78ce323a7923fc424ea94bc15 to your computer and use it in GitHub Desktop.
Save jaklinger/6bf635e78ce323a7923fc424ea94bc15 to your computer and use it in GitHub Desktop.
Example of dissolving lat/lon points in different projection systems, using UK TTWA shape files
# Geometry, shapefiles and projections
import fiona
from shapely.geometry import shape
from shapely.geometry import Point
import pyproj
# Generate a function to create a UK East/North point from Lon/Lat
wgs84 = pyproj.Proj(init = 'epsg:4326')
ukgrid = pyproj.Proj(init = 'epsg:27700')
EnPoint = lambda lon, lat : Point(*pyproj.transform(wgs84, ukgrid, lon, lat))
# File containing lat lon points
df = pd.read_json('gtr_academic_private.json')
# Open TTWA shapes
shp_filename = 'Travel_to_Work_Areas_December_2011_Full_Clipped_Boundaries_in_United_Kingdom.shp'
ttwa_counts = {}
with fiona.open(shp_filename) as ttwas:
# Iterate through ttwas
for ttwa in ttwas:
id_ = ttwa['id'] # TTWA id
if id_ not in ttwa_counts:
ttwa_counts[id_] = 0
# Iterate through table looking for lat/lon points in TTWA
for _,row in df.iterrows():
lat = row['lat']
lon = row['lon']
try:
point = EnPoint(lon,lat)
except:
continue
if not point.within(shape(ttwa['geometry'])):
continue
ttwa_counts[id_] += 1
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment