1

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.

  • You have to know the projection your shape files are in. No way around it, at either import time or specify on the column after import. – Mark Cupitt Jul 27 '14 at 3:14
  • I know the projections, my issue is that PostGIS doesn't support all of them and I can't see how to add it. – alexmojaki Jul 27 '14 at 18:18
2

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

  • It appears you've posted the same link twice. – Chris W Jul 26 '14 at 18:00
  • Thanks. 1 and 2 were points I was not aware of and they helped me solve my problem. I'm not sure what you were saying in 3 regarding supporting a new reference system in proj4, but I found that I could extract the proj4 text of a layer using QGIS. – alexmojaki Sep 11 '14 at 2:29
  • @alexmojaki The thing is you can use the proj4 text of a spatial reference system to update the spatial_ref_sys table in postgis with an insert statement (like this, for example: insert statement). In this way postgis can handle the new projection even if it was not originally supprted. Something similar applies for proj4 and gdal – plablo09 Sep 12 '14 at 14:42
  • @plablo09 Yes, that's what I did. My problem was that I only had the WKT and didn't know how to get the proj4 text from it. That's essentially what my question was about. The python script was meant to do that but apparently it can't handle everything. – alexmojaki Sep 12 '14 at 19:01
0

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 ..

0

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.

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.