cov$getCoverageStack() returns file with "The coordinate reference system must be the same for all objects." exception

[tphilippi]

I am not clear if this is my ignorance, a problem with this specific OGC data service, or something about getting a set of temporal coverages via getCoverageStack(), but I'm still floundering trying to pull a stack of annual rasters. I have read all of the discussions tagged as OGC WCS, which were helpful but did not solve my problem.

This data service provides rasters of land cover classification for each year from 1985-2023. I'm guessing at the serviceVersion from the xml metadata: xmlns:ows="http://www.opengis.net/ows/2.0" but have tried not just 2.0.0, but 2.0.1, and 1.0.0.
I can successfully use getCoverage() to pull single rasters from their other data services that vend single rasters, which also use epsg:5070, so that crs alone is not the problem.

Reproducible example:

library(terra)
library(ows4R)
# US National Land Cover Data Services.  See "Annual NLCD" in:
#  https://www.mrlc.gov/data-services-page

bbox5070 <- array(c(107392, 779125, 191640, 862413), 
                  dim = c(2,2), 
		  dimnames = list(c("x", "y"), c("min", "max")))
base <- "https://dmsdata.cr.usgs.gov/geoserver/mrlc_Land-Cover-Native_conus_year_data/wcs?"				  
wcs <- WCSClient$new(base,
                     serviceVersion = "2.0.0", 
                     logger = "INFO")

caps <- wcs$getCapabilities()
cov  <-caps$getCoverageSummaries()[[1]]
cov$getDimensions()
wcs$describeCoverage(cov$CoverageId)

tmp <- cov$getCoverageStack(bbox = bbox5070,
                            crs = "EPSG:5070")
sessionInfo()

The above code fails with an error that it cannot open the tiff file.
R error message:

[ows4R][INFO] WCSCoverageSummary - Fetching coverage stack with 'temporal' dimension
[ows4R][INFO] WCSCoverageSummary - 1985-01-01T00:00:00
[ows4R][WARN] WCSCoverageSummary - Coverage without vertical dimension: 'elevation' argument is ignored

....|-- lowerCorner: 107392 779125 "1985-01-01T00:00:00"
....|-- upperCorner: 191640 862413 "2023-01-01T00:00:00"[ows4R][INFO] WCSGetCoverage - Fetching
https://dmsdata.cr.usgs.gov/geoserver/mrlc_Land-Cover-Native_conus_year_data/wcs?service=WCS&version=2.0.0&coverageId=mrlc_Land-Cover-Native_conus_year_data__Land-Cover-Native_conus_year_data&subset=X(107392,191640)&subset=Y(779125,862413)&subset=time(%221985-01-01T00:00:00%22)&format=image/tiff&request=GetCoverage

...
======| 100%
Error: [rast] cannot open this file as a SpatRaster:
C:\Users\TPHILI\AppData\Local\Temp\1\RtmpqESuGX/mrlc_Land-Cover-Native_conus_year_data__Land-Cover-Native_conus_year_data_1985-01-01T00_00_00_107392,779125,191640,862413.tif
In addition: Warning message:
`C:\Users\TPHILI\AppData\Local\Temp\1\RtmpqESuGX/mrlc_Land-Cover-Native_conus_year_data__Land-Cover-Native_conus_year_data_1985-01-01T00_00_00_107392,779125,191640,862413.tif' not recognized as being in a supported file format. (GDAL error 4)

The returned tiff file is identical to the returned xml file:

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">
<ows:Exception exceptionCode="NoApplicableCode">
ows:ExceptionTextFailed to read the coverage mrlc_Land-Cover-Native_conus_year_data__Land-Cover-Native_conus_year_data
The coordinate reference system must be the same for all objects.</ows:ExceptionText>
</ows:Exception>
</ows:ExceptionReport>

• Is there something obviously wrong in my ows4R code?
• Is the crs exception because something else has a different crs (if so, how do I find out what it is?), or is the time dimension somehow triggering it?
• Is there a different parameter to getCoverageStack() that I need to be using? I tried and failed to find the source for getCoverageStack in this repo, which suggests I might be in over my head.
• Should I be asking the USGS MRLC folks behind this data service some specific questions about their WCS services?

After all this, I hope this can serve as another example or case study for the WCS vignette.

R version 4.4.2 (2024-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.utf8 LC_CTYPE=English_United States.utf8
[3] LC_MONETARY=English_United States.utf8 LC_NUMERIC=C
[5] LC_TIME=English_United States.utf8

time zone: America/Los_Angeles
tzcode source: internal

attached base packages:
[1] stats graphics grDevices utils datasets methods base

other attached packages:
[1] sf_1.0-19 ows4R_0.4 terra_1.7-83

loaded via a namespace (and not attached):
[1] geometa_0.8-0 jsonlite_1.8.9 compiler_4.4.2 crayon_1.5.3
[5] Rcpp_1.0.13-1 parallel_4.4.2 assertthat_0.2.1 mime_0.12
[9] readr_2.1.5 R6_2.5.1 curl_6.0.1 classInt_0.4-10
[13] XML_3.99-0.17 tibble_3.2.1 units_0.8-5 openssl_2.2.2
[17] DBI_1.2.3 pillar_1.9.0 tzdb_0.4.0 rlang_1.1.4
[21] utf8_1.2.4 cli_3.6.3 magrittr_2.0.3 class_7.3-22
[25] grid_4.4.2 keyring_1.3.2 fortunes_1.5-4 hms_1.1.3
[29] askpass_1.2.1 lifecycle_1.0.4 vctrs_0.6.5 KernSmooth_2.23-24
[33] proxy_0.4-27 glue_1.8.0 codetools_0.2-20 fansi_1.0.6
[37] e1071_1.7-16 httr_1.4.7 tools_4.4.2 pkgconfig_2.0.3

[tphilippi]

This may be (distantly) related to issue #133 Unsupported file format error in get coverage, although in my case the unsupported file format was that the xml exception report was also returned as the .tif fie.
However, note that getDimensions() works in my example, which suggests I am correct in using 2.0.0 or 2.0.1. And, describeCoverage() works for me, too.

[eblondel]

I had a look to your issue, but still no solution. It would be interesting to get some feedback from the service provider on the exception got "The coordinate reference system must be the same for all objects.".

There is no additional parameter for getCoverageStack . The latter iterates over getCoverage calls (for different time instants here) and builds a timeseries stack on it.

I keep investigating on my side and let you know if I find something

[tphilippi]

Thanks! I sent a detailed query to the usgs.gov customer service email address and the specific usgs contractor named in the metadata. These data and services are new enough that I expect that she is still working there (and perhaps actively adding more datasets to their data services), and will respond.

Thank you for having your functions display the full composed url of the GET or POST! That allowed me to provide exactly what worked (getCapabilities and DescribeCoverage, and getCoverage for a different one of their services), and what threw the xml exception (which google didn't find anything about). I also manually tested by removing the subsetting on time &subset=time(%221985-01-01T00:00:00%22), and that failed with the exact same exception.

As soon as I hear anything back, I will test it and update here.

[eblondel]

Hi @tphilippi did you get some answer from the service provider?