38

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?

| improve this question | |
9

Borrowing the idea from @iant, here is a PL/Python3 module that will look up the EPSG SRID integer codes from a PRJ file using the http://prj2epsg.org web service.

First, install PL/Python3:

CREATE LANGUAGE plpython3u;

an now add the SQL function, which has code written for Python 3:

CREATE OR REPLACE FUNCTION prj2epsg(prj_file text) RETURNS integer AS
$BODY$

import json
from urllib.parse import urlencode
from urllib.request import urlopen

with open(prj_file, 'r') as fp:
    prj_txt = fp.read()

query = urlencode({
    'exact': True,
    'error': True,
    'mode': 'wkt',
    'terms': prj_txt})

webres = urlopen('http://prj2epsg.org/search.json', query.encode())
jres = json.loads(webres.read().decode())

return int(jres['codes'][0]['code'])

$BODY$ LANGUAGE plpython3u VOLATILE COST 100;

To use it from PostgreSQL:

SELECT prj2epsg(E'C:\\Temp\\countries.prj');

returns 4326 for my test Shapefile.

| improve this answer | |
  • I am going to mark this as the solution. While the others are all excellent, I love this idea. Now if we can just get someone with the coding ability to include this type of functionality into the pgAdmin PostGIS shapefile loader so that it automatically determines the correct SRID upon reading the SHP. I will keep my fingers crossed. – RyanKDalton Mar 28 '11 at 23:07
  • 2
    The caveats of course is it requires an Internet connection and depends on an external web service that needs to up and running. – Mike T Mar 28 '11 at 23:33
58

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):
   prj_file = open(shapeprj_path, 'r')
   prj_txt = prj_file.read()
   srs = osr.SpatialReference()
   srs.ImportFromESRI([prj_txt])
   print 'Shape prj is: %s' % prj_txt
   print 'WKT is: %s' % srs.ExportToWkt()
   print 'Proj4 is: %s' % srs.ExportToProj4()
   srs.AutoIdentifyEPSG()
   print 'EPSG is: %s' % 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
| improve this answer | |
  • I came across two issues this approach: (1) +proj=longlat +datum=WGS84 +no_defs isn't in the spatial_ref_sys table, so you can't use the output to look up the SRID; and (2) I can't seem to find any SRID property or method (there's a handy ImportFromEPSG(SRID) method, but not the other way around) – Mike T Mar 27 '11 at 21:44
  • 4
    I have updated the script with a call to the AutoIdentifyEPSG() method that will do the magic ;) – capooti Apr 6 '11 at 14:41
  • Very cool. Excellent work! – RyanKDalton Apr 6 '11 at 15:46
  • When gdalsrsinfo and ogrinfo fail you, this is the way to go! – kontextify Jun 4 '17 at 12:07
  • 1
    Notice that srs.GetAuthorityCode(None) might still return none if no close SRID was identified. – astrojuanlu Jul 26 '17 at 16:10
20

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.

| improve this answer | |
  • This is a really smart website. Looking at the API, you could write a nice server-side script to automate the procedure. – Mike T Mar 27 '11 at 9:27
4

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)
| improve this answer | |
3

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?

| improve this answer | |
1

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.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.