[MapProxy] shift when reprojecting WMS

[Travis Kirstine]

Greg,

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
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.

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.





On Thu, 16 Jun 2022 at 11:54, Greg Troxel <gdt at lexort.com> wrote:

>
> Travis Kirstine <traviskirstine at gmail.com> writes:
>
> > We have a fully seeded cache in a Google Maps / EPSG:3857 SRS and offer
> the
> > data up through a WMS with various SRS supported,  MapProxy is handling
> the
> > reprojection.  We've noticed that with our latest upgrade to mapproxy
> v14.0
> > with pyproj 3.0.1 / proj 7.2.1 there is an error with some SRS being
> > shifted incorrectly, specifically EPSG:2958 / UTM Zone17 CSRS.  How does
> > mapproxy handle the image reprojection - is it using pyproj or some other
> > library?
>
> What do you think is happening, and why?
> What do you think shoudl be happening?
>
> CSRS is presumably NAD83(CSRSv?).
>
> 3857 refers to the WGS84 ensemble, which includes WGS84(TRANSIT) which
> is approximately equal to NAD83(1986) and presumably then to NAD83(CSRS)
> (the unsuffixed version).
>
> So I'd expect no transform.  Arguably, there should be one, as the
> imagery probably ought to be in spherical mercator from WGS84(G2139)
> which is distinct from NAD83(CSRSvRECENT).
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/mapproxy/attachments/20220617/747eb412/attachment.htm>