2010-11-21

Converting WKT projection info to Proj4 (and back)

I've had many occasions where I've had WKT-style projection text (e.g., the contents of the *.prj files that normally accompany ESRI Shapefiles, or the srtext field in the spatial_ref_sys table in PostGIS), and I needed to get that to an equivalent Proj4-compatible projection string (and vice versa).  I've known for some time there had to be an automated way to do this...as the GDAL/OGR libraries have been doing this internally for some time.  I ran into the issue again today, so I decided to poke around a bit, and found the API functions that are available to do this...I wrote my first ever Python scripts that use the GDAL/OGR python bindings (specifically, OSR).  These scripts simply convert WKT projection text to Proj4 (wkt2proj.py) and back (proj2wkt.py).  Maybe there are some better examples elsewhere...but I didn't quite find what I was looking for.  Since my work is mostly in PHP right now, I'd love to see the GDAL/OGR PHP bindings maintained someday, but this will do for now:

================== wkt2proj.py =======================
#!/usr/bin/env python

import os
import sys
import string
import osgeo.osr

if (len(sys.argv) <> 2):
        print 'Usage: wkt2proj.py [WKT Projection Text]'
else:
        srs = osgeo.osr.SpatialReference()
        srs.ImportFromWkt(sys.argv[1])
        print srs.ExportToProj4()
======================================================

================== proj2wkt.py =======================
#!/usr/bin/env python

import os
import sys
import string
import osgeo.osr

if (len(sys.argv) <> 2):
        print 'Usage: proj2wkt.py [Proj4 Projection Text]'
else:
        srs = osgeo.osr.SpatialReference()
        srs.ImportFromProj4(sys.argv[1])
        print srs.ExportToWkt()
======================================================

On Ubuntu, provided you have the python-gdal package installed, the above files executed by python will accept the first argument as a wkt/proj projection string, and return the proj/wkt equivalent.  For proj2wkt.py, you can also take advantage of known EPSG numbers.  For example, the WKT for the Spherical Mercator projection used by most free online mapping services (e.g., Google Maps, OpenStreetMap, etc) will be returned by: proj2wkt.py '+init=epsg:3857'

I haven't tested this on Windows - but as long as the Proj library is installed and the GDAL Python modules are installed and accessible to the Python environment, I imagine they will work the same. 

No comments:

Post a Comment