Re: Different Projection Results in AWS Managed Services

Pat Blair <[email protected]> writes:

> Below I am pasting some SQL that I hope will illustrate what we are
> seeing.  When we transform a coordinate from EPSG:26860 (NAD83 /
> Minnesota Central (ftUS)) to WGS84 in two different environments (“AWS
> RDS” and “AWS Aurora”) we get a different result with a displacement
> of about 0.12 meters.

Multiple things going on here.  You need two things to start:

  Geodetic clarity on what you are doing.  You are calling 0.12m a
  discrepancy which means the allowable fuzz in understanding is zero.
  And, unless you are working unusually closely with the US DoD, you
  can't possibly be obtaining positions in WGS84 to anywhere near the
  0.1 m level, because there are no mechanisms for accurate access to
  that datum.  (You'd need carrier phase data from the official
  stations, which as far as I know they don't publish.)

  The ability to compare test vectors to authoritative sources such as
  NGS.

> The proj4 definitions for the projected
> coordinate system look to be the same on both systems, though since
> these are managed database services we don’t have direct access to any
> of the projection databases, grids, etc.

Hence, these questions need to be should be directed to your hosted
service provider.  But once you achieve clarity, you may not need to.

A key point is whether proj has access to and is using downloadable
grids.  You can't tell from the info string if they are present already.

Note you are using 8.0.1 and 9.1.0.  Install both locally and see if you
can reproduce the differences.  This may have nothing to do with hosted
services.

> Many thanks in advance if you can provide some advice!  We would
> certainly be happy to provide any more information that might be
> helpful in describing what we are seeing.

So....

  EPSG:26860 that you wrote about above is "NAD83(HARN) / Nebraska
  (ftUS)".  But it seems you mean EPSG:26850.  This just says NAD83,
  which strikes me as surprising today.
      https://epsg.org/crs_4269/NAD83.html
  You are concerned about 0.12 meter, but it's hard to believe you have
  coordinates in the original NAD83 datum that are better than 1 m.
  It's likely there is a code for your SPCS with NAD83(2011) epoch
  2010.0 instead.  But first the question is what frame are your data in
  really.  I have found that data from a GIS environment often lacks
  clarity about these fine points -- doesn't have meter-level accuracy
  anyway.  (With MassGIS some data is sub-meter, but it, perhaps not
  surprising, is recent and is labeled more carefully.  An example is
  15cm imagery which has been ground-controlled to NAD83(2011) epoch
  2010.0).

  WGS84 is an ensemble.  There are a number of members, and
  unfortunately one of them is the original, WGS84(TRANSIT), about 2m
  different than modern members.  When you use an ensemble you mean that
  **your data is in one of the members but you don't know which**.  In
  my view, transforming **to** an ensemble is non-sensical and I
  recommend that you come up with a new plan where you don't do that.
  It does not make sense to transform to WGS84 and consider any
  discrepancy under 2m any kind of error.  Generally people want WGS84
  because some standard refers to it, and my guess at what the standard
  authors would have wanted had they really understood is the most
  recent realization.  Either that or they do not believe their standard
  will be used with data that is accurate to better than a few meters.
  (Which you might or might not be doing :-)

  Because NAD83 is an ensemble and particularly because WGS84 is an
  ensemble, and because of how proj currently works (which I think is
  suboptimal but have not offered a patch to address), proj will use a
  null transform from NAD83 to WGS84.  If you expect your coordinates to
  be in WGS84(G2159), that's wrong.  You can instead transform to some
  ITRF, or perhaps now you can transform to a specific WGS84.  But you
  have to run proj to look at the transforms and understand what ought
  to be.

  WGS84 realizations and ITRF realizations (that you might use as WGS84
  proxies) are dynamic datums.  Points have velocities (2-3 cm/y in NA),
  and at the 0.12m level this cannot be neglected.  Coordinates in WGS84
  realizations and ITRF realizations are typically expressed at some
  reference epoch, or a specific epoch.  That's hard to start with, and
  the reference epoch of an ensemble puts you through the thin ice, not
  just on it.

> /**
> AWS RDS
> */

> -- POINT(-95.1331380928918 45.773407072930254)

FWIW, running your transform on my system (NetBSD 10 amd64, proj 9.5.1,
postgis 3.5.3), I got the same result as RDS.  But, here I ask to
transform to NAD83(1986), to NAD83(2011) epoch 2010.0 (3d, but that's
ignored I think), to WGS84(ensemble) as you did, and to ITRF2014.

First, let's transform from the CRS you said the data is in to the
modern (but still ftUS :) flavor with NAD83(2011) epoch 2010.0:

gdt=> SELECT ST_AsText(ST_Transform(ST_SetSRID('POINT(2399319.7374562174 611365.0841263086)'::geometry, 26850), 6501)) as crs6501;
                   crs6501                   
---------------------------------------------
 POINT(2399319.3357411437 611365.1332752574)
(1 row)

Just eyeballing, that's 45 cm right there.


Let's transform to NAD83(1986), NAD83(2011) epoch 2010.0,
WGS84(ensemble) and ITF2014.

gdt=> SELECT ST_AsText(ST_Transform(ST_SetSRID('POINT(2399319.7374562174 611365.0841263086)'::geometry, 26850), 4269)) as crs4269;
                   crs4269                   
---------------------------------------------
 POINT(-95.13313652618976 45.77340681030725)
(1 row)

gdt=> SELECT ST_AsText(ST_Transform(ST_SetSRID('POINT(2399319.7374562174 611365.0841263086)'::geometry, 26850), 6319)) as crs6319;
                   crs6319                   
---------------------------------------------
 POINT(-95.13313810253163 45.77340693280185)
(1 row)

gdt=> SELECT ST_AsText(ST_Transform(ST_SetSRID('POINT(2399319.7374562174 611365.0841263086)'::geometry, 26850), 4326)) as crs4326;
                   crs4326                   
---------------------------------------------
 POINT(-95.1331380928918 45.773407072930254)
(1 row)

gdt=> SELECT ST_AsText(ST_Transform(ST_SetSRID('POINT(2399319.7374562174 611365.0841263086)'::geometry, 26850), 7912)) as crs7912;
                   crs7912                    
----------------------------------------------
 POINT(-95.13267275410031 45.773464527082645)
(1 row)


Now the same, but relabeling the source as NAD83(2011) epoch 2010.0:

gdt=> SELECT ST_AsText(ST_Transform(ST_SetSRID('POINT(2399319.7374562174 611365.0841263086)'::geometry, 6501), 4269)) as crs4269;
                   crs4269                    
----------------------------------------------
 POINT(-95.13313494984804 45.773406687812546)
(1 row)

gdt=> SELECT ST_AsText(ST_Transform(ST_SetSRID('POINT(2399319.7374562174 611365.0841263086)'::geometry, 6501), 6319)) as crs6319;
                   crs6319                   
---------------------------------------------
 POINT(-95.13313652618976 45.77340681030725)
(1 row)

gdt=> SELECT ST_AsText(ST_Transform(ST_SetSRID('POINT(2399319.7374562174 611365.0841263086)'::geometry, 6501), 6319)) as crs4326;
                   crs4326                   
---------------------------------------------
 POINT(-95.13313652618976 45.77340681030725)
(1 row)

gdt=> SELECT ST_AsText(ST_Transform(ST_SetSRID('POINT(2399319.7374562174 611365.0841263086)'::geometry, 6501), 7912)) as crs7912;
                  crs7912                   
--------------------------------------------
 POINT(-95.13267117775811 45.7734644045764)
(1 row)


Use NCAT or other NGS tools to transform your point to NAD83(2011) epoch
2010.0. Also do it in proj.  You might need to use proj to transform to
the base CRS of your SPC, EPSG:4269.  That's just a projection, and not
a datum shift.  You might see how different "NAD83(1986)" is from
"NAD83(2011) epoch 2010.0".

  https://www.ngs.noaa.gov/NCAT/

Then, you can use various tools to transform that to WGS84(G2296) or
WGS84(G2139).  NGS doesn't deal as much in WGS84.  As I understand it
G2296 is equivalent to ITRF2020 and G2139 is equivalent to ITRF2014.
This is probably useful

  https://www.ngs.noaa.gov/cgi-bin/HTDP/htdp.prl?f1=4&f2=1

This is old, but it gives you some sense  of the difficulties.
  https://www.ngs.noaa.gov/CORS/Articles/WGS84NAD83.pdf


> /**
> AWS Aurora.
> */

> -- POINT(-95.13313652618976 45.77340681030725)

This matches my transform to NAD83(1986), which is a projection and not
a datum shift.

proj 9 had major improvements in US national datum transforms, and seems
to be transforming between NAD83(1986) and NAD83(2011) epoch 2010.0.  My
impression is that it is close to what NCAT produces.  Such transforms
are almost certainly gridded as they are really about distortions in the
early frames from errors, greatly reduced later.

What's surprising me and I need to dig in more is that it appears the
transformation from 6319 to 4326 is not null, but it's also quite
different from 7912 (ITRF2014).

This may raise more questions than answers, but you really need to drive
definitional uncertainty to near zero if you want to talk about 12 cm.

Please report back with a summary of your understanding after
investigating, as I suspect it is helpful to others, as well as in
understanding the current state of proj.  Dealing with ensembles and
with older datums is hard.

Greg