Understanding the difference between transform() and transform_bounds()

[tlarcher]

Hello community !

Problem description
Despite the documentation and what I could read online I'm having a hard time understanding how to interpret the difference between transform() and transform_bounds() methods of a pyproj.Transformer object.

I have 2 pyproj.CRS objects :

• "WGS 84 / UTM zone 31N" : crs32361 = pyproj.CRS.from_epsg(32361)
• "WGS 84" : crs4326 = crs32631.geodetic_crs
crs32361 has a bounds attribute which gives the WGS84 bounds of the CRS, which is :
• east: 6.0
• north: 84.0
• south: 0.0
• west: 0.0
  I need to project (east, north) and (west, south) in the EPSG:32631 for various computations. To do so I could:

transformer = pyproj.Transformer.from_crs(crs32631.geodetic_crs, crs32631)

1. Use transform() on each point:

bounds_32631 = transformer.transform([0, 6], [0, 84])  #([array_of_x], [array_of_y])
print(bounds_32631)
>>> ([166021.44308054057, 534994.6550611363], [0.0, 9329005.182447437])  # ([xmin, xmax], [ymin, ymax])

2. Use transform_bounds() which differs by "densifying the edges to account for nonlinear transformations along these edges and extracting the outermost bounds":

bounds_32631_tb = transformer.transform_bounds(0, 0, 6, 84)  #(left, bottom, right, top)
print(bounds_32631_tb)
>>> (166021.44308054057, 0.0, 833978.5569194613, 9329005.182447437)  #(xmin, ymin, xmax, ymax)

$\rightarrow$ matches with epsg.io and spatial reference.

---

Furthermore, when I compute the inverse transforms, more inconsistencies appear:

1. Using inverse transform() on points projected via transform():

bounds_32631_inv = transformer.transform(*bounds_32631, direction="INVERSE", errcheck=True)  #([array_of_x], [array_of_y])
print(bounds_32631_inv)
>>> ([3.975693351829396e-16, 6.0], [0.0, 84.0])  # ([xmin, xmax], [ymin, ymax])

OK (up to numerical approximations)

2. Using inverse transform_bounds() on points projected via transform_bounds():

bounds_32631_tb_inv = transformer.transform_bounds(*bounds_32631_tb, direction="INVERSE", errcheck=True)  #(left, bottom, right, top)
print(bounds_32631_tb_inv)
>>> (-23.58287765881726, 0.0, 29.58287765881738, 84.00816352991251)  #(xmin, ymin, xmax, ymax)

KO ($x$ values are wrong)

3. Using inverse transform_bounds() on points projected via transform():

bounds_32631_inv_tb = transformer.transform_bounds(bounds_32631[0][0], bounds_32631[1][0], bounds_32631[0][1], bounds_32631[1][1], direction="INVERSE", errcheck=True)  #(left, bottom, right, top)
print(bounds_32631_inv_tb)
>>> (-23.58287765881726, 0.0, 6.0, 84.00814947301849)  #(xmin, ymin, xmax, ymax)

KO ($x_{min}$ value is wrong)

4. Using inverse transform() on points projected via transform_bounds():

bounds_32631_tb_inv = transformer.transform([bounds_32631_tb[0], bounds_32631_tb[2]], [bounds_32631_tb[1], bounds_32631_tb[3]], direction="INVERSE", errcheck=True)  #([array_of_x], [array_of_y])
print(bounds_32631_tb_inv)
>>> ([3.975693351829396e-16, 29.58287765881738], [0.0, 83.30591198713375])  # ([xmin, xmax], [ymin, ymax])

KO ($x_{max}$ value is wrong, others are OK up to numerical approximations)

Questions
Why does the $x$ component of the (east, north) point differ in the FORWARD transformation ?
Does that mean that the bounding box in the EPSG:32361 is actually a trapeze ?
Will it be impossible to find points within the crs32361.geodetic_crs bounds falling between $534994.6550611363$ and $833978.556919461$ on the projected $x$ axis ?
Considering what happens in INVERSE transforms, how will using both transform for points and transform_bounds for bbox NOT cause further problems ?

Installation
Python 3.10.6
pyproj==3.4.0
odc-geo==0.3.2

[snowman2]

Related: #809

This is likely a helpful reference to visually demonstrate what is going on: ref

[tlarcher]

Thank you for the provided resources. I understand well that non-linear transformations are implied in some CRS projections.

However I'm also interested in the underlying reasons of these transformations.
Could one imagine that a CRS would not be linearly transformable to include land (e.g.: islands) outside of the "main" area ? Or for the CRS boundaries to follow natural boundaries (e.g.: mountains) ?

[tlarcher]

Also, since the correct projected bounds obtained using transform_bounds() give a larger or equal area than if we were to project each bound, this means that in order to project to this "extra space" from the CRS' geodetic CRS, one would need to use a point outside of the CRS' area_of_use geodetic bounds.
Isn't that bothering when working with a set of projections like, for instance, EPSG:32601 to 32661 (northern hemisphere split by 6° UTM zones) ?

[snowman2]

Could one imagine that a CRS would not be linearly transformable to include land (e.g.: islands) outside of the "main" area ? Or for the CRS boundaries to follow natural boundaries (e.g.: mountains) ?

That is not necessarily what is happening. The shape of the boundary in one projection can be quite different in another projection while covering the same locations. The main difference is that the boundary of the shape cannot always be estimated using the corners. All transform_bounds is doing is taking the boundary, adding additional points along the edges of the boundary, re-projecting all of the points and taking the envelope of the new shape.

Here is a function to add points between the corners of a bounding box:

import numpy

def densify_box(minx, miny, maxx, maxy):
    """
    Add intermediate points along the edges of the box.
    """
    xx = numpy.arange(minx, maxx, 1)
    yy =  numpy.arange(miny, maxy, 1)

    xbox = numpy.concatenate((
        numpy.ones(len(yy)) * minx,
        xx,
        numpy.ones(len(yy)) * maxx,
        xx[::-1],
    ))

    ybox = numpy.concatenate((
        yy,
        numpy.ones(len(xx)) * maxy,
        yy[::-1],
        numpy.ones(len(xx)) * miny,    

    ))
    return xbox, ybox

Using this function, you can generate a polygon of the bounding box after it has transformed to different projections:

import geopandas
from pyproj import CRS, Transformer
from shapely.geometry import Polygon, box

proj_crs = CRS("EPSG:3035")
min_lon, min_lat, max_lon, max_lat = proj_crs.area_of_use.bounds
transformer = Transformer.from_crs("EPSG:4326", proj_crs, always_xy=True)
xbox, ybox = transformer.transform(*densify_box(min_lon, min_lat, max_lon, max_lat))
projected_bounding_polygon = Polygon(zip(xbox, ybox))

Then, you can compare calculating the bounds using the original corners:

xx, yy = transformer.transform(
    [min_lon, max_lon, min_lon, max_lon],
    [min_lat, min_lat, max_lat, max_lat]
)
simple_bounding_polygon = box(min(xx), min(yy), max(xx), max(yy))

Comparing against the projected bounding polygon, you can see that the simple bounds do not contain the projected polygon:

gs = geopandas.GeoSeries([simple_bounding_polygon, projected_bounding_polygon])
gs.plot(cmap="viridis")

That is where transform_bounds comes in:

real_bounds = box(*transformer.transform_bounds(min_lon, min_lat, max_lon, max_lat))

You can see that it contains the projected bounding shape:

gs = geopandas.GeoSeries([real_bounds, projected_bounds])
gs.plot(cmap="viridis")

[tlarcher]

Ok, I understand much better now with this example. I had mixed up the purposes of transform and transform_bounds as the former projects points in the destination CRS space, whereas the later simply reconstructs the correct projected boundaries of the projected space, taking deformations (non-linear transformations) into account (as indicated in its name, obviously enough...).

I also realize now that in order to test whether a point or an area in WGS84 belongs a given projected space, it doesn't make much sense to use the projected boundaries; rather for this purpose proj_crs.area_of_use.bounds is better suited (provided that WG84 is the geodetic crs of the projected space in question).

Finally, to assess membership of a point in an unknown CRS to a set of unknown 2D-projected CRS like EPSG:3035, one would have to first project the source point to each 2D-projected CRS' geodetic CRS and compare the values of the projected point with each area_of_use.bounds. Correct ?

[snowman2]

Finally, to assess membership of a point in an unknown CRS to a set of unknown 2D-projected CRS like EPSG:3035, one would have to first project the source point to each 2D-projected CRS' geodetic CRS and compare the values of the projected point with each area_of_use.bounds. Correct ?

That sounds correct to me.