removal of +datum names from nad/epsg in r1874 has bad repercussions

[hamish]

Hi,

While debugging grass ticket # 1452 I tracked the problem down to a broken epsg file. OSGeo4w seems to be shipping the unreleased 4.7.1 version of proj4 which includes that file.

 https://trac.osgeo.org/grass/ticket/1452

In r1874 (see also #18) the +datum names were removed from the nad/epsg file (!?@@?!) and were replaced by hardcoded +towgs84 terms which do not come from EPSG.

The problem is that for any given epsg code (AKAIK) EPSG does not define which transform terms should be used, and it seems that proj4 is here silently picking one on our behalf. And the one it picks is not necessarily the correct one for the end-user to use!

Also, with the removal of the +datum term there is no way to undo or work around this (even if there were it would still require an unpleasant hack). Consider a parallel epsg.proj4 file if the standard epsg data is going to be modified/subsetted. (hmm, I'm not even sure if their data license allows to continue using their name if it has be modified..)

Especially 7-term params will be very localized, while the projection terms themselves may be more regional -- this is a problem for anyone with multiple +towgs84 and/or grid files to chose from within a particular epsg definition.

In my case it is a problem for epsg:27200 (NZ Map Grid + NZGD49), whereas before I could pick between 3-term, 7-term, or a NTv2 grid file to go with the epsg code, but now I get the 7-term imposed upon me with no way out.

please don't take our datum names away. :-(

thanks, Hamish

[warmerdam]

Hamish,

First keep in mind that many datums do not have PROJ.4 names.

I presume you were able to select a particular method for the datum shift because the internal expansion of the +datum field was added to the end of the definition instead of the beginning. Is that possible? Have you tried adding the override before the +init=epsg:xxx?

Perhaps we need to refine the rules for precidence of definitions?

Lets not through the baby out with the bath water. There was a reason for the change even though I don't exactly recall what it was.

[hamish]

Replying to warmerdam:

First keep in mind that many datums do not have PROJ.4 names.

Hi Frank,

sure. but the ones I'm struggling with right now are NAD83 and NZGD49 which are both known to proj -ld and gdal's ogr/ogr_srs_proj4.cpp.

the NAD83 one may be a GRASS+proj 4.7.1 specific issue, I'm still trying to figure out what's going on there.

a little grepping shows that the only surviving +datum= in the 4.7.1svn epsg file are WGS84, NAD27, and NAD83.

is the idea to seed some +towgs84 terms in the epsg file, where otherwise all unknown datums would just be ignored? (if so, where do the seeded +towgs84 terms come from?) and if so should all datums known to proj -ld be theoretically kept intact in the epsg file? +ellps= does not seem to be likewise expanded.

I presume you were able to select a particular method for the datum shift because the internal expansion of the +datum field was added to the end of the definition instead of the beginning. Is that possible? Have you tried adding the override before the +init=epsg:xxx?

the code I'm working on is in GRASS's new (and dare I say really really nice) location setup wizard. You can give it an EPSG code, +proj4 terms, georef'd file, or named projection + datum/ellps from a list, and if it sees that datum transform options exist for that datum (cross-checked from a GRASS table) it prompts you for which transform flavour you'd like to use. In the case of nzgd49 you get offered the 3-term, 7-term, or grid file to choose from.

• the relevant bit for this ticket is that for the EPSG code we parse the $PROJ_LIB/epsg file directly, and for the selected entry if we see a +datum= which matches one listed in our datum transform table file we check/prompt for the desired transform. No +datum in the epsg file means the user is locked into the pre-selected +towgs84 that happened to be in the epsg file. (and no way for us to identify that this has happened and unwind the decision)

as mentioned above we're also having problems with GRASS+pj4.7.1 & nad83,nzgd49 failing to detect available transforms (although nad27 seems ok for some reason), but I'll leave that problem for the GRASS ticket unless I figure out otherwise.

Perhaps we need to refine the rules for precidence of definitions?

Personally, I'm pretty happy with first-come first-served - if there are conflicting terms given at least it is easy to figure out what happened in the postmortem. If the precedence were to change I don't think it would be a very big deal to write some code in grass which builds up the final +proj4 string from key pairs in some predetermined order instead of just tacking stuff on the end as may happen now; I'd be more worried about the backwards compatibility issues that changing the order might bring.

Lets not through the baby out with the bath water.

ok, but bear with me as I'm just seeing flash before my eyes the possibility of the entire osgeo stack locked-in to a single transform for New Zealand for the lifetime of 4.7.1 (which if it makes it into debian/stable could be some years), and minorly freaking out at the prospect of having to support/explain some yet to be discovered workaround to my countrymen. :-)

There was a reason for the change even though I don't exactly recall what it was.

.. hunting the mailing list archive for +/- July 31st, 2010 .. I'm seeing some noises about TMSO, HARN, threadsafe, and on 2010-03-15 a note about getting public access to the structs in proj_api.h, but nothing which rings any bells for me.

thanks, Hamish

[hamish]

just looking through r1874, maybe it was done to homogenize the NAD83/HARN state plane stuff??

 http://thread.gmane.org/gmane.comp.gis.proj-4.devel/4786

If so I wonder if an expansion script got a little bit greedy in the process? and also, as in the above comment, would be rather sad that end users in the USA no longer get the obvious opportunity to decide the local transform terms they'd like to use. Well reasoned defaults are of course nice, but only if they come with the opportunity to easily override & let you know what has happened on your behalf.

just to double check what happens if you give +datum but do no specify explicit transform options, with a known test point:  http://apps.linz.govt.nz/gdb/index.aspx?mode=gmap&code=OUSD

45° 52' 10.20573" S  170° 30' 39.31470" E  NZGD2k  (~wgs84)
317720.88 698906.84 	N. Taieri TM  NZGD49 epsg:27231

proj 4.7.0 with no datum transform opts:

# nzgd2k (wgs84) to epsg:27231
echo "170d30'39.31470\"E 45d52'10.20573\"S" | \
  cs2cs +proj=longlat +datum=WGS84 +to \
    +proj=tmerc +lat_0=-45.86151336111111 +lon_0=170.2825891111111 \
    +k=0.99996 +x_0=300000 +y_0=700000 +datum=nzgd49

317721.81       698908.40 1.24

same with proj -ld's 7-term +towgs84 (identical result to above)

   +towgs84=59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993
317721.81       698908.40 1.24

same using 3-term transform: (a few meters away)

   +towgs84=54.400,-20.100,183.10
317717.61       698908.87 -32.92

same using the NTv2 grid file: (within 4cm of published value)

   +nadgrids=nzgd2kgrid0005.gsb
317720.84       698906.83 0.00

so it would seem with 4.7.0 that not explicitly listing the transform opts with a (known) datum incorporates the listed default transform & all works as expected.

cheers, Hamish

[Alister]

Hi Hamish, Do you think there is a problem with gdal/ogr as well, or just Proj? After building Proj 4.7.0 and rebuilding gdal I would have thought I could do something like this, and it would automatically use the natgrid, but it doesn't:

ogr2ogr -f "ESRI Shapefile" /c/coord_test/test /c/coord_test/NZTM.shp -T_SRS 
EPSG:27200

(BTW I'm actually trying to figure out how to make QGIS use the NTV2 grid for NZMG. In the past I *think* their CRS database was manually edited to add the +natgrids options, but they are now automatically generated from gdal).

[Alister]

they are now

Sorry, I mean "it is now"

[Alister]

Hmmm. I was a little confused when I posted that. Sorry if it seems like I'm spamming this ticket :(

What I'm concerned about is actually a different issue, i.e. rather than "why does proj only allow one default datum transform?", "why do proj and gdal use the 7-term +towgs transform by default, when the NTV2 grid is more accurate?". Wouldn't it be better to use the NTV2 grid by default if it is available? That's what QGIS did for some time, unfortunately now it automatically uses the gdal defaults, and this is causing all sorts of confusion because the results are so much less accurate than other software that people might also be using.

[neteler]

Not even EUR50 Datum is recognized any more with PROJ 4.7.1 :( The removal of +datum= breaks the GRASS location wizard among other things. What should we as PROJ users do?

Markus

[neteler]

For the record, some related postings from the list:

 http://osgeo-org.1803224.n2.nabble.com/European-Datum-1950-is-not-saved-in-the-region-information-td6958890.html#a6962707

[warmerdam]

For the record, the change in OGRSpatialReference was done 20 months ago and is revision:

 http://trac.osgeo.org/gdal/changeset/19903/trunk/gdal/ogr/ogr_srs_proj4.cpp

Related to ticket:

 http://trac.osgeo.org/gdal/ticket/3450

I see there is already a CPL config options to use +datum, so I should be able to regenerate the PROJ.4 epsg init file with no change to GDAL. I need to run now, I'll try to address this tonight.

[warmerdam]

Regenerated using:

  epsg_tr.py --config OVERRIDE_PROJ_DATUM_WITH_TOWGS84 FALSE -proj4 -skip -list gcs.csv > epsg
  epsg_tr.py --config OVERRIDE_PROJ_DATUM_WITH_TOWGS84 FALSE -proj4 -skip -list pcs.csv >> epsg

This is also now noted in the libgeotiff/csv/README (see  http://trac.osgeo.org/geotiff/changeset/2171).

I'm closing this for now, and beta2 will have this reverted file.

[warmerdam]

I should have noted that the epsg update was r2172.