GDAL Translate WCS – Unable to Export GeoTIFF Files with Zero Bands Using gdal_translate

gdalgdal-translatewcs

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)

Best Answer

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.

Related Solutions

[GIS] Fetching GeoTiff data from WCS via GDAL in code

GeoServer includes a helpful WCS request builder (on the demos page) that produced the following XML for me:

Version 1.0.0

<?xml version="1.0" encoding="UTF-8"?>
<GetCoverage version="1.0.0" service="WCS" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns="http://www.opengis.net/wcs" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:gml="http://www.opengis.net/gml" xmlns:ogc="http://www.opengis.net/ogc" xsi:schemaLocation="http://www.opengis.net/wcs http://schemas.opengis.net/wcs/1.0.0/getCoverage.xsd">
  <sourceCoverage>training:output</sourceCoverage>
  <domainSubset>
    <spatialSubset>
      <gml:Envelope srsName="EPSG:3006">
        <gml:pos>372500.0 6164500.0</gml:pos>
        <gml:pos>373000.0 6165000.0</gml:pos>
      </gml:Envelope>
      <gml:Grid dimension="2">
        <gml:limits>
          <gml:GridEnvelope>
            <gml:low>0 0</gml:low>
            <gml:high>0 0</gml:high>
          </gml:GridEnvelope>
        </gml:limits>
        <gml:axisName>x</gml:axisName>
        <gml:axisName>y</gml:axisName>
      </gml:Grid>
    </spatialSubset>
  </domainSubset>
  <output>
    <crs>EPSG:3006</crs>
    <format>GeoTIFF</format>
  </output>
</GetCoverage>

Version 1.1.0

<?xml version="1.0" encoding="UTF-8"?>
<GetCoverage version="1.1.1" service="WCS" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns="http://www.opengis.net/wcs/1.1.1" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:gml="http://www.opengis.net/gml" xmlns:ogc="http://www.opengis.net/ogc" xsi:schemaLocation="http://www.opengis.net/wcs/1.1.1 http://schemas.opengis.net/wcs/1.1.1/wcsAll.xsd">
  <ows:Identifier>training:output</ows:Identifier>
  <DomainSubset>
    <ows:BoundingBox crs="urn:ogc:def:crs:EPSG::3006">
      <ows:LowerCorner>6164500.000000001 372499.99999999994</ows:LowerCorner>
      <ows:UpperCorner>6165000.000000001 372999.99999999994</ows:UpperCorner>
    </ows:BoundingBox>
  </DomainSubset>
  <Output store="true" format="image/tiff">
    <GridCRS>
      <GridBaseCRS>urn:ogc:def:crs:EPSG::3006</GridBaseCRS>
      <GridType>urn:ogc:def:method:WCS:1.1:2dSimpleGrid</GridType>
      <GridOffsets>NaN NaN</GridOffsets>
      <GridCS>urn:ogc:def:cs:OGC:0.0:Grid2dSquareCS</GridCS>
    </GridCRS>
  </Output>
</GetCoverage>

In both cases you would want to POST these files to the server.

If you want to go down the GET route then the best place to look is the WCS spec (page 40) where you will find a table with all the options. You are required to use SERVICE=WCS, VERSION=1.0.0, REQUEST=GetCoverage, COVERAGE=name,CRS=crs_identifier,BBOX=minx,miny,maxx,maxy or TIME=time1,time2, WIDTH=w,HEIGHT=h or RESX,RESY, FORMAT=format

Best Answer

Related Solutions

[GIS] Fetching GeoTiff data from WCS via GDAL in code

Related Question