Skip to content

Instantly share code, notes, and snippets.

@sgillies
Created February 22, 2012 19:27
Show Gist options
  • Save sgillies/1886782 to your computer and use it in GitHub Desktop.
Save sgillies/1886782 to your computer and use it in GitHub Desktop.
Functional style geoprocessing with Fiona, pyproj, Shapely
# Example of using Fiona, pyproj, and Shapely together in a functional
# style.
import functools
import itertools
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 single records.
def transform_coords(func, rec):
# Transform record geometry coordinates using the provided function.
# Returns a record or None.
try:
assert rec['geometry']['type'] == "Polygon"
new_coords = []
for ring in rec['geometry']['coordinates']:
# Here is where we call func() to transform the coords.
x2, y2 = func(*zip(*ring))
new_coords.append(zip(x2, y2))
rec['geometry']['coordinates'] = new_coords
return rec
except Exception, e:
# Writing untransformed features to a different shapefile
# is another option.
logging.exception(
"Error transforming record %s:", rec)
return None
def clean_geom(rec):
# Ensure that record geometries are valid and "clean".
# Returns a record or None.
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)
return rec
except Exception, e:
# Writing uncleanable features to a different shapefile
# is another option.
logging.exception(
"Error cleaning record %s:", rec)
return None
# 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:
# Do the geoprocessing in a functional style. The transform_coords
# function operates on and returns a single record. Its first argument
# is a function that transforms a single coordinate pair between two
# coordinate reference systems. We'll compose partial functions of
# these two into one function.
func = functools.partial(
transform_coords,
functools.partial(
transform,
Proj(source.crs),
Proj(sink.crs) )
)
# Transform the input records, using imap to loop/generate.
results = itertools.imap(func, source)
# The transformed record generator is used as input for a
# "cleaning" generator.
results = itertools.imap(clean_geom, results)
# Remove null records from the results
results = itertools.ifilter(bool, results)
# Finally we write records from the last generator to the output
# "sink" file.
sink.writerecords(results)
# 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)
# The "sink" output file is closed here, no matter what occurs above.
# Use ``with`` like this when you can.
# The "source" input file is closed here, no matter what occurs above.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment