Another WCS question

[maelle]

@eblondel in the reprex below

• I download the whole coverage from a service. This works.
• I try to subset it using a bounding box that intersects with both the metadata and data bounding boxes. It does not work.
• I try to subset it using a bounding box, specifying the CRS. It does not work either.

I wonder what I am doing wrong. Looking at the GML enveloppe messages I wonder whether there's an issue with how I'm passing the bounding box / CRS?

library("ows4R")
#> Loading ISO 19139 XML schemas...
#> Loading ISO 19115-3 XML schemas...
#> Loading ISO 19139 codelists...

# Connect to service
client <- ows4R::WCSClient$new(
  "https://ows.emodnet-humanactivities.eu/wcs",
  serviceVersion = "2.0.1",
  headers = c(
    "User-Agent" = "emodnet.wcs R package https://github.com/EMODnet/emodnet.wcs"
  )
)

# Download coverage
summary <- client$getCapabilities()$findCoverageSummaryById(
  "emodnet__vesseldensity_all"
)
cov <- summary$getCoverage(
  bbox = NULL,
  crs = NULL,
  time = NULL,
  elevation = NULL,
  format = NULL,
  rangesubset = NULL,
  filename = "reprex.tif"
)
#> Loading required package: sf
#> Linking to GEOS 3.12.2, GDAL 3.10.3, PROJ 9.4.1; sf_use_s2() is FALSE
#> No encoding supplied: defaulting to UTF-8.
#> <GMLEnvelope>
#> ....|-- lowerCorner: -9763441.13028937 1684718.73495521 "2017-01-01T00:00:00"
#> ....|-- upperCorner: 10868458.8067861 22822359.9215428 "2023-12-01T00:00:00"

# Metadata bbox
summary$getBoundingBox()$BoundingBox$getBBOX()
#>     xmin     ymin     xmax     ymax 
#> -9763441  1684719 10868459 22822360
# Actual bbox
sf::st_bbox(cov)
#>      xmin      ymin      xmax      ymax 
#> -13053174  -1887757  13640295  27500483

# Metadata bbox (WFS84)
summary$getWGS84BoundingBox()$WGS84BoundingBox$getBBOX()
#>      xmin      ymin      xmax      ymax 
#> -87.70648  14.96111  97.63303  86.80081
# Actual bbox (WFS84)
sf::st_bbox(cov) |>
  sf::st_transform(crs = "EPSG:4326")
#>       xmin       ymin       xmax       ymax 
#> -117.25866  -16.71571  122.53286   88.46330


# Does it intersect user's bbox
user_bbox <- c(xmin = 4.2, ymin = 53, xmax = 8.8, ymax = 54) |>
  sf::st_bbox() |>
  sf::st_as_sfc()

sf::st_crs(user_bbox) <- "EPSG:4326"

cov_bbox <- sf::st_bbox(cov) |>
  sf::st_transform(crs = "EPSG:4326") |>
  sf::st_as_sfc()

metadata_bbox <- summary$getWGS84BoundingBox()$WGS84BoundingBox$getBBOX() |>
  sf::st_as_sfc()

sf::sf_use_s2(TRUE)
#> Spherical geometry (s2) switched on
sf::st_intersects(user_bbox, metadata_bbox)
#> Sparse geometry binary predicate list of length 1, where the predicate
#> was `intersects'
#>  1: 1
sf::st_intersects(user_bbox, cov_bbox)
#> Sparse geometry binary predicate list of length 1, where the predicate
#> was `intersects'
#>  1: 1

# Two different bounding boxes
# metadata_bbox: in the metadata
# cov_bbox: in the data
# both intersect the user's bbox

# However no data if we use the user's bbox
# Why?!

bbox <- c(xmin = 4.2, ymin = 53, xmax = 8.8, ymax = 54)
ows_bbox <- ows4R::OWSUtils$toBBOX(
  xmin = bbox["xmin"],
  xmax = bbox["xmax"],
  ymin = bbox["ymin"],
  ymax = bbox["ymax"]
)
ows_bbox
#>    min  max
#> x  4.2  8.8
#> y 53.0 54.0
summary$getCoverage(
  bbox = ows_bbox,
  crs = NULL,
  time = NULL,
  elevation = NULL,
  format = NULL,
  rangesubset = NULL,
  filename = "reprex.tif"
)
#> <GMLEnvelope>
#> ....|-- lowerCorner: 4.2 53 "2017-01-01T00:00:00"
#> ....|-- upperCorner: 8.8 54 "2023-12-01T00:00:00"
#> Warning in new_CppObject_xp(fields$.module, fields$.pointer, ...): GDAL Error
#> 4: `/tmp/RtmpkT1950/reprex-66263506e32c-wax-mouse/reprex.tif' not recognized as
#> being in a supported file format.
#> Error: [rast] cannot open this file as a SpatRaster: /tmp/RtmpkT1950/reprex-66263506e32c-wax-mouse/reprex.tif
brio::read_lines("reprex.tif")
#> [1] "<?xml version=\"1.0\" encoding=\"UTF-8\"?><ows:ExceptionReport xmlns:xs=\"http://www.w3.org/2001/XMLSchema\" xmlns:ows=\"http://www.opengis.net/ows/2.0\" xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\" version=\"2.0.0\" xsi:schemaLocation=\"http://www.opengis.net/ows/2.0 http://schemas.opengis.net/ows/2.0/owsExceptionReport.xsd\">"
#> [2] "  <ows:Exception exceptionCode=\"InvalidSubsetting\" locator=\"\">"                                                                                                                                                                                                                                                                                  
#> [3] "    <ows:ExceptionText>Empty intersection after subsetting</ows:ExceptionText>"                                                                                                                                                                                                                                                                      
#> [4] "  </ows:Exception>"                                                                                                                                                                                                                                                                                                                                  
#> [5] "</ows:ExceptionReport>"

summary$getCoverage(
  bbox = ows_bbox,
  crs = "EPSG:4326",
  time = NULL,
  elevation = NULL,
  format = NULL,
  rangesubset = NULL,
  filename = "reprex.tif"
)
#> <GMLEnvelope>
#> ....|-- lowerCorner: 4.2 53 "2017-01-01T00:00:00"
#> ....|-- upperCorner: 8.8 54 "2023-12-01T00:00:00"
#> Warning in new_CppObject_xp(fields$.module, fields$.pointer, ...): GDAL Error
#> 4: `/tmp/RtmpkT1950/reprex-66263506e32c-wax-mouse/reprex.tif' not recognized as
#> being in a supported file format.
#> Error: [rast] cannot open this file as a SpatRaster: /tmp/RtmpkT1950/reprex-66263506e32c-wax-mouse/reprex.tif
brio::read_lines("reprex.tif")
#> [1] "<?xml version=\"1.0\" encoding=\"UTF-8\"?><ows:ExceptionReport xmlns:xs=\"http://www.w3.org/2001/XMLSchema\" xmlns:ows=\"http://www.opengis.net/ows/2.0\" xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\" version=\"2.0.0\" xsi:schemaLocation=\"http://www.opengis.net/ows/2.0 http://schemas.opengis.net/ows/2.0/owsExceptionReport.xsd\">"
#> [2] "  <ows:Exception exceptionCode=\"InvalidSubsetting\" locator=\"\">"                                                                                                                                                                                                                                                                                  
#> [3] "    <ows:ExceptionText>Empty intersection after subsetting</ows:ExceptionText>"                                                                                                                                                                                                                                                                      
#> [4] "  </ows:Exception>"                                                                                                                                                                                                                                                                                                                                  
#> [5] "</ows:ExceptionReport>"

^{Created on 2025-09-19 with reprex v2.1.1}