1

In accordance with GDAL's WCS examples, I am trying to use gdal_translate to import a GTiff raster from a coverage.

The server in question is a GeoServer instance of the German Weather Service (DWD) at https://cdc.dwd.de/geoserver/.

The coverage ID is CDC__GRD_DEU_P1D_BF-GRB. Mind the double underscore. Cf. GetCapabilities URL.

Running gdalinfo:

gdalinfo \
"WCS:https://cdc.dwd.de/geoserver/ows?version=2.0.1&coverage=CDC__GRD_DEU_P1D_BF-GRB" \
--debug ON

against it gives me a valid response (see bottom of text).

I understand the layer has two additional dimensions: time and elevation. As per the examples, I am trying to subset the response by those two dimensions using:

gdal_translate \
-oo Subset="time(2022-02-20T00:00:00.000Z)" \
-oo Subset="elevation(10.0)" \
-oo CACHE=wcs_cache \
-oo CLEAR_CACHE \
-oo INTERLEAVE=PIXEL \
"WCS:https://cdc.dwd.de/geoserver/ows?version=2.0.1&coverage=CDC__GRD_DEU_P1D_BF-GRB" \
out.tif \
--debug ON

However, I get the following response:

HTTP: Fetch(https://cdc.dwd.de/geoserver/ows?version=2.0.1&SERVICE=WCS&REQUEST=GetCapabilities)
HTTP: libcurl/7.80.0 OpenSSL/1.1.1l zlib/1.2.11 brotli/1.0.9 nghttp2/1.46.0
WCS: Requesting https://cdc.dwd.de/geoserver/ows?SERVICE=WCS&REQUEST=DescribeCoverage&VERSION=2.0.1&COVERAGEID=CDC__GRD_DEU_P1D_BF-GRB&FORMAT=text/xml
HTTP: Fetch(https://cdc.dwd.de/geoserver/ows?SERVICE=WCS&REQUEST=DescribeCoverage&VERSION=2.0.1&COVERAGEID=CDC__GRD_DEU_P1D_BF-GRB&FORMAT=text/xml)
GDAL: GDALOpen(WCS:https://cdc.dwd.de/geoserver/ows?version=2.0.1&coverage=CDC__GRD_DEU_P1D_BF-GRB, this=0x7fdc19e40800) succeeds as WCS.
GDAL: Using GTiff driver
Input file size is 866, 654
GDAL: Using GTiff driver
ERROR 1: out.tif: Unable to export GeoTIFF files with zero bands.
GDAL: GDALClose(wcs_cache/uZCYr.xml, this=0x7fdc19e40800)
GDAL: In GDALDestroy - unloading GDAL shared library.

As per this gdal-dev mailing list thread I would assume that I'm not using the Subset parameter correctly. However, e.g. using different combinations of time and elevation formatting, e.g. 10 instead of 10.0 did not succeed.

Debug hint: For the time dimension, the coverage does not not contain a raster for the current day, or some days prior. ~10 days prior is a safe bet.

I am wondering what I'm missing.



gdalinfo output:

/ # gdalinfo \
> "WCS:https://cdc.dwd.de/geoserver/ows?version=2.0.1&coverage=CDC__GRD_DEU_P1D_BF-GRB" \
> --debug ON
HTTP: Fetch(https://cdc.dwd.de/geoserver/ows?version=2.0.1&SERVICE=WCS&REQUEST=GetCapabilities)
HTTP: libcurl/7.80.0 OpenSSL/1.1.1l zlib/1.2.11 brotli/1.0.9 nghttp2/1.46.0

WCS: Requesting https://cdc.dwd.de/geoserver/ows?SERVICE=WCS&REQUEST=DescribeCoverage&VERSION=2.0.1&COVERAGEID=CDC__GRD_DEU_P1D_BF-GRB&FORMAT=text/xml
HTTP: Fetch(https://cdc.dwd.de/geoserver/ows?SERVICE=WCS&REQUEST=DescribeCoverage&VERSION=2.0.1&COVERAGEID=CDC__GRD_DEU_P1D_BF-GRB&FORMAT=text/xml)
GDAL: GDALOpen(WCS:https://cdc.dwd.de/geoserver/ows?version=2.0.1&coverage=CDC__GRD_DEU_P1D_BF-GRB, this=0x7f2fd9faf000) succeeds as WCS.
Driver: WCS/OGC Web Coverage Service
Files: /root/.gdal/wcs_cache/UUhpN.xml
       /root/.gdal/wcs_cache/UUhpN.xml.aux.xml
Size is 866, 654
Coordinate System is:
PROJCRS["DHDN / 3-degree Gauss-Kruger zone 3",
    BASEGEOGCRS["DHDN",
        DATUM["Deutsches Hauptdreiecksnetz",
            ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4314]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",9,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",3500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    ID["EPSG",31467]]
Data axis to CRS axis mapping: 2,1
GeoTransform =
  3281415, 0, -1000
  6102501, 1000, 0
Metadata:
  Abstract=Generated from ImageMosaic
  CoverageSubtype=RectifiedGridCoverage
  DIMENSION_0_AXIS=X
  DIMENSION_0_INTERVAL=3280415,3934415
  DIMENSION_0_UOM=m
  DIMENSION_1_AXIS=Y
  DIMENSION_1_INTERVAL=5237501,6103501
  DIMENSION_1_UOM=m
  DIMENSION_2_AXIS=elevation
  DIMENSION_2_INTERVAL=10.0,200.0
  DIMENSION_2_UOM=m
  DIMENSION_3_AXIS=time
  DIMENSION_3_UOM=s
  DOMAIN=X,Y
  FIELD_1_DESCR=nFK
  FIELD_1_INTERVAL=-32768 32767
  FIELD_1_NAME=nFK
  Keywords=GRD_DEU_P1D_BF-GRB,WCS,ImageMosaic
  TimeDomain=
  Title=Tägliche Raster der mittleren Bodenfeuchte unter Gras für Deutschland
  WCS_GLOBAL#crsSupported=http://www.opengis.net/def/crs/EPSG/0/2000,http://www.opengis.net/def/crs/EPSG/0/2001,

[...]

  WCS_GLOBAL#formatSupported=application/gml+xml,application/x-gzip,application/x-netcdf,application/x-netcdf4,image/jpeg,image/png,image/tiff,text/plain
  WCS_GLOBAL#InterpolationSupported=nearest-neighbor,linear,cubic
  WCS_GLOBAL#Profiles=WCS/2.0/conf/core,WCS_protocol-binding_get-kvp/1.0.1,WCS_protocol-binding_post-xml/1.0,WCS_service-extension_crs/1.0/conf/crs-gridded-coverage,WCS_geotiff-coverages/1.0/conf/geotiff-coverage,GMLCOV/1.0/conf/gml-coverage,GMLCOV/1.0/conf/special-format,GMLCOV/1.0/conf/multipart,WCS_service-extension_scaling/1.0/conf/scaling,WCS_service-extension_crs/1.0/conf/crs,WCS_service-extension_interpolation/1.0/conf/interpolation,WCS_service-extension_interpolation/1.0/conf/interpolation-per-axis,WCS_service-extension_interpolation/1.0/conf/nearest-neighbor,WCS_service-extension_interpolation/1.0/conf/linear,WCS_service-extension_interpolation/1.0/conf/cubic,WCS_service-extension_range-subsetting/1.0/conf/record-subsetting
  WCS_GLOBAL#ServiceIdentification.AccessConstraints=NONE
  WCS_GLOBAL#ServiceIdentification.Fees=NONE
  WCS_GLOBAL#ServiceProvider.ProviderName=Deutscher Wetterdienst
  WCS_GLOBAL#ServiceProvider.ServiceContact.ContactInfo.Address.AdministrativeArea=Hessen
  WCS_GLOBAL#ServiceProvider.ServiceContact.ContactInfo.Address.City=Offenbach
  WCS_GLOBAL#ServiceProvider.ServiceContact.ContactInfo.Address.Country=Germany
  WCS_GLOBAL#ServiceProvider.ServiceContact.ContactInfo.Address.DeliveryPoint=Frankfurter Strasse 135
  WCS_GLOBAL#ServiceProvider.ServiceContact.ContactInfo.Address.ElectronicMailAddress=info@dwd.de
  WCS_GLOBAL#ServiceProvider.ServiceContact.ContactInfo.Address.PostalCode=63067
  WCS_GLOBAL#updateSequence=4778
  WCS_GLOBAL#version=2.0.1
Corner Coordinates:
Upper Left  ( 3281415.000, 6102501.000) (  5d34'57.92"E, 55d 0'18.10"N)
Lower Left  ( 2627415.000, 6102501.000) (  4d26'42.32"W, 54d18' 0.81"N)
Upper Right ( 3281415.000, 6968501.000) (  4d43' 5.35"E, 62d45'46.76"N)
Lower Right ( 2627415.000, 6968501.000) (  7d41'10.87"W, 61d48'46.17"N)
Center      ( 2954415.000, 6535501.000) (  0d24'13.74"W, 58d35'50.92"N)
GDAL: GDALClose(/root/.gdal/wcs_cache/UUhpN.xml, this=0x7f2fd9faf000)
Improve this question

1 Answer 1

Reset to default
3

Unfortunately it is hard to say why GDAL fails in finding a band from this WCS coverage. The DescribeCoverage request is https://cdc.dwd.de/geoserver/ows?service=WCS&version=2.0.1&request=DescribeCoverage&coverageID=CDC__GRD_DEU_P1D_BF-GRB

and it returns the range type that looks good to me

<gmlcov:rangeType>
 <swe:DataRecord>
  <swe:field name="nFK">
   <swe:Quantity>
    <swe:description>nFK</swe:description>
     <swe:nilValues>
      <swe:NilValues>
       <swe:nilValue reason="http://www.opengis.net/def/nil/OGC/0/unknown">-9999.0</swe:nilValue>
      </swe:NilValues>
     </swe:nilValues>
     <swe:uom code="‰ nFK"/>
     <swe:constraint>
      <swe:AllowedValues>
       <swe:interval>-32768 32767</swe:interval>
      </swe:AllowedValues>
     </swe:constraint>
   </swe:Quantity>
  </swe:field>
 </swe:DataRecord>
</gmlcov:rangeType>

However, if you look inside the xml file that the WCS driver generates by GetCapabilities and DescribeCoverage that is in your case /root/.gdal/wcs_cache/UUhpN.xml as printed in the beginning of the gdalinfo report, there the band count is reported as zero <BandCount>0</BandCount>.

I fear that for making GDAL to work with this WCS service would require that a developer who knows the WCS driver well would have a look and fix it. Meanwhile you must write the GetCoverage requests manually but with WCS 2.0.1 it is not too difficult. With your subset the request is

https://cdc.dwd.de/geoserver/ows?service=WCS&version=2.0.1&request=GetCoverage&coverageID=CDC__GRD_DEU_P1D_BF-GRB&SUBSET=elevation(10)&SUBSET=time("2022-02-20T00:00:00.000Z")

The elevation subset does not have an effect because 10 is the default <wcsgs:ElevationDomain uom="cm" default="10.0"> so any other value is better for testing. The subset syntax is slightly different than in the open options of the GDAL WCS driver. The syntax above comes from the standard.

Improve this answer
2
  • Thank you for the enlightenment on the band count! I learned something new for the future! My backup idea is indeed to retrieve the raster via GetCoverage. I was only recently introduced to GDAL's interfaces for "external" data and it looked like the more elegant solution. E: Apparently I cannot even upvote your answer, but I will gladly set it as accepted should no working solution with GDAL WCS turn up.
    – user18171348
    Commented Mar 1, 2022 at 12:55
  • Right, GDAL can make some operations easier because it is made to behave in the same way with any data sources. Unfortunately there is some issue with this WCS server. The server that is used in the WCS documentation examples seems to be still up and the subset examples work.
    – user30184
    Commented Mar 1, 2022 at 13:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.