Functional style geoprocessing with Fiona, pyproj, Shapely

even-more-functional.py

# 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)

with-pyproj-shapely-functional.py

# 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.