[GIS] Using GDAL GetProjection information to make coordinate conversion in pyproj

coordinate systemgdalpyprojpython

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?

Best Answer

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.

Related Question