3

I am looking for a Python package that would allow me to input a UTM zone code and output different geo_referencing strings (for example, WKT or Proj4). Are there any lightweight, widely used packages out there that perform this function?

  • GDAL is big but with it you can do this gdal.org/gdalsrsinfo.html. Study the source code and think how much of GDAL and Proj4 you really need. – user30184 Mar 27 '17 at 19:13
2

I'd use the Python bindings for GDAL and OGR. If GDAL/OGR is installed on our machine, the pcs.csv file in gdal-data directory contains an extract of the EPSG Registry about projected coordinate reference systems, so we can filter this file to obtain the result. Suppose we are interested in zone 33, here's a quick and dirty solution:

from osgeo import ogr, osr

in_ds = ogr.GetDriverByName('CSV').Open('pcs.csv')
layer = in_ds.GetLayer()
zone = '33' #'33N'
layer.SetAttributeFilter("COORD_REF_SYS_NAME LIKE '%UTM zone {0}%'".format(zone))
for feature in layer:
    code = feature.GetField("COORD_REF_SYS_CODE")
    name = feature.GetField("COORD_REF_SYS_NAME")
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(int(code))
    print 'EPSG:' + code + ' - ' + name
    print 'Proj4: ', srs.ExportToProj4()
    print 'WKT: ', srs.ExportToWkt()
    print

Note: we can refine the search choosing also the hemisphere (e.g. 33N) changing the value of the zone variable in the code.

This is an extract of the result:

EPSG:2078 - ELD79 / UTM zone 33N
Proj4:  +proj=utm +zone=33 +ellps=intl +towgs84=-115.8543,-99.0583,-152.4616,0,0,0,0 +units=m +no_defs 
WKT:  PROJCS["ELD79 / UTM zone 33N",GEOGCS["ELD79",DATUM["European_Libyan_Datum_1979",SPHEROID["International 1924",6378388,297,AUTHORITY["EPSG","7022"]],TOWGS84[-115.8543,-99.0583,-152.4616,0,0,0,0],AUTHORITY["EPSG","6159"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4159"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","2078"]]

EPSG:2312 - Garoua / UTM zone 33N
Proj4:  +proj=utm +zone=33 +ellps=clrk80 +units=m +no_defs 
WKT:  PROJCS["Garoua / UTM zone 33N",GEOGCS["Garoua",DATUM["Garoua",SPHEROID["Clarke 1880 (RGS)",6378249.145,293.465,AUTHORITY["EPSG","7012"]],AUTHORITY["EPSG","6197"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4197"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","2312"]]

EPSG:2313 - Kousseri / UTM zone 33N
Proj4:  +proj=utm +zone=33 +ellps=clrk80 +units=m +no_defs 
WKT:  PROJCS["Kousseri / UTM zone 33N",GEOGCS["Kousseri",DATUM["Kousseri",SPHEROID["Clarke 1880 (RGS)",6378249.145,293.465,AUTHORITY["EPSG","7012"]],AUTHORITY["EPSG","6198"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4198"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","2313"]]
  • I am having trouble, using your above code I get the error:"'NoneType' 'object has no attribute 'GetLayer' 'getting thrown on: ' layer = in_ds.GetLayer(). Any ideas? – Bryce Frank Mar 30 '17 at 17:11
  • 1
    Look at the 3rd row of the script: the pcs.csv file should be in the same directory where the script is located or alternatively type there its original path. – Antonio Falciano Mar 30 '17 at 17:19
2

GDAL contains the most complete open source implementation and I'm not aware of any port to Python.

Rasterio does the same kind of thing as GDAL's Python bindings and calls the same C library functions.

>>> from rasterio.crs import CRS
>>> CRS.from_epsg(4326).wkt
'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]

Rasterio is not lightweight by any measure.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.