# Geometry Projection Speed Test

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?

## 1 Answer

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`.
``````
• +1. Also Shapely Points, LinearRings and LineStrings have a numpy array interface, so you can do something like `projected_coords = numpy.vstack(pyproj.transform(origin, target, *numpy.array(geom).T)).T` Jul 16, 2013 at 1:23
• This is awesome @sgillies. For some reason my version of shapely doesn't have transform? shapely.__version__: '1.2.17'. I'll try grabbing the source directly. Jul 16, 2013 at 16:41
• Oops, sorry. Coming in version 1.2.18 (this weekend, probably). Jul 19, 2013 at 16:03