Dealing with custom projections (from prj file) in PostGIS (also GeoServer and OpenLayers)

Question

Basically I want to perform spatial queries like distance calculations in PostGIS using data from shapefiles with different projections, some of which have no SRID code that I can find and also which the python osr package can't seem to handle. What do I do? I have little understanding of projections and how to deal with them.

Full story:

I'm creating a system using OpenLayers, GeoServer, and PostGIS. My data is from shapefiles that I loaded into PostGIS using shp2psql. Most of the shapefiles didn't have a .prj file associated with them and seem to just be long/lat / EPSG:4326. However for some of the files the coordinates have numbers in the millions which are clearly not degrees (fortunately there were .prj files here). In PostGIS the data appears as a geometry column, not geography. Until now I've been mostly concerned with displaying the data in OpenLayers, so I took care of projection in GeoServer by specifying an EPSG code when I added a new layer. Usually this was 4326, once it was a standard code I found using http://www.prj2epsg.org/search, and once I created a custom projection by adding the WKT to the epsg.properties files in the data directory. Everything seems to display fine in OpenLayers then; it's hard to be sure since I don't know the data personally but it looks like it all lines up. I haven't had to write any code on the client side (OpenLayers) relating to projections, everything is in terms of longitude and latitude and very easy.

Now, suppose I specify in the client (e.g. by clicking on the OpenLayers map) the coordinates of a point using longitude and latitude, and a distance such as 100 km, and I want to display all the features from a certain layer that are within said distance from said point. So, in PostGIS I'd like to say something like:

SELECT gid, geom
FROM layer
WHERE ST_Distance(geom, ST_GeomFromText('POINT(24 -30)')) < 100000;

Right now this returns everything (for most layers) since the geom column is a geometry type instead of geography. As far as I can see, I need to turn everything into a geography type, transforming into the EPSG:4326 where necessary. Then I should just be able to declare everything in GeoServer as EPSG:4326 and no longer worry about anything there, and ST_Distance should return results in metres calculated on the spheroid.

But for some of my data I can't see how to do this. For example, here is the text from the .prj file for one of the layers:

PROJCS["SA Albers Equal Area",GEOGCS["GCS_Hartebeesthoek_1994",DATUM["D_Hartebeesthoek_1994",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Albers"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",24.0],PARAMETER["Standard_Parallel_1",-24.0],PARAMETER["Standard_Parallel_2",-33.0],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

This is a projection for South Africa, and none of the results from http://www.prj2epsg.org/search look relevant. So I want to define it as a custom projection in PostGIS. As far as I understand I can't use WKT for this and need the parameters to be expressed in proj4. I tried running it through this python script seen around here:

#!/usr/bin/env python

import osr
import sys

def main(prj_file):
    prj_text = open(prj_file, 'r').read()
    srs = osr.SpatialReference()
    if srs.ImportFromWkt(prj_text):
        raise ValueError("Error importing PRJ information from: %s" % prj_file)
    print srs.ExportToProj4()
    #print srs.ExportToWkt()

if __name__=="__main__":
    main(sys.argv[1])

and I got the following:

ERROR 6: No translation for Albers to PROJ.4 format is known.

Now I don't know how to proceed. I'd also like confirmation that my reasoning about everything else I've mentioned is sensible.

Answer 1

1) You seem to have a mix of projections in your original files (shapes) so, how are you loading the files to Postgis? Specifically, what projection are you passing as a parameter to shp2pgsql? Because if I remember correctly, if you don't specify a srid then Postgis assumes it to be in local coordinates (srid 0) and planar units so that maybe why you're seeing your layers in the geometry columns (projected units).

2) ST_GeomFromText without an srid argument (the variant you are using) returns a geometry with no spatial reference so you can't expect to have meaningful distance calculations using that variant. It would be better to use ST_GeomFromText(geom,srid) and then operate on a specific projection.

3)I'ts not really hard to update the spatial reference table in postgis and the list of supported reference systems in proj.4 (the one gdal uses) to include the projection you are using: For postgis: http://postgis.refractions.net/documentation/manual-1.4/ch04.html#spatial_ref_sys For gdal: http://www.gdal.org/osr_tutorial.html

Summarizing, you need to be sure your layers have (correct) spatial reference in Postgis and then use this spatial reference to calculate distance and create points from text. It wouldn't hurt to do a little reading in projections and transformations.

*edit fixed wrong link

Answer 2

Information on creating Projections is available here from the Esri Site. See the part on Manual Creation down the page.

There is an excellent discussion on ShapeFile Projections and PostGis here and an example of the SQL to insert the Custom Projection into PostGis can be seen here

The SQL to insert the Custom Projection is probably of the most interest to you ..

Answer 3

There is a bug in gdal. You have to chgange "Albers" to "Albers_Conic_Equal_Area" in your WKT then you can use whatever tool you are using to create the proj4.