Shapefile PRJ to PostGIS SRID lookup table?

srsly. I want one too.

Many people seem to look them up at http://spatialreference.org

When you import shapefiles using PostGIS (and the PostGIS loader for PGAdmin), it looks up the proj information in a table called spatial_ref_sys.

From what I understand, the standard spatial_ref_sys table packaged with PostGIS includes only OGC WKT (Open Geospatial Consortium Well Known Text) representations of some Spatial Reference Systems and NOT the ESRI Spatial reference systems.

From the PostGIS 1.5.2 documentation: >

The spatial_ref_sys table is a PostGIS included and OGC compliant database table that lists over 3001 known spatial reference systems and details needed to transform/reproject between them.

Although the PostGIS spatial_ref_sys table contains over 3000 of the more commonly used spatial reference system definitions that can be handled by the proj library, it does not contain all known to man and you can even define your own custom projection if you are familiar with proj4 constructs. Keep in mind that most spatial reference systems are regional and have no meaning when used outside of the bounds they were intended for.

An excellent resource for finding spatial reference systems not defined in the core set is http://spatialreference.org/ Some of the more commonly used spatial reference systems are: 4326 - WGS 84 Long Lat, 4269 - NAD 83 Long Lat, 3395 - WGS 84 World Mercator, 2163 - US National Atlas Equal Area, Spatial reference systems for each NAD 83, WGS 84 UTM zone - UTM zones are one of the most ideal for measurement, but only cover 6-degree regions.

Various US state plane spatial reference systems (meter or feet based) - usually one or 2 exists per US state. Most of the meter ones are in the core set, but many of the feet based ones or ESRI created ones you will need to pull from spatialreference.org.

However ogr2ogr contains ESRI spatial ref systems as a I recently learned through the generosity of others.

In both ogr2ogr and spatial_ref_sys, it seems that text contained in the .proj file is compared with a table of OGC WKT, which is a slightly different text format from the ESRI WKT format that you often find in a .proj file. Also, I'm not sure how PostGIS looks up each SRS, but the small differences between ESRI WKT and OGC WKT might result in failed matches.

It seems like it would be simple to attach ESRI spatial ref systems to the default spatial_ref_sys table in PostGIS. Maybe someone already has, with some patch or a script.

I could be wrong, because I've just been running into this for the past few days, and I've been frustrated with the same thing. Maybe someone else knows a great resource?

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

It's been a while since I used POSTGIS srids but if they are just EPSG codes then you can use http://prj2epsg.org/search to look them up from (broken) ESRI.prj files.

It's been a while since I needed to, but as I recall, http://spatialreference.org/ in addition to allowing you to search, also gives you the option of uploading a prj file.

Then it will as one of the output options give you the equivalent postgis insert to insert into spatial_ref_sys table.

For the insert statement it gives, I replace the generated srid it creates with the EPSG or ESRI one. If you get a primary key violation, then you know most likely its already in the table.

As a mix of solutions I have created a script to help me load arbitrary shapefiles into postgis. It also tries to detect the encoding of the DBF.

from chardet.universaldetector import UniversalDetector
import os.path
import sys
import dbfUtils
import sys
from osgeo import osr
from urllib import urlencode
from urllib2 import urlopen
import json

shp_file = sys.argv[1]
dbf_file = shp_file[0:-4] + '.dbf'
prj_file = shp_file[0:-4] + '.prj'

#Try detecting the SRID, by default we set to 4326 and hope the best
srid=4326
if os.path.isfile(prj_file):
    prj_filef = open(prj_file, 'r')
    prj_txt = prj_filef.read()
    prj_filef.close()
    srs = osr.SpatialReference()
    srs.ImportFromESRI([prj_txt])
    srs.AutoIdentifyEPSG()
    code = srs.GetAuthorityCode(None)
    if code:
        srid= code
    else:
        #Ok, no luck, lets try with the OpenGeo service
        query = urlencode({
            'exact' : True,
            'error' : True,
            'mode' : 'wkt',
            'terms' : prj_txt})
        webres = urlopen('http://prj2epsg.org/search.json', query)
        jres = json.loads(webres.read())
        if jres['codes']:
            srid = int(jres['codes'][0]['code'])

#Try to detect the encoding
dbf = open(dbf_file, 'rb')
db = dbfUtils.dbfreader(dbf)

detector = UniversalDetector()
for row in db:
    detector.feed(str(row))
    if detector.done: break
detector.close()
dbf.close()

encoding = detector.result["encoding"]
if encoding=="ascii":
    encoding="LATIN1"

print "shp2pgsql -s %s -k -i -I -W %s %s.shp public.importing_table" %(srid,encoding,shp_file)