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