Inaccurate import of images and tracks that use geographic coordinates

[sfroyen]

Steps to reproduce

1. Create a map using an accurate CRS for your area. (I used EPSG:6342 which is the current NAD83 datum for Colorado).
2. Go out with an accurate GNSS device and add some point and line features that are easily seen in images, using the Android version of OOM.
3. Copy the files from the Android device to the your non-Android OOM.
4. Import the GPX track data from the Android session. Notice that the tracks are offset compared to point and line features that were added in step 2. In my case the offset is just over 1 meter.
5. Use QGIS to capture a georeferenced Google areal image of the area that was used in step 2.
6. Add the image as a georeferenced template in OOM. Notice that the image is offset compared to point and line features that were added in step 2.

Actual behaviour

Any georeferenced data that use geographic coordinates is offset when added as a template.

Expected behaviour

Georeferenced data using geographic coordinates are not offset when added as templates.

Root cause

It appears that OOM is not specifying a default EPSG for geographic coordinates. Without a default, GDAL/PROJ uses EPSG:4326. This is an ensemble EPSG that assumes all geographic coordinate datums are the same and ignores the shifts that have accumulated over time. I have verified that this is the cause of the shifts by using QGIS to convert the images and tracks from EPSG:7667 (the 2019 datum revision) to EPSG:6342. When importing the converted files as templates, the shifts are gone.

Fix

Hardcoding EPSG:7667 in the initialization of GDAL/PROJ would almost all of the cases. Adding a configuration setting for the EPSG code would solve any issue with (very) old data (pre 2012).

PS There are several other issue reports that discuss the same issue. I filed this issue because I have a simple solution that will fix most of the other cases.

Configuration

Mapper Version: 0.95
Operating System: MacOS Sonoma

[pkturner]

I agree that EPSG:4326 is poor for geographic coordinates, and that something like EPSG:7667 is an improvement.

It's been a long time since when I looked at this issue. Initializing PROJ differently did not appear to be an option.

[krticka]

I think that accurate GNSS receivers and data will be more and more used, so it seems to me logical to stop using low precision ensemble which is based on dynamic datums. For example if you want to save any data to 4326, QGIS now warns you about limited accuracy. Thinking about this, not sure what is the best solution. What about to add parameter in Mapper settings where you can define WGS84 epoch for importing data?
BTW latest WGS 84 revision should be from 2021 (G2139).
BTW2 importing/opening data from different CRS Mapper is like a black box for me as a user. You never know what transformation was used, at least some information about this would be useful.

[krticka]

Just for the record
EPSG 5514
purple buildings - added as a template from shapefile
grey buildings - added as a template from the same shapefile, but with prj file deleted. Imported using real coordinates

[sfroyen]

pkturner - If your tests were before ~2019, EPSG files did not include shifts. So your tests with GDAL/PROJ would be as you describe (I had the same experience). Now, most of the newer ones do.

[pkturner]

@sfroyen
Your #1264 got me started experimenting with the need for shifts back in 2020. Using EPSG:7665 for GPX tracks works nicely for PROJ. With GDAL it's problematic. In a discussion with @dg0yt, he and I considered a method of forcing EPSG:7665 on GDAL, and he had a reason for not doing that.

Hardcoding EPSG:7667 in the initialization of GDAL/PROJ would almost all of the cases.

I think you mean EPSG:7665.

For PROJ, I know how EPSG:7665 could be substituted for +datum=wgs84. If I knew a good way to make the same substitution when setting up GDAL, I would create a PR.

[sfroyen]

@pkturner
You're correct, EPSG:7667 is incorrect. Not sure where that came from.

Here's a table of EPSGs for WGS84:

1994 G730 EPSG:7656 (x,y,z), EPSG:7657 (lat,lon,h), EPSG:9053 2D (lat,lon)
1997 G873 EPSG:7658 (x,y,z), EPSG:7659 (lat,lon,h), EPSG:9054 2D (lat,lon)
2002 G1150 EPSG:7660 (x,y,z), EPSG:7661 (lat,lon,h), EPSG:9055 2D (lat,lon)
2012 G1674 EPSG:7662 (x,y,z), EPSG:7663 (lat,lon,h), EPSG:9056 2D (lat,lon)
2013 G1762 EPSG:7664 (x,y,z), EPSG:7665 (lat,lon,h), EPSG:9057 2D (lat,lon)
2019 G1762' EPSG:7664 (x,y,z), EPSG:7665 (lat,lon,h), EPSG:9057 2D (lat,lon) ???
2021 G2139 EPSG:9753 (x,y,z), EPSG:9754 (lat,lon,h), EPSG:9755 2D (lat,lon)

For my location we currently* use NAD83(2011) UTM 13N / EPSG:6342
*in process of being replaced

Running the shell script below (location ~40N 105W) produces a table of EPSG:6342 coordinates for various WGS84 realizations.

for i in 4326 7657 7659 7661 7663 7665 9754
do
echo EPSG:$i to EPSG:6342; echo 40.02020698N 105W | cs2cs EPSG:$i EPSG:6342
done

EPSG:4326 to EPSG:6342
500000.00 4430000.00 0.00
EPSG:7657 to EPSG:6342
500000.00 4430000.00 0.00
EPSG:7659 to EPSG:6342
500000.00 4430000.00 0.00
EPSG:7661 to EPSG:6342
500000.00 4430000.00 0.00
EPSG:7663 to EPSG:6342
500000.86 4429999.28 0.89
EPSG:7665 to EPSG:6342
500000.86 4429999.28 0.89
EPSG:9754 to EPSG:6342
500001.08 4429999.37 0.87

Note that EPSG:4326 as well as all WGS84 realizations before 2012 (EPSG:7663) are missing the shifts.

[pkturner]

PROJ and GDAL now support combining epoch with CRS.

$ echo  40.02020698 -105 0.0 | cs2cs EPSG:7665 +to EPSG:6342
500000.86       4429999.28 0.89
$ echo  40.02020698 -105 0.0 | cs2cs "WGS 84 (G1762)" +to EPSG:6342
500000.86       4429999.28 0.89
$ echo  40.02020698 -105 0.0 | cs2cs "WGS 84 (G1762)@1997" +to EPSG:6342
500000.86       4429999.28 0.89
$ echo  40.02020698 -105 0.0 | cs2cs "WGS 84 (G2139)@1997" +to EPSG:6342
500000.87       4429999.28 0.88

Notice that WGS 84(G1762) is synonymous with EPSG:7665, and PROJ implicitly assigns an epoch of 1997.

$ echo  40.02020698 -105 0.0 | cs2cs EPSG:9754 +to EPSG:6342
500001.08       4429999.37 0.87
$ echo  40.02020698 -105 0.0 | cs2cs "WGS 84 (G2139)" +to EPSG:6342
500001.08       4429999.37 0.87
$ echo  40.02020698 -105 0.0 | cs2cs "WGS 84 (G1762)@2010" +to EPSG:6342
500001.08       4429999.37 0.87
$ echo  40.02020698 -105 0.0 | cs2cs "WGS 84 (G2139)@2010" +to EPSG:6342
500001.08       4429999.37 0.87

Notice that WGS 84 (G2139) is synonymous with EPSG:9754, and PROJ implicitly assigns an epoch of 2010.

The point of these examples is that for the recent realizations of WGS 84 that PROJ handles dynamically, the calculation of shift is dependent on the epoch much more than the specific realization. The epoch should be determined by the date of the observation, e.g. by the date when a track was recorded.

And the relationship has shifted another 23 cm since the epoch assumed by PROJ for EPSG:9754.

$ echo  40.02020698 -105 0.0 | cs2cs "WGS 84 (G2139)@2024" +to EPSG:6342
500001.30       4429999.46 0.85

The fix suggested by @sfroyen, if it can be implemented, would beg for Mapper to also support the epoch of observations for tracks and for imagery that is georeferenced to WGS 84.

[pkturner]

By the way, OpenStreetMap uses WGS 84 latitude and longitude for its data. It is disturbing to think of all of OSM's data on North America slowly becoming erroneous as the continental plate drifts away.

[sfroyen]

@pkturner, I believe that high accuracy archived GNSS data is usually recorded with positions transformed back to the frame epoch. For the coordinates of a location today one therefore has to transform the archived coordinates using a Helmert transformation with the time elapsed since the epoch.

For autonomous GPS measurements, the positions are provided in present time realization of WGS84. For archival accuracy these should be transformed back to the frame epoch. When using SBAS corrections, the positions would use the realization of SBAS satellite. (I read somewhere that WAAS satellites are updated every 6 months.) Similarly, for RTK corrected data, the realization is whatever the RTK provider uses.

For orienteering mapping, the differences caused by time are not important. For updates of older maps, the map's georeferencing could be adjusted to match the current realization (keeping map coordinates as is and shifting the geographic and reference coordinates to match).

My original post was meant as simple fix for the incorrect shifts of images and gps tracks. These tend to be recent data with coordinates very close to the current WGS84 realization, so using G2139 (EPSG:9754) should be optimal in almost every case. And for other cases, it should be no worse than the default (EPSG:4326). If needed, we could add a UI for CRS selection later.

[sfroyen]

The workaround steps using QGIS are as follows (I'm hoping this description might suggest a OOM fix):

1. Open QGIS (create a new project if needed)
2. Select Web->QuickMapServices, search for Google and use, e.g., Google Satellite Hybrid
3. Select an image area
4. Choose EPSG:7665 (WGS84 G1762) as the CRS (my QGIS does not have the latest realization)
5. Export area as image (Project->Import/export->Image)
6. Open OOM and load map. Navigate to image location
7. Open the image as a template (georeferenced). Choose "by EPSG code" and select EPSG 6775
8. Verify image is correctly placed correctly
9. Repeat steps 4, 5 and 7 using your map EPSG in steps 4 and 7 (I used 6342).
10. Verify image is shifted.

Image step 7:

Image step 10:

The purple point symbols were added to the map using Android OOM and a high accuracy GPS (using RTK with current WGS84 realization as the correction CRS). The shifts are ~1m. The location is an unused runway, northern end at the US Air Force Academy.

[sfroyen]

Following up my previous comment...

EPSG:3857 - WGS 84 / Pseudo-Mercator - Spherical Mercator - used by Google Maps, OpenStreetMap, Bing, ArcGIS, ESRI, works just as well as EPSG:7665 (see step 7 above). As long as you specify the CRS of the template image.

It seems that OOM is working correctly (at least to ~20 cm) and that the only issue is with the QGIS conversion to non-geographic coordinates.

[sfroyen]

... and, of course, the gpx track templates

They are positioned incorrectly both when added as templates and when imported as map objects. Could those cases be handled like the image templates, with an option to select a CRS and a default to use EPSG:7665 (or later)?

[pkturner]

@sfroyen said

It seems that OOM is working correctly (at least to ~20 cm) and that the only issue is with the QGIS conversion to non-geographic coordinates.

... and, of course, the gpx track templates

I'm only now fully analyzing your report. Originally I accepted your prescription of of hardcoding EPSG:7667 into Mapper almost everywhere. I even put an effort into implementing that fix (#2205). Then I had to look again, because the fixed Mapper places GPX tracks in the same position as v0.9.5 (as long as GDAL is disabled).

The shift that your images show -- I agree it is in the QGIS creation of the image in EPSG:6342. I have analyzed the shift which @dg0yt encountered in #2008 (NAD83 changes in PROJ 8.2.0), and the two appear the same. What I think happened is

1. Your Android Mapper is v0.9.5. This release was built with an earlier PROJ, that does not have the # 2008 bug.
2. Your desktop Mapper is equivalent to the OBS build openorienteering-mapper-unstable_0.9.20211224-0+544.1. This build had the # 2008 problem.
3. Your QGIS is affected by the same flaw in PROJ 8.

I suggest you check how Mapper v0.9.5 handles GPX tracks on your EPSG:6342 maps. If it's not shifted, then this report is a duplicate of # 2008.

[pkturner]

@krticka
Your example images appear to show a shift of more than 5 meters. Using a better WGS84 realization, or specifying the epoch, would only change the alignment by 1 or 2 meters, so that's not the main problem.

Your comment about "You never know what transformation is being used" is an interesting point. PROJ has an ability to dump the transformation, at least on the command line. As a feature in Mapper I expect it would rarely be used, but for your case it would be worthwhile.