[MapProxy] shift when reprojecting WMS

Sat Jun 18 05:36:28 PDT 2022

Travis Kirstine <traviskirstine at gmail.com> writes:

> These are good questions. Hopefully my response will answer your questions.
>
> In my old install of pyproj v2.2.2 the transformation from EPSG:3857 to UTM
> Z17N NAD83 Original (EPSG:26917) and UTMZ17N CSRS (EPSG:2958) would produce
> the same coordinate results.
>
> # python2.7 pyproj v2.2.2
> from pyproj import CRS, Transformer
> crs_3857 = CRS("EPSG:3857")
> crs_26917 = CRS("EPSG:26917")
> crs_2958 = CRS("EPSG:2958")
> wmerc2utmOrg =  Transformer.from_crs(crs_3857, crs_26917, always_xy=True)
> wmerc2utmCSRS =  Transformer.from_crs(crs_3857, crs_2958, always_xy=True)
> wmerc2utmOrg.transform(-8851000, 5433000)
> (619889.7678080541, 4849626.937111458)
> wmerc2utmCSRS.transform(-8851000, 5433000)
> (619889.7678080541, 4849626.937111458)
>
> In this case the web mercator to utm nad 83 original is correct, however

See [1] below.

> the CSRS result is off as both transformations are treated the same.  I can
> determine a more accurate position for the CSRS coordinates by applying a
> gridshift using the proj cs2cs utility:
>
> cs2cs +proj=utm +zone=17 +lat_0=0 +lon_0=-81 +k=0.9996 +x_0=500000 +y_0=0
> +ellps=GRS80 +units=m +nadgrids=/home/grids/ON83CSv1.gsb +to
> +init=EPSG:2958 <<EOF
> 619889.76 4849626.93
> EOF
> 619889.73       4849626.65 0.00
>
> This is a shift of -0.03m X and 0.28m Y from utm nad83 original to utm
> csrs.

When you have data in web mercator, what realization of WGS84 do you
think that data is in?  The point is that WGS84(TRANSIT) and
WGS84(G2159) are both members of the ensemble, yet are very different,
far more so that the differences you are talking about.  The error
estimate from the ensemble is large and causes proj to decide to use
null transformations.

[1] If you are talking about sub-meter, then talking about a correct web
mercator to NAD83(1986) transform requires being clear about which WGS84
realization you are using.

This is hard because
  - EPSG does not have codes for web mercator on the various WGS84
    realizations.
  - The web mercator world does not have codepoints to specify the
    realization in use.
  - The web mercator spec, as far as I have found, is silent/fuzzy on
    this entire issue.

WGS84(TRANSIT) is equal to NAD83(1986), more or less by construction.
So correct would be just the projection with no shift.  The shift from
WGS84(TRANSIT) to WGS84(G2159) is huge compared to the shifts between
any two NAD83 realizations.

Your grid shifts from NAD83(1986) to NAD83(CSRS) are small, and that
seems like the usual improved datum with reduced internal distortion.

So, converting from the WGS84 ensemble to any particular NAD83 with a
null transform is defensible (to the extent null transforms of ensembles
are ok in the first place).

Given that you are saying "more accurate" I am guessing you have
accuracy information about your data in "web mercator" that is outside
the definition of 3857, and thus you don't really have 3857.

> After I upgraded to pyproj v3.0.1 I'm seeing a different result for the Web
> Mercator to UTM CSRS transform
>
> # python.3.6.8 pyproj v3.0.1
> from pyproj import CRS, Transformer
> crs_3857 = CRS("EPSG:3857")
> crs_26917 = CRS("EPSG:26917")
> crs_2958 = CRS("EPSG:2958")
> wmerc2utmOrg =  Transformer.from_crs(crs_3857, crs_26917, always_xy=True)
> wmerc2utmCSRS =  Transformer.from_crs(crs_3857, crs_2958, always_xy=True)
> wmerc2utmOrg.transform(-8851000, 5433000)
> (619889.7678080541, 4849626.937111458)
> wmerc2utmCSRS.transform(-8851000, 5433000)
> (619890.0493627792, 4849625.986608229)
>
> So I'm seeing a shift of -0.28m X and 0.95m Y from utm nad83 original to
> utm csrs which is off by -0.25m X and 0.67m Y from my expected result
>
> My guess this is due to how pyproj utilizing a helmert transformation
>
> # python.3.6.8 pyproj v3.0.1
> from pyproj.transformer import TransformerGroup
> wmerc2utmOrg_tg = TransformerGroup("epsg:3857", "epsg:26917")
> wmerc2utmOrg_tg.transformers[0].description
> 'Inverse of Popular Visualisation Pseudo-Mercator + Inverse of NAD83 to WGS
> 84 (1) + UTM zone 17N'
> wmerc2utmOrg_tg.transformers[0].definition
> 'proj=pipeline step inv proj=webmerc lat_0=0 lon_0=0 x_0=0 y_0=0
> ellps=WGS84 step proj=utm zone=17 ellps=GRS80'
> wmerc2utmCSRS_tg = TransformerGroup("epsg:3857", "epsg:2958")
> wmerc2utmCSRS_tg.transformers[0].description
> 'Inverse of Popular Visualisation Pseudo-Mercator + Inverse of NAD83(CSRS)
> to WGS 84 (2) + UTM zone 17N'
> wmerc2utmCSRS_tg.transformers[0].definition
> 'proj=pipeline step inv proj=webmerc lat_0=0 lon_0=0 x_0=0 y_0=0
> ellps=WGS84 step proj=push v_3 step proj=cart ellps=WGS84 step inv
> proj=helmert x=-0.991 y=1.9072 z=0.5129 rx=-0.0257899075194932
> ry=-0.0096500989602704 rz=-0.0116599432323421 s=0
> convention=coordinate_frame step inv proj=cart ellps=GRS80 step proj=pop
> v_3 step proj=utm zone=17 ellps=GRS80'
>
> This should probably be posted to the proj community but any help would be
> great.

It is a proj question, and you have reduced your mapproxy case to a
small proj case.

I am far from a proj expert, but I think the issue is that you asked for
a transform between two CRSes, and did not specify a geographic area
that would allow proj to use the grid shift file.   I don't know how to
do that.

But I think your real problem is likely that all your coordinates are a
meter off, or perhaps your imagery is too.

You might find this interesting:
  https://webapp.csrs-scrs.nrcan-rncan.gc.ca/geod/tools-outils/nad83-docs.php

As a datapoint, MassGIS releases imagery which is natively in
"NAD83(2011) epoch 2010.0" EPGS:6319 (in UTM, but that is not a datum
change).  They document checking alignment with ground control and an
RMS error in E and N of about 0.1m (with 0.15m pixel size).

They also release this imagery as TMS, and hence in 3857.  In the proj
world, there is no transform from 6319 to the WGS84 ensemble.  But,
transforming 6319 to ITRF2014 (as a proxy for WGS84(G2159) has a shift
of about a meter.  The TMS data has in fact been transformed.  I asked
about this (did you transform, and to what) and have not gotten an
answer.

For me, this is not academic: I have actually compared ground positions
from RTK (in EPSG:6319) and imagery, and treating 6319 and ITRF2014
(7912) as distinct is necessary to get them to line up in both an NAD83
GIS view and in a "web mercator" openstreetmap tools view.

I deal with this by converting to ITRF2014 as an intermediate step, to
avoid the null transform from 6319.

There's a further wrinkle that WGS84 is a dynamic datum and crust-fixed
points have velocities in the 2-3 cm/y range.  That's a whole 15 cm
pixel in 5 years.  So one could ask about the epoch in TMS, and I
haven't found that documented either.

This is why I am asking you about the actual datum of your data that you
say is in 3857.   I would expect that's a far larger error source, as
data in 3857 now I would argue should be in the most recent realization
of WGS84, and using either epoch of data, 2022.5 (to align with GPS
control segment practice), or the reference epoch of ITRF2014, which is
as far as I can tell the matching ITRF from G2139.

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 194 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/mapproxy/attachments/20220618/a8cfa8c4/attachment.sig>