Geometry Projection Speed Test

Question

Lately I've been using the OGR projection classes that come with ogr/gdal, but pyproj was recommended to me, so I thought I'd give it a try. To help me decide if I should make the switch, I did a speed test. The following is a small (almost) reproducible example that I came up with to test the two. I'm not sure if this test is totally fair, so comments and recommendations are welcome!

Imports first, to make sure we start off with a level playing field:

from pandas import Series # This is what our geometries are stored in
from shapely import wkb
import functools
from shapely.geometry import shape
from osgeo import ogr
# The following two imports are the important ones
from pyproj import Proj, transform
from osgeo.osr import SpatialReference, CoordinateTransformation

Because I'm storing the shapely geometries in a pandas 'Series', the functions need to work with Series.apply(). Here I define two functions (one using 'ogr.osr' and one using 'pyproj') to perform coordinate transformations inside a call to Series.apply():

def transform_osr(geoms, origin, target):
    target_crs = SpatialReference()
    target_crs.ImportFromEPSG(origin)
    origin_crs = SpatialReference()
    origin_crs.ImportFromEPSG(origin)
    transformer = CoordinateTransformation(origin_crs, target_crs)
    def trans(geom):
        g = ogr.CreateGeometryFromWkb(geom.wkb)
        if g.Transform(transformer):
            raise Exception("Ahhhh!")
        g = wkb.loads(g.ExportToWkb())
        return g
    return geoms.apply(trans)

def transform_pyproj(geoms, origin, target):
    target = Proj(init="epsg:%s" % target)
    origin = Proj(init="epsg:%s" % origin)
    transformer = functools.partial(transform, origin, target)
    def trans(geom):
        interface = geom.__geo_interface__
        geom_type = interface['type']
        coords = interface['coordinates']
        result = apply_to_coord_pairs(transformer, coords)
        return shape({'coordinates':result, 'type':geom_type})
    def apply_to_coord_pairs(fun, geom):
        return [not all([hasattr(y, "__iter__") for y in x]) and \
                fun(*x) or apply_to_coord_pairs(fun, x) for x in geom]
    return geoms.apply(trans)

Each of these functions takes an EPSG code as input for the origin and destination coordinate reference systems. Both libraries offer alternative ways to define projection information, but EPSG codes keeps the code pretty simple for now.

Here are the results, using the %timeit magic function in ipython:

In [1]: %timeit transform_pyproj(l, 29903, 4326)
1 loops, best of 3: 641 ms per loop

In [2]: %timeit transform_osr(l, 29903, 4326)
10 loops, best of 3: 27.4 ms per loop

Clearly the 'ogr.osr' version is faster, but is it a fair comparison? The 'pyproj' version is done on individual points, and is mostly run in Python, whereas the 'ogr.osr' version operates directly on the OGR geometry object. Is there a better way to compare these?

Answer 1

Pyproj is a Python C extension which uses the PROJ4 library and osgeo.ogr is a Python C extension which uses the PROJ4 library. If you're only considering coordinate projection, in the fairest test they'd be almost equal.

Pyproj's transform can operate on arrays of coordinate values, so you only need to call it once per line or ring instead of for every pair. This should speed things up quite a bit. Example: https://gist.github.com/sgillies/3642564#file-2-py-L10.

(Update) Also, Shapely provides a function that transforms geometries in 1.2.16:

Help on function transform in module shapely.ops:

transform(func, geom)
    Applies `func` to all coordinates of `geom` and returns a new
    geometry of the same type from the transformed coordinates.

    `func` maps x, y, and optionally z to output xp, yp, zp. The input
    parameters may iterable types like lists or arrays or single values.
    The output shall be of the same type. Scalars in, scalars out.
    Lists in, lists out.

    For example, here is an identity function applicable to both types
    of input.

      def id_func(x, y, z=None):
          return tuple(filter(None, [x, y, z]))

      g2 = transform(id_func, g1)

    A partially applied transform function from pyproj satisfies the
    requirements for `func`.

      from functools import partial
      import pyproj

      project = partial(
          pyproj.transform,
          pyproj.Proj(init='espg:4326'),
          pyproj.Proj(init='epsg:26913'))

      g2 = transform(project, g1)

    Lambda expressions such as the one in

      g2 = transform(lambda x, y, z=None: (x+1.0, y+1.0), g1)

    also satisfy the requirements for `func`.