Suspicious bounds in

[schwehr]

Describe the bug

-280,-90,180,90 is suspicious with long going from -280 to 180

sf/R/transform.R

Line 92 in a5f9f0f

#' try(st_transform(sfc, 3857, aoi = c(-280,-90,180,90)))

has

#' st_transform(st_sf(a=2:1, geom=sfc), "+init=epsg:3857")
#' try(st_transform(sfc, 3857, aoi = c(-280,-90,180,90)))
#' if (sf_extSoftVersion()["GDAL"] >= "3.0.0") {
#'   st_transform(sfc, pipeline =
#' 	  "+proj=pipeline +step +proj=axisswap +order=2,1") # reverse axes
#'   st_transform(sfc, pipeline =
#' 	  "+proj=pipeline +step +proj=axisswap +order=2,1", reverse = TRUE) # also reverse axes
#' }

Got this: GDAL Error 1: Invalid dfWestLongitudeDeg

  1. st_transform.Rd
    Output prior to failure:
      > ### Name: st_transform
      > ### Title: Transform or convert coordinates of simple feature
      > ### Aliases: st_transform st_transform.sfc st_transform.sf st_transform.sfg
      > ###   st_wrap_dateline st_wrap_dateline.sfc st_wrap_dateline.sf
      > ###   st_wrap_dateline.sfg sf_proj_info
      > 
      > ### ** Examples
      > 
      > p1 = st_point(c(7,52))
      > p2 = st_point(c(-30,20))
      > sfc = st_sfc(p1, p2, crs = 4326)
      > sfc
      Geometry set for 2 features 
      Geometry type: POINT
      Dimension:     XY
      Bounding box:  xmin: -30 ymin: 20 xmax: 7 ymax: 52
      Geodetic CRS:  WGS 84
      POINT (7 52)
      POINT (-30 20)
      > st_transform(sfc, 3857)
      Geometry set for 2 features 
      Geometry type: POINT
      Dimension:     XY
      Bounding box:  xmin: -3339585 ymin: 2273031 xmax: 779236.4 ymax: 6800125
      Projected CRS: WGS 84 / Pseudo-Mercator
      POINT (779236.4 6800125)
      POINT (-3339585 2273031)
      > st_transform(st_sf(a=2:1, geom=sfc), "+init=epsg:3857")
      Warning in CPL_crs_from_input(x) :
        GDAL Message 1: +init=epsg:XXXX syntax is deprecated. It might return a CRS with a non-EPSG compliant axis order.
      Simple feature collection with 2 features and 1 field
      Geometry type: POINT
      Dimension:     XY
      Bounding box:  xmin: -3339585 ymin: 2273031 xmax: 779236.4 ymax: 6800125
      Projected CRS: WGS 84 / Pseudo-Mercator
        a                     geom
      1 2 POINT (779236.4 6800125)
      2 1 POINT (-3339585 2273031)
      > try(st_transform(sfc, 3857, aoi = c(-280,-90,180,90)))
      Warning in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy,  :
        GDAL Error 1: Invalid dfWestLongitudeDeg
      Error : values for area of interest not accepted
      > if (sf_extSoftVersion()["GDAL"] >= "3.0.0") {
      +   st_transform(sfc, pipeline =
      + 	  "+proj=pipeline +step +proj=axisswap +order=2,1") # reverse axes
      +   st_transform(sfc, pipeline =
      + 	  "+proj=pipeline +step +proj=axisswap +order=2,1", reverse = TRUE) # also reverse axes
      + }

To Reproduce

Ran the tests using PROJ ~7.0 and GDAL ~3.4

• GDAL at d2f9067ffb15e593e9b826ca939dbd183636c780
• PROJ at cbf2eebf1bcb96408ebb20b8772a079d977d9fa2

I'm upgrading PROJ from 0ea51a6a2bef5c0215bc5f475473fcb5044e01f4 (6.3.2)

If reporting a change from previous versions

Please read https://cran.r-project.org/web/packages/sf/news/news.html first.

Additional context

I'm working in a custom bazel environment, so it's hard to provide a reproducer info.

The test was with sf 1.0-9, but the -280 is in head.

[schwehr]

After converting -280 → -180, the full output is:

  1. st_transform.Rd
    Output prior to failure:
      > ### Name: st_transform
      > ### Title: Transform or convert coordinates of simple feature
      > ### Aliases: st_transform st_transform.sfc st_transform.sf st_transform.sfg
      > ###   st_wrap_dateline st_wrap_dateline.sfc st_wrap_dateline.sf
      > ###   st_wrap_dateline.sfg sf_proj_info
      > 
      > ### ** Examples
      > 
      > p1 = st_point(c(7,52))
      > p2 = st_point(c(-30,20))
      > sfc = st_sfc(p1, p2, crs = 4326)
      > sfc
      Geometry set for 2 features 
      Geometry type: POINT
      Dimension:     XY
      Bounding box:  xmin: -30 ymin: 20 xmax: 7 ymax: 52
      Geodetic CRS:  WGS 84
      POINT (7 52)
      POINT (-30 20)
      > st_transform(sfc, 3857)
      Geometry set for 2 features 
      Geometry type: POINT
      Dimension:     XY
      Bounding box:  xmin: -3339585 ymin: 2273031 xmax: 779236.4 ymax: 6800125
      Projected CRS: WGS 84 / Pseudo-Mercator
      POINT (779236.4 6800125)
      POINT (-3339585 2273031)
      > st_transform(st_sf(a=2:1, geom=sfc), "+init=epsg:3857")
      Warning in CPL_crs_from_input(x) :
        GDAL Message 1: +init=epsg:XXXX syntax is deprecated. It might return a CRS with a non-EPSG compliant axis order.
      Simple feature collection with 2 features and 1 field
      Geometry type: POINT
      Dimension:     XY
      Bounding box:  xmin: -3339585 ymin: 2273031 xmax: 779236.4 ymax: 6800125
      Projected CRS: WGS 84 / Pseudo-Mercator
        a                     geom
      1 2 POINT (779236.4 6800125)
      2 1 POINT (-3339585 2273031)
      > try(st_transform(sfc, 3857, aoi = c(-180,-90,180,90)))
      Geometry set for 2 features 
      Geometry type: POINT
      Dimension:     XY
      Bounding box:  xmin: -3339585 ymin: 2273031 xmax: 779236.4 ymax: 6800125
      Projected CRS: WGS 84 / Pseudo-Mercator
      POINT (779236.4 6800125)
      POINT (-3339585 2273031)
      > if (sf_extSoftVersion()["GDAL"] >= "3.0.0") {
      +   st_transform(sfc, pipeline =
      + 	  "+proj=pipeline +step +proj=axisswap +order=2,1") # reverse axes
      +   st_transform(sfc, pipeline =
      + 	  "+proj=pipeline +step +proj=axisswap +order=2,1", reverse = TRUE) # also reverse axes
      + }
      Geometry set for 2 features 
      Geometry type: POINT
      Dimension:     XY
      Bounding box:  xmin: -30 ymin: 20 xmax: 7 ymax: 52
      CRS:           NA
      POINT (7 52)
      POINT (-30 20)
      > # local modification(chiricom): Requires 'lwgeom'
      > ## Not run: 
      > ##D nc = st_read(system.file("shape/nc.shp", package="sf"))
      > ##D st_area(nc[1,]) # area from long/lat
      > ##D st_area(st_transform(nc[1,], 32119)) # NC state plane, m
      > ##D st_area(st_transform(nc[1,], 2264)) # NC state plane, US foot
      > ##D library(units)
      > ##D set_units(st_area(st_transform(nc[1,], 2264)), m^2)
      > ## End(Not run)
      > st_transform(structure(p1, proj4string = "+init=epsg:4326"), "+init=epsg:3857")
      Warning in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy,  :
        GDAL Error 1: PROJ: proj_create_operations: At least one of the operation lacks a source and/or target CRS
      Warning in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy,  :
        GDAL Error 6: Cannot find coordinate operations from `GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]],CS[ellipsoidal,2],AXIS["longitude",east,ORDER[1],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],AXIS["latitude",north,ORDER[2],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],USAGE[SCOPE["unknown"],AREA["World"],BBOX[-90,-180,90,180]]]' to `+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crs'
      Error: OGRCreateCoordinateTransformation(): transformation not available

[edzer]

Thanks - I think I wanted to trigger that error (hence try()), but it doesn't make much sense.