[postgis-devel] SRID discovery

[maplabs at light42.com]

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 
it.. 

  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... 

  hth
   -Brian

#!/usr/bin/python
## 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'
  sys.exit()

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]
       gdal.PushErrorHandler('CPLQuietErrorHandler')
       err = sr.ImportFromEPSG(int(epsg_code))
       gdal.PopErrorHandler()
       if err == 0:
           gdal.PushErrorHandler('CPLQuietErrorHandler')
           proj4 = sr.ExportToProj4()
           gdal.PopErrorHandler()
           if proj4 == proj4_ref:
               print 'EPSG:',epsg_code
               tFound = True
           #print proj4
       line = f.readline()
   f.close()
   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'