[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
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'
More information about the postgis-devel
mailing list