CRS issue when downloading WCS SoilGrid ; EPSG code not found

[ecor]

Dear Emmanuel @eblondel ,

thanks for your package. I tried to download Soligurd 250 m dataset (WCS service) https://maps.isric.org/. See R reprodubile code below. I found out the following error:

> cov$getDimensions()
[ows4R][INFO] WCSCoverageSummary - Fetching Coverage envelope dimensions by CRS interpretation 
failed to load HTTP resource
Error : 1: failed to load HTTP resource

[ows4R][ERROR] WCSCoverageSummary - Error during CRS interpretation for srsName = 'http://www.opengis.net/def/crs/EPSG/0/152160' 
NULL

> out <- cov$getCoverage()
[ows4R][INFO] WCSCoverageSummary - Fetching Coverage envelope dimensions by CRS interpretation 
failed to load HTTP resource
Error : 1: failed to load HTTP resource

[ows4R][ERROR] WCSCoverageSummary - Error during CRS interpretation for srsName = 'http://www.opengis.net/def/crs/EPSG/0/152160' 
<GMLEnvelope>
....|-- lowerCorner: -19949750 -6147500
....|-- upperCorner: 19861750 8361000[ows4R][INFO] WCSCoverageSummary - Fetching Coverage envelope dimensions by CRS interpretation 
failed to load HTTP resource
Error : 1: failed to load HTTP resource

[ows4R][ERROR] WCSCoverageSummary - Error during CRS interpretation for srsName = 'http://www.opengis.net/def/crs/EPSG/0/152160' 
[ows4R][INFO] WCSGetCoverage - Fetching https://maps.isric.org/mapserv?map=/map/clay.map&service=WCS&version=2.0.0&coverageId=clay_0-5cm_mean&subset=x(-19949750,19861750)&subset=y(-6147500,8361000)&format=image/tiff&request=GetCoverage 
Downloading: 600 B     [ows4R][ERROR] WCSGetCoverage - Error while executing request 'GetCoverage' 
Error during wrapup: [rast] cannot open this file as a SpatRaster: /tmp/RtmpRqVvk8/clay_0-5cm_mean.tif
Error: no more error handlers available (recursive errors?); invoking 'abort' restart

It look like that the projections EPSG 152160 was created by ESRI and not recognized (see https://www.isric.org/explore/soilgrids/faq-soilgrids#How_can_I_use_SoilGrids_in_a_different_projection) . Is there a way to insert manually the CRS / SRDS ?

Thank you best, Emanuele @ecor

library(ows4R)
library(terra)

#get WCS client connection (be careful, it's WCS, not CSW)
#WCS = Web Coverage Service (for raster data access)
#CSW = Catalogue Service for the Web (for geographic metadata cataloguing)

## https://maps.isric.org/ 

## CLAY MAP
WCS = WCSClient$new(
  url = "https://maps.isric.org/mapserv?map=/map/clay.map",
  serviceVersion = "2.0.0",
  logger = "INFO"
)

#find the coverage
cov = WCS$capabilities$findCoverageSummaryById("clay_0-5cm_mean")
cov$getDescription()
#get data dimensions (here's a spatial 2D coverage)

## https://www.isric.org/explore/soilgrids/faq-soilgrids#How_can_I_use_SoilGrids_in_a_different_projection
cov$getDimensions()

out <- cov$getCoverage()

[eblondel]

I see, indeed it deals with an ESRI projection for which a pseudo EPSG code was given, the reason why getting http://www.opengis.net/def/crs/EPSG/0/152160 fails, there is no GML CR definition associated to it (see on contrary this example https://www.opengis.net/def/crs/EPSG/0/4326). Indeed I make a non-silent try call so you see an Error appearing. The getDimensions should not be a blocker for getting the coverage, especially here it's a 2D grid, and we already get properly the GML envelope associated with the coverage.

However I see that getting coverage fails, and it's due to the default grid envelope that is not accepted. You can see the message by browsing this request on the web: https://maps.isric.org/mapserv?map=/map/clay.map&service=WCS&version=2.0.0&coverageId=clay_0-5cm_mean&subset=x(-19949750,19861750)&subset=y(-6147500,8361000)&format=image/tiff&request=GetCoverage In ows4R, service exceptions are not handled, so user is not well informed about the problem. I will support this in #116

[eblondel]

@ecor For the get coverage, I've improved the handling of service exceptions in #116, and also took the opportunity to support a custom print function for the ows4R objects. Consequently, when trying to download the full coverage above, you will get this:

<OWSException>
....|-- ExceptionText: msWCSGetCoverage20(): WCS server error. Raster size out of range, width and height of resulting coverage must be no more than MAXSIZE=16384.

These exceptions were hidden before, and the default print didn't help in understanding.

Let me know

[ecor]

Thanks @eblondel ,

Here is the output I obtained, afrer installing latest version of ows4R package from this Github Repo (remotes::install_github("eblondel/ows4R")) :

> library(ows4R)
Loading ISO 19139 XML schemas...
Loading ISO 19115 codelists...
> library(terra)
terra 1.7.69
> 
> 
> #get WCS client connection (be careful, it's WCS, not CSW)
> #WCS = Web Coverage Service (for raster data access)
> #CSW = Catalogue Service for the Web (for geographic metadata cataloguing)
> 
> ## https://maps.isric.org/ 
> 
> ## CLAY MAP
> WCS = WCSClient$new(
+     url = "https://maps.isric.org/mapserv?map=/map/clay.map",
+     serviceVersion = "2.0.0",
+     logger = "INFO"
+ )
[ows4R][INFO] OWSGetCapabilities - Fetching https://maps.isric.org/mapserv?map=/map/clay.map&service=WCS&version=2.0.0&request=GetCapabilities 
  |======================================================================================================================| 100%
> 
> #find the coverage
> cov = WCS$capabilities$findCoverageSummaryById("clay_0-5cm_mean")
> cov$getDescription()
[ows4R][INFO] WCSDescribeCoverage - Fetching https://maps.isric.org/mapserv?map=/map/clay.map&service=WCS&version=2.0.0&coverageId=clay_0-5cm_mean&request=DescribeCoverage 
  |======================================================================================================================| 100%
Loading required package: sf
Linking to GEOS 3.10.2, GDAL 3.4.3, PROJ 8.2.0; sf_use_s2() is FALSE
<WCSCoverageDescription>
....|-- boundedBy [srsName=http://www.opengis.net/def/crs/EPSG/0/152160,axisLabels=x y,uomLabels=m m,srsDimension=2] <GMLEnvelope>
........|-- lowerCorner: -19949750 -6147500
........|-- upperCorner: 19861750 8361000
....|-- domainSet [dimension=2,gml:id=grid_clay_0-5cm_mean] <GMLRectifiedGrid>
........|-- limits  <GMLGridEnvelope>
............|-- low  
................|-- value: 0 0
............|-- high  
................|-- value: 159245 58033
........|-- axisLabels  
............|-- value: x y
........|-- origin [gml:id=grid_origin_clay_0-5cm_mean,srsName=http://www.opengis.net/def/crs/EPSG/0/152160] <GMLPoint>
............|-- pos: -19949625 8360875
........|-- offsetVector [srsName=http://www.opengis.net/def/crs/EPSG/0/152160]
............|-- value: 250.000000 0
........|-- offsetVector [srsName=http://www.opengis.net/def/crs/EPSG/0/152160]
............|-- value: 0 -250.000000
....|-- rangeType  <SWEDataRecord>
........|-- field <SWEQuantity>
............|-- nilValues  
............|-- uom [code=W.m-2.Sr-1] 
............|-- constraint  
................|-- AllowedValues  
....................|-- interval  
........................|-- value: 0 65535
....................|-- significantFigures  
........................|-- value: 5
....|-- CoverageId  
........|-- value: clay_0-5cm_mean
....|-- ServiceParameters <ISOElementSequence>
........|-- CoverageSubtype: RectifiedGridCoverage
........|-- nativeFormat: image/tiff> #get data dimensions (here's a spatial 2D coverage)
> 
> ## https://www.isric.org/explore/soilgrids/faq-soilgrids#How_can_I_use_SoilGrids_in_a_different_projection
> cov$getDimensions()
[ows4R][INFO] WCSCoverageSummary - Fetching Coverage envelope dimensions by CRS interpretation 
No encoding supplied: defaulting to UTF-8.
[ows4R][ERROR] WCSCoverageSummary - Error during CRS interpretation for srsName = 'https://www.opengis.net/def/crs/EPSG/0/152160' 
NULL
> 
> out <- cov$getCoverage()
[ows4R][INFO] WCSCoverageSummary - Fetching Coverage envelope dimensions by CRS interpretation 
No encoding supplied: defaulting to UTF-8.
[ows4R][ERROR] WCSCoverageSummary - Error during CRS interpretation for srsName = 'https://www.opengis.net/def/crs/EPSG/0/152160' 
<GMLEnvelope>
....|-- lowerCorner: -19949750 -6147500
....|-- upperCorner: 19861750 8361000[ows4R][INFO] WCSCoverageSummary - Fetching Coverage envelope dimensions by CRS interpretation 
No encoding supplied: defaulting to UTF-8.
[ows4R][ERROR] WCSCoverageSummary - Error during CRS interpretation for srsName = 'https://www.opengis.net/def/crs/EPSG/0/152160' 
[ows4R][INFO] WCSGetCoverage - Fetching https://maps.isric.org/mapserv?map=/map/clay.map&service=WCS&version=2.0.0&coverageId=clay_0-5cm_mean&subset=x(-19949750,19861750)&subset=y(-6147500,8361000)&format=image/tiff&request=GetCoverage 
Downloading: 600 B     
> cov$getDimensions()
[ows4R][INFO] WCSCoverageSummary - Fetching Coverage envelope dimensions by CRS interpretation 
No encoding supplied: defaulting to UTF-8.
[ows4R][ERROR] WCSCoverageSummary - Error during CRS interpretation for srsName = 'https://www.opengis.net/def/crs/EPSG/0/152160' 
NULL
> out
<OWSException>
....|-- ExceptionText: msWCSGetCoverage20(): WCS server error. Raster size out of range, width and height of resulting coverage must be no more than MAXSIZE=16384.
> 

[eblondel]

@ecor this is expected behavior:

• for the CRS EPSG2160, we still get an error, because there is no official EPSG code for this, so 'https://www.opengis.net/def/crs/EPSG/0/152160' resource doesn't exist, but in your case it only prevents to inherit proper CRS definitions associated to the dimensions.
• for getting the coverage as terra object, you need see explicitely the exception, in your case it says the raster size is out of range. This is set by the server to prevent people getting the full coverage. You will have to set a filter to reduce the bounding box of the area you want.

[ecor]

Thank you @eblondel for the clarification. Now it works, setting a correct BBOX: For me ticket can be closed.

library(ows4R)
library(terra)

#get WCS client connection (be careful, it's WCS, not CSW)
#WCS = Web Coverage Service (for raster data access)
#CSW = Catalogue Service for the Web (for geographic metadata cataloguing)

## https://maps.isric.org/ 

## CLAY MAP
WCS = WCSClient$new(
  url = "https://maps.isric.org/mapserv?map=/map/clay.map",
  serviceVersion = "2.0.0",
  logger = "INFO"
)

#find the coverage
cov = WCS$capabilities$findCoverageSummaryById("clay_0-5cm_mean")
cov$getDescription()
#get data dimensions (here's a spatial 2D coverage)

## https://www.isric.org/explore/soilgrids/faq-soilgrids#How_can_I_use_SoilGrids_in_a_different_projection
cov$getDimensions()

 ## https://git.wur.nl/isric/soilgrids/soilgrids.notebooks/-/blob/master/markdown/webdav_from_R.md
bb=c(-337500.000,1242500.000,152500.000,527500.000) # Example bounding box (homolosine) for Ghana

out <- cov$getCoverage(OWSUtils$toBBOX(-437500.000,2242500.000,152500.000,527500.000))

[eblondel]

Thanks @ecor for your feedback, great to see it works!