Ticket #4345 (assigned defect)
reprojection using ESRI .prj files is not accurate - missing TOWGS84
|Reported by:||etourigny||Owned by:||etourigny|
|Severity:||normal||Keywords:||reprojection esri shapefile|
|Cc:||rouault, warmerdam, jmckenna|
When re-projecting ESRI shapefiles that have a .prj. file with ESRI WKT, the re-projection can be inaccurate as .prj files lack TOWGS84 parameters.
This problem has been discussed on the gdal-dev mailing list subject "ogr2ogr reprojection, features are not transformed".
This does not happen if the SRS has WGS84 datum, but happens for other datums (e.g. SAD69).
$ cat brazil.prj GEOGCS["SAD69",DATUM["D_South_American_1969",SPHEROID["GRS_1967_Modified",6378160,298.25]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] $ gdalsrsinfo -o proj4 brazil.prj '+proj=longlat +ellps=aust_SA +no_defs ' $ gdalsrsinfo -o wkt -p brazil.prj GEOGCS["SAD69", DATUM["South_American_Datum_1969", SPHEROID["GRS_1967_Modified",6378160,298.25]], PRIMEM["Greenwich",0], UNIT["Degree",0.017453292519943295]] $ gdalsrsinfo -o proj4 EPSG:4618 '+proj=longlat +ellps=aust_SA +towgs84=-57,1,-41,0,0,0,0 +no_defs ' $ gdalsrsinfo -o wkt -p EPSG:4618 GEOGCS["SAD69", DATUM["South_American_Datum_1969", SPHEROID["GRS 1967 Modified",6378160,298.25, AUTHORITY["EPSG","7050"]], TOWGS84[-57,1,-41,0,0,0,0], AUTHORITY["EPSG","6618"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4618"]]
The workaround is to set -s_srs to the desired EPSG code (in the case of SAD69 EPSG:4618).
For example, the following produces the incorrect transformation (there is in fact no transformation at all) of a file from SAD69 to WSG84:
ogr2ogr -t_srs EPSG:4632 tmp1.shp brazil.shp
Specifying the EPSG code gives a proper transformation:
ogr2ogr -s_srs EPSG:4618 -t_srs EPSG:4632 tmp2.shp brazil.shp
There are 2 problems with this
1) This does not fix the underlying problem, and it can creep up elsewhere. For example, when the original file (in SAD69 datum) is viewed in QGIS, there is a shift (relative to the file in WGS84) due to absent TOWGS84 in the .prj files (it uses CRS '+proj=longlat +ellps=aust_SA +no_defs' - without +TOWGS84)
2) It would be better if one did not have to set the EPSG code every time.
A solution to 1) and 2) would be to automatically translate the ESRI WKT definition to the one given by the appropriate EPSG code (with TOWGS84).
This can be done by searching all EPSG codes that GDAL/OGR knows (in pcs.csv, gcs.csv and others), and finding the one that matches the input ESRI WKT.
From the mailing list, here are some ideas:
I have been playing around a bit and here is what I did that works (first try): - take a given CRS definition (from say EPSG or .prj file) and find it's ESRI WKT or "simple" WKT. - for all the EPSG codes in pcs.csv and gcs.csv, get it's ESRI (or simple WKT), and compare that to the target WKT - if you've a matching WKT, then get the full WKT corresponding to the EPSG code that matches. The problem is that it's pretty inefficient as you can imagine, taking a few seconds to find one single target. A second iteration: - generate full WKT, ESRI WKT and "simple" (StripCT) WKT for all EPSG codes in pcs.csv and gcs.csv - save these to a flat (gzipped) file in csv form - use these tables to find the EPSG code that matches a given WKT (in whatever WKT flavor you need) This is rather efficient in terms of processing time. [...] This works for all EPSG codes I tried (think of it as a reverse EPSG lookup), and also a few .prj files. A problem I encountered was the differences in significant digits in the ESRI-WKT and OGC-WKT, so for now it works best if warping to ESRI WKT. [...] Should I create a sandbox for an experimental gdalsrsinfo util implementing this idea?
Problems to solve are:
- efficient scanning requires a pre-built database of known WKT definitions, can this reside in the GDAL data tree?
- accurate conversion between ESRI and OGC WKT, as done in morphToESRI/morphFromESRI
- dealing with differences in significant digits for numerical values.
Also, I have written code in gdalsrsinfo that does this lookup (when given option '-o epsg'). Could I commit that to trunk for testing or should I setup a sandbox?
The final solution would be to incorporate this into the ESRI shapefile driver (with a new function in OGRSpatialReference) so it would be used when re-projecting (in ogr2ogr but also when viewing the data).