PROPOSAL: proj_transform_bounds

[snowman2]

Related:
pyproj4/pyproj#809
https://lists.osgeo.org/pipermail/gdal-dev/2021-July/054450.html

Any interest in porting the logic from pyproj into PROJ?

[rouault]

I guess there is an interest for that, but there's a number of subtelties to consider. Would we want to return just a rectangle, or a more complex polygon ? What about of a rectangle in a CRS that crosses the antimeridian (let's say in EPSG:32601 or 32660) and transformed to EPSG:4326 ? Or a rectangle in a polar stereographic projection that includes the pole and is transformed to EPSG:4326 etc.

[snowman2]

This is how pyproj handles it: https://github.com/pyproj4/pyproj/blob/dd84df30c0ccee371334a8c89f7c14f86b98c49c/test/test_transformer.py#L1279

The box coordinates are always returned, but for the antimeridian, minx>maxx.

[snowman2]

Then users can create two polygons if they desire: rasterio/rasterio#2222 (comment)

[snowman2]

A more complex polygon returned would likely be more user friendly.

[snowman2]

Or a rectangle in a polar stereographic projection that includes the pole and is transformed to EPSG:4326

How can you know if a rectangle in the polar stereographic projection includes the pole?

[rouault]

How can you know if a rectangle in the polar stereographic projection includes the pole?

Completely opaque solution: project the pole (0,90) and see if the projected coordinate is in the projected rectangle.
If you know it is polar stereographic projection, I guess checking if (false_easting, false_northing) is in the projected rectangle should be OK since the pole is the projection center

[snowman2]

Thanks for your response @rouault 👍

project the pole (0,90) and see if the projected coordinate is in the projected rectangle

This was the solution I tried. Saw strange behavior, so I am continuing the investigation for this approach.

[snowman2]

Here is the part that I am uncertain about. I am hoping you could help me understand this better.

If I create an area in the range 80-85 latitude and generate points in WGS 84 / UPS North (N,E):

import xarray
import numpy
from pyproj import Transformer

transformer = Transformer.from_crs("EPSG:4326", "EPSG:32661", always_xy=True)

xpole, ypole = transformer.transform(0, 90)  # (2000000.0, 2000000.0)
lon = numpy.arange(-180, 180, 1)
lat = numpy.arange(80, 85, 1)
lon2d, lat2d = numpy.meshgrid(lon, lat)
xx, yy = transformer.transform(lon2d, lat2d)
xds = xarray.Dataset(
    coords={"x": lon, "y": lat},
    data_vars={
        "xx": (('y', 'x'), xx),
        "yy": (('y', 'x'), yy)
    }
)
xds.xx.min().item(), xds.xx.max().item()  # (887048.8630450945, 3112951.1369549055)
xds.yy.min().item(), xds.yy.max().item()  # (887048.8630450945, 3112951.1369549055)

The pole appears to be within the min/max values for the projected X and Y coordinates.

(xds.xx.min() < xpole < xds.xx.max()).item()   # True
(xds.yy.min() < ypole < xds.yy.max()).item()   # True

[rouault]

The resulting shape must be a circular "polygon" whose external edge is latitude=80 with a hole starting at latitude=85. So if you only consider the bounding box of that circle, the pole fits into it, but not in the polygon since it is in its hole.

[snowman2]

@rouault you are indeed correct. I was thinking about it in 3D lat/lon space and it didn't make sense. But, if you flatten it out in 2D space and change the perspective, it definitely makes sense.

import numpy
import geopandas
from shapely.geometry import Point, Polygon

min_lon = -180
max_lon = 180
min_lat = 80
max_lat = 85
lons = numpy.arange(min_lon, max_lon, 1)
lats =  numpy.arange(min_lat, max_lat, 1)

lonbox = numpy.concatenate((
    numpy.ones(len(lats)) * min_lon,
    lons,
    numpy.ones(len(lats)) * max_lon,
    lons[::-1],
))
    
    
latbox = numpy.concatenate((
    lats,
    numpy.ones(len(lons)) * max_lat,
    lats[::-1],
    numpy.ones(len(lons)) * min_lat,    
    
))
xbox, ybox = transformer.transform(lonbox, latbox)
xpole, ypole = transformer.transform(0, 90)
gs = geopandas.GeoSeries([Point(xpole, ypole), Polygon(zip(xbox, ybox))])
gs.plot()

[snowman2]

Thanks again, @rouault 👍 Based on that input, I have a PR with updates for polar support in transform_bounds. Any suggestions welcome.

Would the implementation in pyproj for transform_bounds be of interest for PROJ or do you have an alternate solution you would prefer?

[rouault]

I guess the approach of pyproj would be a reasonable start. I would probably drop the always_xy / is_radians logic I can see from the function prototype and would expect input and output to be with the axis order and units of the input & target CRS. The user can "normalize" its CRS if wanted outside of this function.