referencing.transform

This is a proposed module for transforming coordinates between reference systems.

I have a customer that needs the ability to transform points between spatial reference systems, but does not want/need the entirety of mapscript, ogr, or PCL. To this end, I'm going to write a new PROJ.4 based point transformer class that will transform a list of coordinate tuples.

Am I reinventing the wheel? Isn't this already done by Thuban's PyProjection?? Not quite. PyProjection? appears to only do the forward and inverse transforms between lat/long and projected systems. My goal is the full monty, the same capabilities as PROJ.4's cs2cs utility.

My intent is that the referencing module should be able to stand alone and be used without requiring the rendering module of PCL. I'm pretty sure that my original implementation will borrow heavily from MapServer's mapproject.c code.

implementation

In fact, that's basically how I did it. One C extension module, referencing.transform._proj4 which uses functions from the proj API. This extension module knows nothing about other PCL classes, just Python strings, lists, and tuples. A referencing.transform.proj4 module provides the main API:

>>> from cartography.referencing.srs import SpatialReference
>>> from cartography.referencing.transform.proj4 import ProjTransform
>>> utm13 = SpatialReference(epsg=26913)
>>> wgs84 = SpatialReference(epsg=4326)
>>> t = ProjTransform(utm13, wgs84)
>>> t.transform([(-105, 39), (-106, 40)])
>>> [(499999.99999999889, 4316776.5830979366), (414639.53815663763, 4428236.0645195534)]

usage

As shown above the proj4 module has a ProjTransform? class, initialized from two instances of SpatialReference?. One is the source SRS, the other is the destination SRS.

t = ProjTransform(destination_srs, source_srs)

This class has a transform() method which takes a list of X,Y coordinate tuples (see above) and returns another such list

transformed_coords = t.transform(source_coords)

The order of the parameters goes a bit against the grain if you're accustomed to mapscript and OGR, but I think it makes more sense. Chain it all out like

transformed_coords = ProjTransform(destination_srs, source_srs).transform(source_coords)

and you see that on the right we have all the sources, and on the left all the outputs. I see it as a built-in hint.

installation

The only external dependency of the referencing module is the PROJ.4 library, available from http://proj.maptools.org/. To build and install the module

$ cd PCL/referencing
$ sudo python setup.py install

testing

There are unit tests and coverage measures

[sean@lenny referencing]$ ./runcoverage 
running coverage tests ...
...............
----------------------------------------------------------------------
Ran 15 tests in 0.068s

OK
Name    Stmts   Exec  Cover
---------------------------
srs        19     18    94%
proj4       7      7   100%
---------------------------
TOTAL      26     25    96%

and the tests run leak-free under valgrind

[sean@lenny referencing]$ ./rungrind 
==18242== Memcheck, a memory error detector for x86-linux.
==18242== Copyright (C) 2002-2004, and GNU GPL'd, by Julian Seward et al.
==18242== Using valgrind-2.2.0, a program supervision framework for x86-linux.
==18242== Copyright (C) 2000-2004, and GNU GPL'd, by Julian Seward et al.
==18242== For more details, rerun with: -v
==18242== 
...............
----------------------------------------------------------------------
Ran 15 tests in 3.586s

OK
==18242== 
==18242== ERROR SUMMARY: 0 errors from 0 contexts (suppressed: 1208 from 7)
==18242== malloc/free: in use at exit: 681366 bytes in 119 blocks.
==18242== malloc/free: 2822 allocs, 2703 frees, 1599734 bytes allocated.
==18242== For counts of detected errors, rerun with: -v
==18242== searching for pointers to 119 not-freed blocks.
==18242== checked 4106828 bytes.
==18242== 
==18242== LEAK SUMMARY:
==18242==    definitely lost: 0 bytes in 0 blocks.
==18242==    possibly lost:   0 bytes in 0 blocks.
==18242==    still reachable: 681166 bytes in 118 blocks.
==18242==         suppressed: 200 bytes in 1 blocks.
==18242== Reachable blocks (those to which a pointer was found) are not shown.
==18242== To see them, rerun with: --show-reachable=yes