Projections and transformation (2)

In a comment on my previous post about projections and transformations Gene pointed out that I could use Ogr with it its Python bindings. I already used the C# bindings from Ogr in the past but I thought it would be a nice exercise to use Ogr from Python.

If you are on windows and need to install gdal/ogr and its Python bindings I suggest you to use the following installation guideline.

Below is the short code I came up with but there are numerous other samples on the web sample 1 sample 2 and off course the documentation on the gdal site. Attention I got an error when trying to use the TransformPoints method.

from osgeo import osr

def get_spatialref(epsg_code):
    spatialref = osr.SpatialReference()
    spatialref.ImportFromEPSG(epsg_code)
    return spatialref

lambert72_spatialref = get_spatialref(31300)
etrs89_spatialref = get_spatialref(4258)

coordinate_transformation = osr.CoordinateTransformation(lambert72_spatialref, etrs89_spatialref)

## start location : random point in Belgium
x = 155335.20410
y = 166096.08270
print coordinate_transformation.TransformPoint(x, y)

Another solution is to install and use the Python bindings for the Proj4 library.
The following sample code comes from its documentation.

import pyproj

# projection 1: UTM zone 15, grs80 ellipse, NAD83 datum
# (defined by epsg code 26915)
p1 = pyproj.Proj(init='epsg:26915')
p1_2 = get_spatialref(26915)
# projection 2: UTM zone 15, clrk66 ellipse, NAD27 datum
p2 = pyproj.Proj(init='epsg:26715')
p2_2 = get_spatialref(26715)
# find x,y of Jefferson City, MO.
x1, y1 = p1(-92.199881,38.56694)
# transform this point to projection 2 coordinates.
x2, y2 = pyproj.transform(p1,p2,x1,y1)
print x2, y2

If you need to look up a spatial reference you might wanna take a look at spatialreference.org.

Succes !