[gdal-dev] gdal-dev Digest, Vol 243, Issue 19

Fox, Shawn D (US) shawn.fox at baesystems.us
Mon Aug 26 09:03:15 PDT 2024


----------------------------------------------------------------------

Message: 1
Date: Sat, 24 Aug 2024 08:15:32 -0400
From: Greg Troxel <gdt at lexort.com>
To: "Fox, Shawn D (US) via gdal-dev" <gdal-dev at lists.osgeo.org>
Subject: Re: [gdal-dev] Upgrading to GDAL 3.4.2 from 1.10.0 shows
	different outputs from gdaltransform and corresponding C++ API
Message-ID: <rmiplpy86cb.fsf at s1.lexort.com>
Content-Type: text/plain

"Fox, Shawn D (US) via gdal-dev" <gdal-dev at lists.osgeo.org> writes:

> I have noticed that from GDAL 1.10 - GDAL 3.1.0 that gdaltransform is 
> producing different results for NAD27 GCS to WGS84 GCS 
> transformations.  The differences in longitude in the output is more 
> significant.  Having tested the differences in feet between the points 
> output from each version I can see a difference of ~5 feet primarily 
> due to the change in the longitude value.

> Would anyone look at the following results and share your thoughts 
> about why the result of this transformation is different in a newer 
> version?  Could this be related to the concept of WGS 84 realization 
> being different in newer versions of GDAL?  I am also linking GDAL
> 3.4.2 with PROJ 8.1.1 at this time.  These results were produced with 
> GDAL 3.4.2, but the same results occur with GDAL 3.1.0 as well.  GDAL
> 3.4.2 is what I am currently trying to upgrade to.  The results from 
> the C++ programming API are similar in difference.  I just used 
> gdaltransform for the purposes of this message because it was a more 
> convenient way of producing something that others might be able to 
> reproduce.

You didn't explain

  - When you were using gdal 1.10 (what year was that??), what version
    of proj were you using?
Sfox: Proj was version 4.4.9, but why does it matter what year I was using GDAL 1.10?  We are using it now in 2024 but it is the same library that was used in 2017, 2018, etc.  The library doesn't produce different results in different years.  I've only seen different results when changing the version.

  - What happens when you just use proj, and not gdal?  It's always best
    to simplify the problem.  Include the results of "projinfo".
Sfox: I haven't tried that because I assume that GDAL is simply using proj to do these transformations, but I will investigate that and try it out. The problem is that we have the GDAL API to work with so I'm not sure if it matters what the proj APIs directly return to us.

  - What you get when using NGS's tools, which are more or less the gold
    standard for the right answer.

Sfox: What does NGS stand far and where can I find those tools?

  - What your approach is to "WGS84", which is an ensemble.  It contains
    members that differ by about 2m, including one that more or less
    matches NAD83(1986) and ones that very closely match recent ITRF.
Sfox: What is an ensemble?  The GDAL API requires me to set the source and target SRS and provide the input coordinates but that is it.
Below is some example source code so that you can see the native C++ APIs being used.  

      /* if source is GCS NAD 27, then OSRSetWellKnownGeogCS */
      errorSource = OSRSetWellKnownGeogCS( sourceSRS, "NAD27");
      /* if target is GCS WGS 84, then OSRSetWellKnownGeogCS */
      errorTarget = OSRSetWellKnownGeogCS( targetSRS, "WGS84");
      
      /* create OGRCoordinateTransformationH object */
      /* returns the obj on success, NULL on failure */
      OGRCoordinateTransformationH ptrCT =
        OCTNewCoordinateTransformation( sourceSRS, targetSRS );

      /* check to make sure that it exists */
      if (NULL != ptrCT)
      {
        /* perform the transform */
        /* returns TRUE on success, FALSE on failure of some or all points */
        status = OCTTransform( ptrCT, NUMBER_OF_POINTS, x_lon, y_lat, NULL );
      }

  - What you are thinking about plate motion -- and if that even makes
    sense given the low accuracy of NAD27.   And after you explain that,
    why you are using WGS84 rather than NAD83(2011) epoch 2010.0.
Sfox: I wasn't thinking about anything related to that.  The results of the transformation haven't changed between GDAL 3.1.0 and GDAL 3.4.2 both using PROJ 8.1.1.  My question is simply why would  the calculations change since changing the PROJ version used with GDAL?  If it is as simple as the default realization has changed then that would be useful information. 

> GDAL 1.10.0 results
> ./gdaltransform -s_srs NAD27 -t_srs WGS84
>
> -117.258911 32.804191
> -117.25978144845 32.8042393819884 1.32247805595398e-07
>
> GDAL 3.4.2 results
> ./gdaltransform -s_srs NAD27 -t_srs WGS84
> -117.258911 32.804191
> -117.259764476223 32.8042388409449 0

If you are seeing less than 2m difference it's really hard to say either is wrong, and it's not really a bug that it changed.

  NAD27 is a low-accuracy datum, with internal distortions.

  Generally, there are grid shift files to use with NAD27-NAD83
  conversions because of this.  They are not included in the main proj
  tarball, and are available in various data tarballs and via the net.
  Having them and not having them is quite different.

  proj treats NAD83 to WGS84 as a null transform because WGS84 has a
  member that matches NAD83.  But if you care about 1m, you have to:

    stop using the WGS84 ensemble

    once explicit about which member you mean, perhaps WGS84(G2139), pay
    attention to NAD83 to ITRF transforms.

Sfox: I have a requirement to use WGS84 which is why it is being used.  What do you mean by the question, "which member you mean"?  The API in GDAL 3.4.2 does not allow me to select a specific WGS84 realization yet.  So how would I know if GDAL is using G2139?  I have not seen any documentation within the GDAL docs that tell me that.  You are using terminology such as "member" and "ensemble" which I am unfamiliar with.  
  
> I have read that the support for specifying epochs related to realizations is not available until GDAL 3.8.0.

And, you aren't specifying them anyway :-)  Your problems are far greater than epochs, which are only going to make sense once you resolve the ensemble issues.

Sfox: Of course not, because GDAL doesn't allow me to specify them.  What do you mean by ensemble vs realization?  You mentioned G1239 which is a realization and yet refer to it as an ensemble.

And, it's nearly impossible to believe that any data which is originally in NAD27 is going to have errors that are as low as 2m in the first place.

Sfox: I don't see the relevance.  Unless the code changed within GDAL and/or PROJ, there is no reason that the exact same transform with the exact same input should produce a different output. I see what you mean that if you start with something that isn't accurate to 2m then a calculation from that isn't going to be accurate either, but the fact remains that from one version of GDAL or PROJ to another the results changed.  It would be nice to know why so that we can be comfortable moving forward with an upgrade to a newer version of GDAL and PROJ.

> Here are some references I've read so far.
> Transforming between WGS84 Realizations | Journal of Surveying 
> Engineering | Vol 148, No 2 
> (ascelibrary.org)<https://ascelibrary.org/doi/10.1061/%28ASCE%29SU.194
> 3-5428.0000389>

That's paywalled so not useful to th rest of us :-(

But the point is that it makes senes to transform between data that is actually in different members of the ensemble, because it was gathered in those members.

It does not make sense to transform *to* this ensemble, unless you declare loudly and then act consistently with a belief that 2m of difference is effectively zero -- in which case you have not presented a problem.


------------------------------

Subject: Digest Footer

_______________________________________________
gdal-dev mailing list
gdal-dev at lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev


------------------------------

End of gdal-dev Digest, Vol 243, Issue 19
*****************************************


More information about the gdal-dev mailing list