Using GDAL GetProjection information to make coordinate conversion in pyproj

Question

I'm opening a geoTIFF file in Python using gdal. Upon reading with the GetProjection() method, I find the following information

PROJCS["NAD83_HARN_UTM_zone_15N",GEOGCS["GCS_North_American_1983_HARN",DATUM["D_North_American_1983_HARN",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["false_easting",500000.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-93.0],PARAMETER["scale_factor",0.9996],PARAMETER["latitude_of_origin",0.0],UNIT["Meter",1.0],VERTCS["NAVD_1988",VDATUM["North_American_Vertical_Datum_1988"],PARAMETER["Vertical_Shift",0.0],PARAMETER["Direction",1.0],UNIT["Meter",1.0]]]

Using this information, and my x,y information for my points, I'd like to convert them to WGS84 to compare with Google maps. So I tried doing the following:

wgs84=pyproj.Proj(init='epsg:4326') p2=pyproj.Proj(init='epsg:3745') print pyproj.transform(p2,wgs84,*upper_left) print pyproj.transform(p2,wgs84,*bottom_right)

where I manually looked up the code for NAD83_HARN_UTM_zone_15N.

Surely me looking up the code for the p2 projection is not the best way to do this, as different files can have different projections listed, and I can't process things automatically if I have to do this step manually.

How do I do this without manual lookup?

Answer 1

You can use the osr module (part of the standard GDAL install, so it should come with your Python bindings) and do something like this:

import osr
import gdal
inDS = gdal.open(r'c:\somedirectory\myRaster.tif')
inSRS_wkt = inDS.GetProjection()  # gives SRS in WKT
inSRS_converter = osr.SpatialReference()  # makes an empty spatial ref object
inSRS_converter.ImportFromWkt(inSRS)  # populates the spatial ref object with our WKT SRS
inSRS_forPyProj = inSRS_converter.ExportToProj4()  # Exports an SRS ref as a Proj4 string usable by PyProj

Of course, if your ultimate goal is to reproject your raster, I would suggest calling the gdalwarp utility from within your code as this will be a lot easier.