Reversing bounding box projection results in different extents

[natemiller]

I'm trying to understand why if I project a bounding box into a new CRS and then use the resulting extents (xmin, ymin, xmax, ymax) to regenerate a bounding box and transform back to the original CRS the resulting extents (x-coordinates) are different than what I started with?

# transform bounding box defined in ESPG:4326 to Equal Earth

bb_ee = sf::st_bbox(c(xmin = -20, xmax = 30, 
                      ymin = 30, ymax = 65),  
                 crs = 4326) %>%
  sf::st_as_sfc() %>%
  sf::st_transform(crs = "+proj=eqearth +datum=WGS84 +wktext") %>%
  sf::st_bbox()

# Using extent of Equal Earth bounding box, transform back to EPSG:4326

bb_4326 = sf::st_bbox(c(xmin = bb_ee[['xmin']], xmax = bb_ee[['xmax']], 
                        ymin = bb_ee[['ymin']], ymax = bb_ee[['ymax']]),  
                 crs = "+proj=eqearth +datum=WGS84 +wktext") %>%
  sf::st_as_sfc() %>%
  sf::st_transform(crs = 4326) %>%
  sf::st_bbox()

# input is
#   xmin      ymin      xmax      ymax 
#   -20       30        30        65 

bb_4326

# result is  
#   xmin         ymin           xmax          ymax 
# -26.12724  30.00000  39.19085  65.00000

[edzer]

No idea, but with

> library(sf)
Linking to GEOS 3.8.1, GDAL 3.1.3, PROJ 7.1.1

I see for the first transformation

Warning messages:
1: In CPL_crs_parameters(x) :
  GDAL Error 1: PROJ: proj_as_wkt: Unsupported conversion method: Equal Earth
2: In CPL_crs_parameters(x) :
  GDAL Error 1: PROJ: proj_as_wkt: Unsupported conversion method: Equal Earth

Didn't you see those?

[natemiller]

No I have not seen those errors.

I'm clearly lagging behind here though, mainly avoiding the time I expect to lose updating GDAL.

> library(sf)
Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0

[edzer]

Please try a similar experiment with GDAL standalone tools, or PROJ standalone tools, in order to identify whether your issue is with sf or with the underlying software (GDAL and/or PROJ).

[rsbivand]

Interesting. The bounding box of the projected object - which is trapezoidal:

ends up being inverse trapezoidal when inverse projected:

The projection seems to be EPSG:8857, and I could run something like the example by coercing to sp and using rgdal with current PROJ/GDAL, yielding the figures. So, the change in the bounding box is expected. I didn't explore why sf is showing the observed warnings.

[natemiller]

Thanks for exploring that. Seems bizarre that a rectangular input object returns to the same projection as a trapezoid. I didn't realize that it wasn't possible to reverse projections in this manner and it throws a wrench into an approach I was planning.

[rsbivand]

The GEOGCRS geometry itself will round-trip project/inverse-project, but the inverse projected bbox/envelope of the projected geometry will in general not equal the bbox/envelope of the input GEOGCRS. There may be some projections that do this, but in general something gets stretched.

[natemiller]

Interesting... so if I have a map (in this case I have one in produced in ggplot) with coordinates defined in one projection and the associated extents (xmin,ymin, xmax, ymax) , I cannot determine what the extents of that map might be in a different projection? That's frustrating.

[edzer]

st_graticule segmentizes the polygon formed by the bounding rectangle, and transforms that to get the coordinate range in the transformed space. That shape, projected, is typically no longer a box. And even that approach is far from fool proof... (just think about covering a pole or antimeridian).

[natemiller]

I certainly understand that the shape in the projected space is not a box. That makes sense (the trapezoid in equal earth makes sense though if we had more vertices, using something like smoothr::densify(n =100L), the lines would be curved). But it seems like the extents of that shape should still be the same, just in a new coordinate system. I'd then expect that given that the new shape should have the same extents, albeit in a different set of coordinates, if I project those extents back to the original coordinate system I should end up where I started.

Clearly this doesn't happen, but it's hard for me to understand how, if the "box" in the new coordinate space represents the same box in the original space, it could have different extents... and interestingly in my example only having shifted extents in the x-axis.

We don't have to continue this here, but I don't feel any of the explanations so far have explained why forward and backward projecting of extents doesn't work (unless the maximum and minimum coordinates in the new projection are not equivalent to the maximum and minimum coordinates in the old projection, but if that's the case then it seems the new projection is inaccurate. Because one of the fundamental features should be, it seems, that a projected shape accurately represents the fundamental features of original shape in its own projection, though the appearance, such as curved lines, may be different).

I likely just need to read more about how these processes work, because some portion of my understanding seems flawed.

[edzer]

Your confusion is caused by recomputing of the bounding box:

library(sf)
# Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
sf::st_bbox(c(xmin = -20, xmax = 30, 
                      ymin = 30, ymax = 65),  
                 crs = 4326) %>%
  sf::st_as_sfc() %>%
  sf::st_transform(crs = "+proj=eqearth +datum=WGS84 +wktext") %>%
  sf::st_transform(crs = 4326)
# Geometry set for 1 feature 
# geometry type:  POLYGON
# dimension:      XY
# bbox:           xmin: -20 ymin: 30 xmax: 30 ymax: 65
# geographic CRS: WGS 84
# POLYGON ((-20 30, 30 30, 30 65, -20 65, -20 30))

shows that the projections is invertible; if you recompute the bbox in the projected space, the min/max values are there recomputed over all points, not just taken from the lower-left and upper right. Roger's post supports that the corner point's roles switch, for the x axis.