[postgis-devel] SRID discovery

maplabs at light42.com maplabs at light42.com
Thu Oct 6 11:15:06 PDT 2011

Hi All-

  I was reading the notes from the PostGIS code sprint in Denver, and 
realized there was something said there about discovering an SRID .. I 
know in my limited experience as a neo-geo, that having the PROJ string 
for a raster or SHP file is not the same as knowing an EPSG code for 

  I asked on IRC (of course) long ago, and the GDAL wizard EvenR wrote 
the following code in response in literally minutes.. 
I have found it useful and modified it a bit.. it uses the python GDAL 
bindings and files supplied with a standard GDAL install.. 
Pass a PROJ string and try to discover an EPSG code for it... 


## original by EvenR  18May10
## rev dbb  27Nov10

import sys
import osr
import gdal

inFileName = sys.argv[1]
tG = gdal.Open(inFileName)

wkt = tG.GetProjectionRef()
sr = osr.SpatialReference()
err = sr.ImportFromWkt(wkt)
if err != 0:
  print 'EPSG:Null'

proj4_ref = sr.ExportToProj4()

def do_proj4_cmp( inFileName ):
   f = open( inFileName, 'r')
   # skip header
   line = f.readline()
   tFound = False
   line = f.readline()
   while line != '':
       epsg_code = (line.split(','))[0]
       err = sr.ImportFromEPSG(int(epsg_code))
       if err == 0:
           proj4 = sr.ExportToProj4()
           if proj4 == proj4_ref:
               print 'EPSG:',epsg_code
               tFound = True
           #print proj4
       line = f.readline()
   return tFound


res1 = do_proj4_cmp('/usr/local/share/gdal/pcs.csv')
res2 = do_proj4_cmp('/usr/local/share/gdal/gcs.csv')

if not (res1 or res2):
  print 'EPSG:Null'

More information about the postgis-devel mailing list