Shapefile PRJ to PostGIS SRID Lookup Table Explained

coordinate systempostgissrid

I was wondering if there is such a thing as an shapefile PRJ to PostGIS SRID lookup table? Something that can translate the most standard shapefile PRJ definitions into the likely SRID.

When using PostGIS and pgAdminIII, if you use the postgisgui to import your shapefiles, the SRID is left as "-1". It seems like the tool should be able to parse the Esri PRJ and determine the correct (or at least a couple of options) that are the likely SRID, rather than just leave the default.

Or does the importer have the capability to reproject on the fly if you choose another SRID?

It may seem lazy on my part, but to me it seems curious that this function hasn't already been put in place. Does anyone know if this concept is in the works, or good reason's why it has been left out?

Best Answer

GDAL has a nice convenient interface to the PROJ4 library.

If you are confident with Python, using the GDAL Python bindings, if you import the osr classes you will have very convenient methods for reading and exporting projection representations to a variety of formats like PROJ4, WKT, Esri .PRJ.

For example this script will convert your .PRJ file of your shapefile to WKT and PROJ4 (the last is used from PostGIS):

#! /usr/bin/env python

import sys
from osgeo import osr

def esriprj2standards(shapeprj_path):
    with open(shapeprj_path, 'r') as f:
        prj_txt = f.read()
    srs = osr.SpatialReference()
    srs.ImportFromESRI([prj_txt])
    print('Shape prj is: ' + prj_txt)
    print('WKT is: ' + str(srs.ExportToWkt()))
    print('Proj4 is : ' + str(srs.ExportToProj4()))
    srs.AutoIdentifyEPSG()
    print('EPSG is: ' + str(srs.GetAuthorityCode(None)))

esriprj2standards(sys.argv[1])

Run this on the command line:

$ python esriprj2standards.py /home/pcorti/data/shapefile/country.prj 
Shape prj is: GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]
WKT is: GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]
Proj4 is: +proj=longlat +datum=WGS84 +no_defs 
EPSG is: 4326
Related Question