GPX tracks display in different position depending on GPX track setting for NAD83 map (regression)

[sfroyen]

Steps to reproduce

1. Open a map on Android OOM and attach a GPS.
2. Go for a walk and sprinkle some point objects using GPS positioning
3. Copy the map from Android to a desktop and open in OOM
4. Observe that yhe point objects added using GPS coordinates no longer align with the GPX track from OOM Android.

Here's a screenshot from Android:

And here's one from the Mac:

Notice that the track position where I added the purple plus sign shifts to the NE by about 1m. This shift is similar to the shift between NAD83 and WGS84 GPS lat/lon coordinates. Note that the shift only happens to the GPX track, not to the map objects that were added.

(I can close and restart Android OOM and the track does not move.)

Caveats:
The map I'm using was created using NAD 83 / UTM zone 13N and then converted to EPSG:6342 / NAD83(2011) when that option was added.

Mapper Version: 0.9.3 and v20200613.1 (I've also tried 0.9.3)
Operating System: Android 7.0 and MacOS 10.15.6

[dg0yt]

There is as a "GPX" setting in GDAL > Import with GDAL/OGR. Is it different on macOS vs. Android?

[pkturner]

Could they be built with different versions of PROJ?

[dg0yt]

Could they be built with different versions of PROJ?

No, not for a single official release (v0.9.3) on Android, macOS, Windows. These binaries come out of Continuous Integration from the same source tarballs.

There is as a "GPX" setting in GDAL > Import with GDAL/OGR. Is it different on macOS vs. Android?

This is important, because this represents different code paths.

[sfroyen]

Android OOM did not have GPX checked under GDAL/OGR. Once I checked that option, Android OOM behaves like the desktop version, i.e., the (saved) track is shifted. I will repeat the exercise with a new capture and report back.

My desktop OOM had the GDAL GPX box checked.

[sfroyen]

FYI - Unselecting the GDAL GPX option on the desktop and importing the GPX track again produces the initial Android result -- as expected.

[dg0yt]

Okay, I was expecting this from line style on the screenshots.
When checked, i.e. when using GDAL, coordinate conversion is done by GDAL. The other option does not use GDAL.

[sfroyen]

Two bugs?

Collecting a track on Android with the GDAL GPX option on behaves as follows. While the track is being collected, it aligns with the position of added point objects. After saving and reopening the map, the GPS track is shifted to the NW.

Re. my previous FYI:
After importing the GPX track on the desktop with the the GDAL GPX option off, there is no way to re-import it with the GPX option on. It seems that OOM caches the previous import and uses the non-GDAL method. Clearing the Closed Template list had no effect.

[dg0yt]

After importing the GPX track on the desktop with the the GDAL GPX option off, there is no way to re-import it with the GPX option on. It seems that OOM caches the previous import and uses the non-GDAL method. Clearing the Closed Template list had no effect.

I have a GPX test file with a named waypoint. Waypoint names are displayed only by the non-GDAL implementation. With this test file, I see that the chosen implementation correctly follows the user setting when opening a map file containing a GPX template, or when adding a new GPX file.

[dg0yt]

I can reproduce a small shift of app. 1 m at 40N 105W with PROJ 7.1.

I started with a map and GPX track created with PROJ 4.9 / GDAL 2 which did not have alignment issues. When I opened the map with Mapper using PROJ 7.1, there was an immediate warning about a georeferencing mismatch of app. 2 m, i.e. the geographic reference point and the projected reference point stored in the georefencing did not match, according to PROJ 7.1. This is what I get the for the same geographic coordinates:

Geographic coordinates	40N 105W
EPSG:6342, PROJ 4.9	500000.00 4427757.22
EPSG:6342, PROJ 7.1	500000.86 4427756.50

(Sidenote: Mapper assumes "WGS 84" for geographic coordinates, but this is not accurate with regard to different realizations. And PROJ switch from "WGS 84" to "GRS 1980" for some purposes.)

Now when I fix the georeferencing, keeping geographic coordinates, I would expect that the georeferenced GPX template keeps its position relative to the map, because GPX data uses geographic coordinates. This seems to work for non-GDAL mode, but not for GDAL mode.

[dg0yt]

The table shows results from Mapper's georeferencing dialog, i.e. from the transformation which Mapper establishes with PROJ. So how does this work?

Mapper constructs a PROJ transformation from +proj=latlong +datum=WGS84 to +init=epsg:6342 (using proj_create_crs_to_crs). Note the source CRS. Changing that to EPSG:4326 (for an accurate modern use of PROJ) makes Mapper produce the original result which we had seen before.

To finally overcome these issues, I would like to fully establish modern usage of PROJ, at least for EPSG.

[sfroyen]

PROJ 4 did not support the necessary transformations and the EPSG files did not have the necessary parameters until late 2019 or early this year. Also, it appears that EPSG:4326 (or wgs84) is some form of generic specification (I have not figured out exactly what that means). Perhaps the current revision at current time?).

EPSG 7657, 7659, 7661, 7663 and the current 7665 have been used by GPS satellites between 1994 and present. A very precise GPS (without correction signals) should have returned coordinates in one of the these systems (EPSG:7665 at present). I believe their parameters include linear corrections for time so, for centimeter precision, we need a timestamp too. If you use SBAS corrections, you should get the same coordinates with parameters updated twice(?) per year. If you do not have a super accurate GPS, you can use RTK corrections from a base station to get cm precision. In this case, the lat/lon coordinates depend on the base station setting. They may be WGS based, e.g., EPSG:7665 or NAD83 based like EPSG:6318. The two are not the same:

echo 40.02020698N 105W | cs2cs -d 8 EPSG:7665 EPSG:6318
40.02020047 -104.99998987 0.88568757

My GPS was using WGS based EPSG:7665, so comments about NAD83 based lat/lon do not apply here. But any solution this issue ought to handle the NAD83 case too. Perhaps geographic coordinates should not be used internally and always be converted from/to on import/export.

I will try to investigate PROJ's behavior when using wgs84 and EPSG:7665.

[dg0yt]

EPSG:4326 is WGS 84 in terms of EPSG.

The GDAL GPX driver uses WGS 84 defined as EPSG:4326 with explicit lat,lon order for geographic 2D coordinates, via macro SRS_WKT_WGS84_LAT_LONG. Many other GDAL drivers do the same. SRS_WKT_WGS84_LAT_LONG is public API:
https://github.com/OSGeo/gdal/blob/e67def79aec4f3354b093b35ddbf956c43812ffb/gdal/ogr/ogr_srs_api.h#L70-L71

By using EPSG:4326 or WGS 84 instead of +proj=latlong +datum=WGS84, we would effectively use the same default as GDAL.

Choosing something else than the default when needed is a different thing. We don't get exaxt information from GPX files or GNSS sources. (But we do have timestamps.)

[pkturner]

Using EPSG:4326 with PROJ is inaccurate in North America, because it sloppily "converts" WGS 84 (the datum of most GPS receivers) to NAD83 with a no-op. The shift is on the order of a meter.

EPSG:7665 (3D) and EPSG:9057 (2D) are more specific and current realizations of WGS 84 that PROJ converts accurately to/from NAD83. I hope Mapper can make it easy to get that accuracy.

[pkturner]

Another consideration -- and at first I thought this was the issue at hand -- is that after Mapper first displays a template with the map, it disrupts the user if the template is opened later and is positioned differently. Should tracks and other templates be defended against changes in PROJ versions, in GDAL options, in Mapper's implementation of Geographic CRS, etc.?

[dg0yt]

Trying to avoid disruption is certainly something what we try to achieve. For this issue, switching from +proj=latlong +datum=WGS84 to EPSG:4326 or WGS 84 is just that - maybe restoring an older shortcoming. I didn't test yet with EPSG:9057, but I will give it a try.

To make the template positioning reproducible, you need to store the exact definition of the CRS of the external data, and the exact definition of the transformation to the internal representation (the projected CRS of the map.) Mapper does this only partially, and still using old PROJ syntax, which leads to a number of issues.

And our unit tests still don't cover the NAD83 world, to detect bugs related to our usage of PROJ or GDAL.

[dg0yt]

I didn't test yet with EPSG:9057, but I will give it a try.

So different options for the (hard-coded) map's geographic CRS in Mapper give the following results:

PROJ	Geographics CRS in Mapper	40N 105W @ EPSG:6342	GPX track offset	Comment
4.9.3	+proj=latlong +datum=WGS84	500000.00 4427757.22	no	original behaviour
7.1.0	+proj=latlong +datum=WGS84	500000.86 4427756.50	yes	reported regression
7.1.0	EPSG:4326	500000.00 4427757.22	no
7.1.0	WGS84	500000.00 4427757.22	no
7.1.0	EPSG:9057	500000.86 4427756.50	yes

So only EPSG:4326 or WGS84 fix the regression directly. Probably, because Mapper use the same geographic CRS as GDAL does. Still, we can't change our geographic CRS easily because there will be a new inconsistency: The data loaded from the file will have a pair of reference points (geographic/projected) which doesn't match the current transformation.

If we want to go with EPSG:9057, we need to use EPSG:9057 for GPX with GDAL, overriding GDAL's default behaviour. And maybe also for other formats which use SRS_WKT_WGS84_LAT_LONG.

Anyway, we have now two deal with two types of maps, both declared as +proj=latlong +datum=WGS84: The ones created with older version of PROJ, and the ones created with 7.1.0. Maybe we need to select either EPSG:4326 or EPSG:9057, depending which one gives the best match for the georeferencing setup, and then saving and using that code, creating the third type of maps...

[sfroyen]

I suggest leaving the choice of geographic CRS up to the user and allowing for multiple CRSes.

Since I have not examined the code, I may be suggesting features that have already been implemented.

It seems to me that map coordinates and projected coordinates are the only CRSes needed inside Mapper. Map coordinates are a (very useful) convenience and and projected coordinate CRS should be established by the data used to create the base map , e.g., LIDAR. If no projected data is used, a common projection should be chosen. In the US, at the moment, that might be UTM NAD83(2011) or NAD83 PA11 plus zone.

A geographic CRS is only needed for import and export of lat/lon data into and out of Mapper. There may be multiple CRSes for sources, including unknown. On import, Mapper should prompt for CRS, both for templates (as is done now) and for live GPS data. If Mapper attaches the CRS to the imported data, an export of, e.g., a GPX track, could use the CRS that was specified on import. This should solve the current issue.

Export of non-georerenced but carefully positioned image templates could be added as a Mapper enhancement.

Reporting/displaying positions in lat/lon is just another export.

[pkturner]

@sfroyen
Thanks for the suggestions. A user's perspective is always helpful.

map coordinates and projected coordinates are the only CRSes needed inside Mapper

When a template has projected coordinates of its own, possibly using a different CRS from the map (e.g. OpenStreetMap), Mapper sets up and keeps a matrix and parameters that transform the template's elements for display.

GPX tracks are lat/lon data using a geographic CRS which is not specified in the data file. If, as you suggest, Mapper were to prompt for the CRS for live GPS data, Mapper would naturally retain that CRS in its template metadata -- it already has a field for that purpose. Perhaps that's exactly what you had in mind. In my opinion, it would be tedious as a user to select a CRS every time I start to record a track. Most days my map and my GPS receiver do not change, so the appropriate CRS doesn't change either.

A geographic CRS is not flat and thus is is not suited to a matrix type of template transform. Mapper potentially manages a variety of projected CRSes, as well as actively transforming GPX tracks based on their geographic CRS.

I like the idea that the user be given an option of what geographic CRS is used to display lat/lon. The would let Mapper off the hook from needing to guess what geographic CRS best fits the map's projected CRS.

[dg0yt]

It seems to me that map coordinates and projected coordinates are the only CRSes needed inside Mapper.

I more or less agree that we only need a well-defined projected CRS.

I like the idea that the user be given an option of what geographic CRS is used to display lat/lon. The would let Mapper off the hook from needing to guess what geographic CRS best fits the map's projected CRS.

Actually I would like to simply take the geographic CRS the projected CRS is based on. And use a direct pipeline between CRS everywhere else (#1598).

Still, we have to carefully deal with old maps, and with GDAL drivers. And again, the NAD83 world is not covered by unit test data yet.

[pkturner]

simply take the geographic CRS the projected CRS is based on

I just now saw that PROJ 7.1 supports that. It would be interesting to see what it produces, e.g. Is the geographic CRS for EPSG:6342 identified as EPSG:6318?

the NAD83 world is not covered by unit test data

I'd be willing to enhance tests in that direction. In georeferencing_t.cpp, an EPSG:6342 example would be added. To satisfy the accuracy discussed above, the max_dist_error would be reduced to 1.0 meter or smaller. Perhaps it should expect the behavior provided by PROJ 4.9, and the PROJ 7.1 behavior would be a bug.

For template_t.cpp, the tests which most need enhancement for NAD83 would be templateTrackTest and ogrTemplateTest. IMO if those tests were given a map based on EPSG:6342 along with data with a recent WGS84 GPS track, the correct result is what PROJ 7.1 provides. I expect that for Mapper to handle this case well, the crs_spec attribute of GPX templates will need to be changed to something like EPSG:9057. Of course Mapper doesn't support that right now.

I'm just throwing out some thoughts on how NAD83 would be covered by tests, to get the ball rolling. @dg0yt please feel free to tell me to hold off, or point me in a different direction.

[sfroyen]

GPX tracks are lat/lon data using a geographic CRS which is not specified in the data file. If, as you suggest, Mapper were to prompt for the CRS for live GPS data, Mapper would naturally retain that CRS in its template metadata -- it already has a field for that purpose. Perhaps that's exactly what you had in mind.

It is.

In my opinion, it would be tedious as a user to select a CRS every time I start to record a track. Most days my map and my GPS receiver do not change, so the appropriate CRS doesn't change either.

I agree. However, it must be possible to change it. Use a reasonable default (I suggest current SBAS* GPS) and allow override in config. Perhaps a popup alert when GPS location is tuned on in order to make the user aware of the setting.

A geographic CRS is not flat and thus is not suited to a matrix type of template transform. Mapper potentially manages a variety of projected CRSes, as well as actively transforming GPX tracks based on their geographic CRS.

When you go through the georeference setup, you select a projected CRS. Any other CRS, projected or geographic, can be converted (by PROJ) to the one you selected. If you are concerned about performance, you can use PROJ to calculate parameters for a simplified transformation (2x2 shift/rotation matrix, 2x2 plus linear position, etc.)

• There is a small difference between SBAS corrected and autonomous GPS positions, but the difference is much smaller than what we can get from autonomous GPS.

[sfroyen]

Actually I would like to simply take the geographic CRS the projected CRS is based on. And use a direct pipeline between CRS everywhere else (#1598).

That seems unnecessarily complicated. Consider the data we encounter. Here are some examples that I've seen.

More common:

1. DXF files from LIDAR data, NAD83(2011) epoch 2010.00 UTM zone 13N / EPSG:6342 (x',y',h)
2. Live SBAS GPS signals, WGS84 G1762 / EPSG: 7664 (x,y,z), 7665 (lat,lon,h), 9057 2D (lat,lon)
3. Pictometry (and similar) areal images, EPSG:4326 (I have to shift these in current Mapper)
4. Google satellite images, Web/Pseudo Mercator / EPSG:3857 (I use QGIS to capture these and convert them to projected coordinates.
5. Open Street Map (same as Google)
6. Recent GPS tracks / GPX files, most likely using WGS84 G1762, user should verify
7. OCAD files, ?

Less common:
6) Leica Smartnet RTK corrected GPS signals, NAD83(NA2011) Epoch 2010.00 / EPSG:6318 (lat,lon,h)?
7) Older GPX files
8) Older georeferenced images

As you can see, almost all the examples use (or approximate) the latest incarnation of WGS84 (G1762). Furthermore it appears that EPSG data for the later realisations are more complete.

Here are results of using PROJ's cs2cs to transform between different realisations of WGS84.

Input coordinates are 40N 105W 0
WGS84 G1762 -> G1762 (2014): 40.00000000 -105.00000000 0.00000000 (identity transformation)
WGS84 G1762 -> G1674 (2013): 39.99999987 -104.99999988 0.04277009
WGS84 G1762 -> G1150 (2002): 39.99999984 -104.99999992 0.01831625
WGS84 G1762 -> G873 (1997): 40.00000000 -105.00000000 0.00000000 (error)
WGS84 G1762 -> G730 (1994): 40.00000000 -105.00000000 0.00000000 (error)
WGS84 EPSG:4326 -> G1762 (2014): 40.00000000 -105.00000000 0.00000000
WGS84 EPSG:4326 -> G1674 (2013): 40.00000000 -105.00000000 0.00000000
WGS84 EPSG:4326 -> G1150 (2002): 40.00000000 -105.00000000 0.00000000
WGS84 EPSG:4326 -> G873 (1997): 40.00000000 -105.00000000 0.00000000
WGS84 EPSG:4326 -> G730 (1994): 40.00000000 -105.00000000 0.00000000

Note that cs2cs produces zero shifts when transforming to WGS84 versions before 2002. It might be specific to cs2cs or it might be missing parameters in the EPSG data files. Also, transformations between EPSG:4326 and any of the others produce no shift.

[pkturner]

Analysis and the testing described in #1757 shows that the track positioning as changed with PROJ 7.1.0 is considerably more accurate than what Mapper was getting with PROJ 4.9.3.

[pkturner]

The tests developed in #1764 check for the specific discrepancies discussed above, namely

• changes in positioning of GPX tracks relative to previous versions of Mapper
• discrepancy between positioning of GPX tracks depending on the setting of GDAL loading of GPX tracks

[sfroyen]

This issue keeps bugging me with v0.9.5. The gpx track file that I get from using OOM on a tablet is still shifted when I use it as template in the desktop version (MacOS if that is relevant). The culprit is (most likely) the default EPSG:4326 which ignores all changes between the original WGS84 and the current realization, used by GNSS coordinates today. Could we please stop using the this default and instead explicitly initialize with the current GPS CRS, EPSG:9753 / WGS84 (G2139)? This should match today's GNSS coordinates exactly and most older geographic data within a few centimeters. Almost all users will receive this type of data. If GDAL or PROJ is missing support for this latest realization, use EPSG:7664 / WGS 84 (G1762) instead.

If someone could build OOM with this default, I would be happy to test the result :-)