Created
February 22, 2012 19:27
-
-
Save sgillies/1886782 to your computer and use it in GitHub Desktop.
Functional style geoprocessing with Fiona, pyproj, Shapely
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# Example of using Fiona, pyproj, and Shapely together in a functional | |
# style. | |
import functools | |
import logging | |
import sys | |
import fiona | |
from pyproj import Proj, transform | |
from shapely.geometry import mapping, shape | |
logging.basicConfig(stream=sys.stderr, level=logging.INFO) | |
# The two functions below operate on sequences of records and are generators | |
# for new sequences. | |
def transform_coords(func, records): | |
# Transform record geometry coordinates using the provided function. | |
for rec in records: | |
try: | |
assert rec['geometry']['type'] == "Polygon" | |
new_coords = [] | |
for ring in rec['geometry']['coordinates']: | |
x2, y2 = func(*zip(*ring)) | |
new_coords.append(zip(x2, y2)) | |
rec['geometry']['coordinates'] = new_coords | |
except Exception, e: | |
# Writing untransformed features to a different shapefile | |
# is another option. | |
logging.exception( | |
"Error transforming record %s:", rec['id']) | |
yield rec | |
def clean_geom(records): | |
# Ensure that record geometries are valid and "clean". | |
for rec in records: | |
try: | |
geom = shape(rec['geometry']) | |
if not geom.is_valid: | |
clean = geom.buffer(0.0) | |
assert clean.is_valid | |
assert clean.geom_type == 'Polygon' | |
geom = clean | |
rec['geometry'] = mapping(geom) | |
except Exception, e: | |
# Writing uncleanable features to a different shapefile | |
# is another option. | |
logging.exception( | |
"Error cleaning record %s:", rec['id']) | |
yield rec | |
# I'm experimenting with changing the name of Fiona's file opening function | |
# from collection() to open(). Feedback greatly appreciated. | |
fiona.open = fiona.collection | |
# Here is where the work is done. | |
with fiona.open("docs/data/test_uk.shp", "r") as source: | |
with fiona.open( | |
"with-pyproj-shapely2.shp", | |
"w", | |
driver=source.driver, | |
schema=source.schema, | |
crs={'init': "epsg:27700", 'no_defs': True} | |
) as sink: | |
# The transform_coords function returns a generator for records. | |
# Its first argument is a function that transforms a single | |
# coordinate pair between two coordinate reference systems. | |
results = transform_coords( | |
functools.partial(transform, Proj(source.crs), Proj(sink.crs)), | |
source ) | |
# The transformed record generator is used as input for the | |
# "cleaning" generator function. | |
results = clean_geom(results) | |
# Finally we write records from the last generator to the | |
# output "sink" file. | |
sink.writerecords(results) | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment