0

I'm looking for the correct method using GDAL (most recent trunk on Ubuntu) or its Python API to get the correct EPSG code for a given prj file WITHOUT relying on the prj2epsg webservice or any external webservice. To clarify, the solution needs to work on PRJ files that do NOT include the EPSG code within the definiton (which is, in my experience, nearly every shapefile out there) such as this one:

PROJCS["NAD_1983_StatePlane_Texas_Central_FIPS_4203_Feet",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",2296583.333333333],PARAMETER["False_Northing",9842500.0],PARAMETER["Central_Meridian",-100.3333333333333],PARAMETER["Standard_Parallel_1",30.11666666666667],PARAMETER["Standard_Parallel_2",31.88333333333333],PARAMETER["Latitude_Of_Origin",29.66666666666667],UNIT["Foot_US",0.3048006096012192]]

I'm connecting to these files using ogr_fdw in Postgres10/Postgis2.4. If there were some way to use ogr_fdw WITHOUT knowing the specific EPSG code for a given shapefile, I would accept that as an answer as well. I don't necessary need to know the source coordinate system -I'm just looking to transform to WGS84 when working with the data in postgis.

THX1138
  • 1,621
  • 3
  • 15
  • 18

2 Answers2

1

If you use GDAL there is no need to know the EPSG code of the coordinate system. Translation from the WKT contents of .prj file goes directly into Proj4 parameters without using EPSG codes in between. With your example you could (in theory, read further to see why it does not work) just have your .prj file available and run ogr2ogr as

ogr2ogr -f [your_outputformat] -s_srs your_prj_file.prj -t_srs epsg:4326 output_dataset input_dataset

You can test the conversion with gdaltransform utility http://www.gdal.org/gdaltransform.html. However, with your example the tool throws an error:

gdaltransform -t_srs epsg:4326 -s_srs unknown.prj
ERROR 6: No translation for Lambert_Conformal_Conic to PROJ.4 format is known.

I do not know projections well enough to say anything about Lambert_Conformal_Conic and why it is unknown to Proj4, but if I edit .prj to use Lambert_Conformal_Conic_SP2 instead I can use gdaltransform

gdaltransform -t_srs epsg:4326 -s_srs unknown.prj
2668812 10416336
-99.1422367719884 31.2389668371 0
user30184
  • 65,331
  • 4
  • 65
  • 118
0

If you are using ogr_fdw, you can change the foreign table geometry definition.

I have this documented in http://postgis.us/presentations/FOSS4G2017_PostGISSpatialTricks.html#/1/12

Repeated here:

ALTER FOREIGN TABLE shps.towns_polym 
ALTER COLUMN geom type geometry(geometry,26986);

Where in this case, ogr_fdw brought the table srid wrong as 4269 instead of what it should be 26986 (MA State plane feet).

Also if you have a mix of polygons/multipolygons you should change the type to something like geometry to handle both cases as shown above.

Regina Obe
  • 10,503
  • 1
  • 23
  • 28
  • But what tool do you use to programatically determine the correct projection without resorting to using prj2epsg? That's the issue I'm having. I know I can change the geometry definition -but need a way to determime the correct srid without resorting to manual lookups. – THX1138 Aug 18 '17 at 05:45
  • Okay misunderstood your question. Do you know a possible set of proj? Normally ogr_fdw will get the DATUM at least right though it screws up on the units if it's not wgs 84. – Regina Obe Aug 18 '17 at 06:13
  • I don't really have a specific set, I deal with data that is spread out across local jurisdictions in the US. I only need to transform the data TO WGS84 though -I don't need to keep the source projections in the fdw tables if that simplifies anything. – THX1138 Aug 18 '17 at 06:43