reprojection using ESRI .prj files is not accurate - missing TOWGS84

[etourigny]

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

For example:

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

[etourigny]

example file

[etourigny]

other example, point data

[etourigny]

found a bug in importFromEPSG(), if the target code is not found this statements resets eErr, therefore the calling function does not know it was not found, and an empty SRS is used.

    eErr = FixupOrdering();

For example:

$ gdalsrsinfo   EPSG:327504

ERROR 6: No translation an empty SRS to PROJ.4 format is known.
PROJ.4 : ''

OGC WKT :

after the fix:

$ gdalsrsinfo   EPSG:327504
ERROR 1: ERROR - failed to load SRS definition from EPSG:327504

fixed in r23396

[etourigny]

Added script and generated files needed by EPSG code detection, these files must reside where CPLFindFile() can find them (e.g. with the other GDAL data files).

[etourigny]

added experimental support (r23397) for EPSG detection in gdalsrsinfo, use either -e or -o epsg to activate search.

should add an autotest using this data when is stable.

One way to test, try to get the following:

$ gdalsrsinfo -o proj4  brazil.prj
'+proj=longlat +ellps=aust_SA +no_defs '

$ gdalsrsinfo -o proj4  -e brazil.prj
'+proj=longlat +ellps=aust_SA +towgs84=-57,1,-41,0,0,0,0 +no_defs '

$ gdalsrsinfo -o proj4  -e 'GEOGCS["SAD69",DATUM["D_South_American_1969",SPHEROID["GRS_1967_Modified",6378160,298.25]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]'
'+proj=longlat +ellps=aust_SA +towgs84=-57,1,-41,0,0,0,0 +no_defs '

[etourigny]

Attaching a patch that adds the functionality to OGRSpatialReference.

New Functions: FindEPSG() and SearchDictForWKT() (also c api versions).

importFromESRI() searches for EPSG code and imports from that if OGR_ESRI_FIND_EPSG=YES .

Even, can you please check this out for blatant errors and legal issues? I am particularly curious to know: if we bundle these files, with some modifications to the EPSG database (coming from *.override.csv), does that violate EPSG? (c.f. PROVENANCE.TXT)

I assume that this would require an RFC (if the default behavior is changed), so how should I proceed? Should I write to the list, go through RFC process then commit if accepted? Or perhaps I can commit the code, then go through RFC process to make it default?

[Even Rouault]

Hum, this is a bit of a tricky piece of code with various implications, so I believe that you should not commit the code before it makes some consensus, being under the form of a formal RFC or not. At least, on such a matter, I would not been confident in proceeding without a +1 coming from FrankW, for several reasons : FWIW, he's the PSC member with the most background on projection issues, and particularly he has regularly updated the EPSG database in the whole proj4/libtiff/GDAL stack.

My observations :

1) We have already a AutoIdentifyEPSG() method with, admittedly, several limitations, but I don't see a strong reason for creating a new one. AutoIdentifyEPSG() could be enhanced with your proposal.

2) The new files are generated from a script, so this raises the question on how and when we produce them. There are several possibilities that come to mind :

• at installation time, in the target release of the makefile (implication for Windows builds to think about) ?
• as a complementary script that lives in GDAL, and that is run each time the EPSG database is refreshed. In which case the resulting files are placed in SVN
• pretty much the same as the previous one, but the script is integrated in the whole machinery used to update the whole proj4/libtiff/GDAL stack.

3) The legal implications are out of my area of expertise. AFAIR, EPSG has never protested - at least publicly - against the use we make of their data. Perhaps Frank has some interesting feedback on this ? It might be perhaps prudent not to put "epsg" in the name of the files, especially the ones that have gone through the "wkt_simple" and "wkt_esri" steps, to avoid implying that their data reflect the EPSG database.

4) I have issues to follow the rationale of the end of the patch.

if( oSRSTemp.importFromEPSG( nCode ) == OGRERR_NONE )
{
     if( oSRSTemp.morphFromESRI() == OGRERR_NONE )

--> weird: importFromEPSG() will create a "OGC WKT" SRS. So it makes non sense to do morphFromESRI() on it, because it is not ESRI WKT.

Similarly, just after

eErr = this->morphFromESRI();
if( eErr == OGRERR_NONE )
    eErr = this->importFromEPSG( nCode );

Why do you start with morphFromESRI(). At that point of the code (since if ( nCode != -1 )), you have identified the SRS as being an EPSG one. Why not directly calling importFromEPSG() ?

So all in all, I believe that the following should be enough :

519	+            if ( EQUAL(CPLGetConfigOption("OGR_ESRI_FIND_EPSG","NO"),"YES" ) )
520	+            {
521	+            int nCode = -1;
522	+            nCode = this->FindEPSG();
523	+            if ( nCode != -1 )
524	+            {
                         eErr = this->importFromEPSG( nCode );
554	+            }

If FindEPSG() succeeded, I can't see how importFromEPSG() could fail.

Am I missing something ?

5) I have not taken enough time to fully understand the implications and limits of how you do the actual matching

6) If you do only pure string equality checks, you could just store some hashing of the string (long enough to avoid collisions, md5 might be far sufficient). This would make the files smaller, and comparison a bit faster. But that's perhaps premature optimization ! And if partial string matching is implemented, this won't fly.

[warmerdam]

I will try to look at this more closely tonight. But briefly I think when we do the FromESRI() conversion we should try to lookup the ESRI datum name in gdal_datum.txt and add the preferred TOWGS84() parameters for the datum if there is a match.

I would not like to have more supporting files if it can be avoided. But we might need to refresh the ESRI datum name generation process a bit for gdal_datum.txt.

[etourigny]

Replying to warmerdam:

I will try to look at this more closely tonight. But briefly I think when we do the FromESRI() conversion we should try to lookup the ESRI datum name in gdal_datum.txt and add the preferred TOWGS84() parameters for the datum if there is a match.

I would not like to have more supporting files if it can be avoided. But we might need to refresh the ESRI datum name generation process a bit for gdal_datum.txt.

This approach would surely work to solve the issue at hand. If I understand correctly, this could work like (using EPSG:4618): fetch datum code from gdal_datum.csv (6618), find matching SRS code in gcs.csv (4618) and associated preferred datums. What would you need to change in gdal_datum.csv?

However, when working on this I thought it was a neat feature to be able to do a "reverse-EPSG" lookup for any SRS (that I implemented into gdalsrsinfo). For a WKT with all authority nodes, it's a trivial task, but when the authority nodes are absent the obvious solution is the one I implemented. The first solution I thought, to loop over all EPSG codes in pcs.csv and gcs.csv and compare various WKT variations to the target, is not efficient (up to 3 seconds in worst-case).

I understand that you want to limit support files, is there a way that you think this can be achieved in an efficient manner without them?

[etourigny]

Replying to rouault:

Hum, this is a bit of a tricky piece of code with various implications, so I believe that you should not commit the code before it makes some consensus, being under the form of a formal RFC or not. At least, on such a matter, I would not been confident in proceeding without a +1 coming from FrankW, for several reasons : FWIW, he's the PSC member with the most background on projection issues, and particularly he has regularly updated the EPSG database in the whole proj4/libtiff/GDAL stack.

My observations :

1) We have already a AutoIdentifyEPSG() method with, admittedly, several limitations, but I don't see a strong reason for creating a new one. AutoIdentifyEPSG() could be enhanced with your proposal.

yes, but I would also like to have this function available, for a "reverse-EPSG" lookup. This function could call a revised AutoIdentifyEPSG() and return the last autority node.

2) The new files are generated from a script, so this raises the question on how and when we produce them. There are several possibilities that come to mind :

• at installation time, in the target release of the makefile (implication for Windows builds to think about) ?
• as a complementary script that lives in GDAL, and that is run each time the EPSG database is refreshed. In which case the resulting files are placed in SVN
• pretty much the same as the previous one, but the script is integrated in the whole machinery used to update the whole proj4/libtiff/GDAL stack.

this would be preferable, as it is rather time-consuming (in its present non-optimized fashion).

3) The legal implications are out of my area of expertise. AFAIR, EPSG has never protested - at least publicly - against the use we make of their data. Perhaps Frank has some interesting feedback on this ? It might be perhaps prudent not to put "epsg" in the name of the files, especially the ones that have gone through the "wkt_simple" and "wkt_esri" steps, to avoid implying that their data reflect the EPSG database.

right. hope I didn't do a blunder in calling these files in the svn version of "gdalsrsinfo"

4) I have issues to follow the rationale of the end of the patch.

if( oSRSTemp.importFromEPSG( nCode ) == OGRERR_NONE )
{
     if( oSRSTemp.morphFromESRI() == OGRERR_NONE )

--> weird: importFromEPSG() will create a "OGC WKT" SRS. So it makes non sense to do morphFromESRI() on it, because it is not ESRI WKT.

Similarly, just after

eErr = this->morphFromESRI();
if( eErr == OGRERR_NONE )
    eErr = this->importFromEPSG( nCode );

Why do you start with morphFromESRI(). At that point of the code (since if ( nCode != -1 )), you have identified the SRS as being an EPSG one. Why not directly calling importFromEPSG() ?

So all in all, I believe that the following should be enough :

519	+            if ( EQUAL(CPLGetConfigOption("OGR_ESRI_FIND_EPSG","NO"),"YES" ) )
520	+            {
521	+            int nCode = -1;
522	+            nCode = this->FindEPSG();
523	+            if ( nCode != -1 )
524	+            {
                         eErr = this->importFromEPSG( nCode );
554	+            }

If FindEPSG() succeeded, I can't see how importFromEPSG() could fail.

Am I missing something ?

I was being overly cautious there, and braindead in writing morphfromesri()

5) I have not taken enough time to fully understand the implications and limits of how you do the actual matching

6) If you do only pure string equality checks, you could just store some hashing of the string (long enough to avoid collisions, md5 might be far sufficient). This would make the files smaller, and comparison a bit faster. But that's perhaps premature optimization ! And if partial string matching is implemented, this won't fly.

I used the simplest method as a first iteration (string comparison).

The idea is to create a simple version (wkt_simple - the WKT you would get with StripCT()) of a WKT, so that we can dumb-down the target WKT and find the closest match. As there seems to be a difference in significant digits between the ESRI KWT and OGC WKT generated by GDAL, I preferred to try both wkt_esri and wkt_simple.

If the precision issue is solved, then using (md5) hashes would be sufficient (I think). Even better, we could add an extra column in pcs.csv and gcs.csv (and others) with the hash of the simple (or esri) wkt. On the other hand, a single file (wkt_simple or wkt_esri) would not be that big (about 100kb compressed). I must stress that this technique is very efficient, it doesn't cause any slow down in the test I did with gdalsrsinfo, so the benefit of using a hashed (smaller) file would not be that great.

[etourigny]

attaching a patch for OGRSpatialReference::morphFromESRI() implementing Frank's suggestion. It checks for matching SPHEROID (SAD69 has 2 datums - 6618 and 6291), and then adds the TOWGS84 node.

Is this approach reasonable?

I need to implement this for PROJCS also, any ideas on where I can find a number of GEOGCS and PROJCS .prj files?

[Even Rouault]

Overall, looks reasonable to me.

1) A suggestion to improve code readability, because I had a hard time to understand it. Don't take it personnaly, you are not responsible for it, but the current practice of using + 0, + 1 et + 2 for using the "members" of papszDatumMapping is awful. I'd suggest to :

/* DM is for Datum Mapping */
#define DM_IDX_EPSG_CODE            0
#define DM_IDX_ESRI_NAME            1
#define DM_IDX_MASSAGED_EPSG_NAME   2
#define DM_ELT_SIZE                 3

#define DM_GET_EPSG_CODE(map, i)          map[(i)*DM_ELT_SIZE + DM_IDX_EPSG_CODE]
#define DM_GET_ESRI_NAME(map, i)          map[(i)*DM_ELT_SIZE + DM_IDX_ESRI_NAME]
#define DM_GET_MASSAGED_EPSG_NAME(map, i) map[(i)*DM_ELT_SIZE + DM_IDX_MASSAGED_EPSG_NAME]

and use them :-) Retrofitting existing code is an extra bonus :-)

2) I couldn't believe that there were different datums with sane name ! But indeed you are true. So the test for if( EQUAL( pszTemp,GetAttrNode( "DATUM|SPHEROID" )->GetChild(0)->GetValue())) is justified. It could be nice to add a comment in that section of the code to explain why it is also necessary to compare spheroid names.

3) Dereferencing GetAttrNode( "DATUM|SPHEROID" ) without NULL-checking could segfault in situations where the string passed to morphFromESRI has no spheroid node. Not sure this is valid WKT, but it is always good if we are robust to invalid input.

As far as collecting .prj files is concerned, perhaps you could launch a "Call to users : request for collecting shapefile .PRJ files" (better worded title welcome ;-)) on the mailing list ? Explain the purpose and invite them to attach their files, or links to datasources, to a wiki page (or another more appropriate place) for example.

Ah, and I wanted to mention that you can join #gdal on irc.freenode.net for more interactive exchanges.

[etourigny]

fixed bug in importFromEPSG(when code is not found) in 1.8 (r23422).

This regression was brought by fixes for bug #4178: Fixed importFromEPSGA() in trunk (r22819) and 1.8 (r22820) to always call FixupOrder?() at the end of the method.

[bishop]

Some prj files for testing

[Even Rouault]

I'm attaching a text file with various URLs that I've found with google searches. There are some .prj that don't match EPSG code in the list, or some unusual formulations of .prj that are likely invalid ESRI WKT, and at the end a link to ESRI resources that are very similar to the ones that Martin Daly mentionned. There are also a few interesting variations of the same SRS with different precision for floating point values...

[Even Rouault]

Various list of URLs with shapefiles, and their .prj files

[etourigny]

new prj file , aea with SAD69 (from IBGE)

[etourigny]

more files from IBGE (​ftp://geoftp.ibge.gov.br/mapas/)

[etourigny]

Attaching a patch that works with the files in ibge.zip and test_cs.zip . Haven't yet tested with Even's suggestions.

The patch uses a config option GDAL_FIX_ESRI_WKT=TOWGS84/DATUM/GEOGCS/NO (default TOWGS84). Based on that option the modifications to the WKT are more or less important.

I made a small modification to support Krassowsky/Krasovsky spheroid (used in some files in test_cs.zip).

I'll make a python test function that reads a file (with one ESRI WKT per line) to test adding TOWGS84 parameters and doing the round-trip conversion. This should help validate the code and also find other errors in this code and also morph to/from ESRI. I think the files from ESRI can serve as a good starting point.

Let me know what you think. This code could be linked in some way to AutoIdentifyEPSG().

[etourigny]

default for GDAL_FIX_ESRI_WKT has to be set to "NO" for autotest osr_esri.py

const char *pszFixWktConfig=CPLGetConfigOption( "GDAL_FIX_ESRI_WKT", "TOWGS84" ); const char *pszFixWktConfig=CPLGetConfigOption( "GDAL_FIX_ESRI_WKT", "NO" );

[etourigny]

r23542 : add optional fixing of TOWGS84, DATUM and GEOGCS

Implemented a slightly modified fix (also checks for identical PRIMEM for bug #4378).

For now is not active by default, must use config option GDAL_FIX_ESRI_WKT=YES/TOWGS84/DATUM/GEOGCS.

My feeling is that we can safely set the default to DATUM, but I don't want to do that without reflexion and consulting the community.

[etourigny]

ESRI prj file for which gdalsrsinfo -e does not work

[etourigny]

Attached 102067.prj - ESRI prj file for which gdalsrsinfo -e does not work (got it from spatialreference.org).

The definition for ESRI code 102067 is in esri_extra.wkt, and even if this file is added to the wkt files to scan for in gdalsrsinfo.cpp, it is not detected, for the following reasons:

• an AUTHORITY node which is present in esri_extra.wkt, but not in the .prj file
• datum names do not match (D_Jednotne_Trigonometricke_Site_Katastralni vs. D_S_JTSK)

It is possible the .prj file I got is not correct - jmckenna please upload your file so I can test it.

[etourigny]

The files needed for 'gdalinfo -e <srs_def>' are: epsg_esri.wkt.gz epsg_ogc_simple.wkt.gz epsg_ogc.wkt.gz .

Attaching a modified gen_epsg_wkt.sh, which fixes '-o wkt_simple' error.

Also detects null strings returned by gdalsrsinfo when EPSG code is not found or not completely supported. Some codes are now ignored when they are not fully supported by OGR, such as EPSG:2221 (EPSG: No WKT support for projection method 9826). In these cases, importFromEPSGA() returns OGRERR_UNSUPPORTED_SRS and we have no way of knowing how serious the error is. This should be fixed eventually in ogr_fromepsg.cpp, to at least create a full WKT even if the projection method is not supported.

[etourigny]

script to generate epsg*.wkt files

[hamish]

Hi,

fwiw, some basic summary/background about the issue (or rather the other way, why ESRI WKT morphing strips off TOWGS84) can be found in this ticket:

​https://trac.osgeo.org/grass/ticket/1988

I still don't know how explicitly defined NTv2 grid files should be handled during conversion of the SRS into (non-ESRI/any) WKT-ese, but would like to.

Hamish

[Even Rouault]

Should be mostly solved with r40207

$ gdalsrsinfo -e 'ESRI::GEOGCS["SAD69",DATUM["D_South_American_1969",SPHEROID["GRS_1967_Modified",6378160,298.25]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]'

EPSG:4618

PROJ.4 : +proj=longlat +ellps=aust_SA +towgs84=-66.87,4.37,-38.52,0,0,0,0 +no_defs 

OGC WKT :
GEOGCS["SAD69",
    DATUM["South_American_Datum_1969",
        SPHEROID["GRS 1967 Modified",6378160,298.25,
            AUTHORITY["EPSG","7050"]],
        TOWGS84[-66.87,4.37,-38.52,0,0,0,0],
        AUTHORITY["EPSG","6618"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4618"]]