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

cgdalwcs

UPDATE: We seem to have found where the trouble lies, but sadly it is a problem that is built into GDAL's handling of WCS requests and it doesn't seem like we can bypass it. Instead we are going to write our own WCS client to fetch the GeoTiff data and then load in into GDAL from memory. See "EDIT 4" at the bottom for details.

The trouble I'm having is in writing the xml string which I send to the server. I either get an error message from the server when I try to fetch the Dataset src or it passes but I've apparently fetched all the data on the server (which can be seen on the size of the raster) despite wanting just a piece of it.

using(OSGeo.GDAL.Dataset src = OSGeo.GDAL.Gdal.Open(xml, OSGeo.GDAL.Access.GA_ReadOnly))
{
    int rasters = src.RasterCount;
    for(int i = 1; i <= rasters; ++i)
    {
        using(OSGeo.GDAL.Band raster = src.GetRasterBand(i))
        {
            //Yet to be written code
        }
    }
}

I've been unable to find a list of which attributes that can be sent to the server. The attributes I've tried are all attributes which I've seen other people try while I've been googling for help.

Some attributes I've seen different alternatives for, and I've tried each. For projection I've tried SRS, CRS and Projection. For image size I've tried WIDTH and HEIGHT and also SizeX and SizeY. For bounding box I've tried BoundingBox and BBOX.

This is one way that I've attempted to construct the xml string:

string xml = "<WCS_GDAL>";
xml += "<ServiceURL>http://example.com/geoserver/wcs?</ServiceURL>";
xml += "<VERSION>1.1.0</VERSION>";*
xml += "<REQUEST>GetCoverage</REQUEST>";
xml += "<CoverageName>coverage_name</CoverageName>";
xml += "<SRS>EPSG:3006</SRS>";
xml += "<CRS>EPSG:3006</CRS>";
xml += "<Projection>EPSG:3006</Projection>";
xml += "<GetCoverageExtra>";
{
    xml += "&amp;CRS=EPSG:3006";
    xml += "&amp;BBOX=372500.0,6164500.0,373000.0,6165000.0" + ",EPSG:3006";**
    xml += "&amp;WIDTH=200&amp;HEIGHT=200";
    xml += "&amp;FORMAT=GTiff";
}
xml += "</GetCoverageExtra>";
xml += "</WCS_GDAL>";

*I've also tried 1.0.0.
**1.1.0 requires a projection to be part of the bounding box while for 1.0.0 GDAL doesn't allow more than 4 entries in the list.

This is another way that I tried to do it. In this case by affecting the final ServiceURL directly. GDAL will automatically append version and request type before the first ampersand, so I had to include them in the XML so that it doesn't append the wrong ones. GDAL will also cast an exception if CoverageName is not included in the XML.

string link = "http://example.com/geoserver/wcs?";
//link += "&VERSION=1.1.0";
//link += "&REQUEST=GetCoverage";
link += "&coverage=coverage_name";
link += "&CRS=EPSG:3006";
link += "&BBOX=372500.0,6164500.0,373000.0,6165000.0,EPSG:3006";
link += "&WIDTH=200";
link += "&HEIGHT=200";
link += "&FORMAT=GTiff";

string xml = "<WCS_GDAL>";
xml += "<ServiceURL>";
xml += link;
xml += "</ServiceURL>";
xml += "<VERSION>1.1.0</VERSION>";*
xml += "<REQUEST>GetCoverage</REQUEST>";
xml += "<CoverageName>coverage_name</CoverageName>";
xml += "</WCS_GDAL>";

If "<REQUEST>GetCoverage</REQUEST>" isn't included it seems that GDAL defaults to DescribeCoverage. (EDIT: This is incorrect. GDAL doesn't allow you to do DescribeCoverage. GDAL automatically does DescribeCoverage before it does GetCoverage.) And if version is not included GDAL seems to default to 1.0.0.

Depending on what attributes I include, what values I've given them and what version I specify I get a few different Java exceptions from the server. Sometimes it says java.util.ArrayList cannot be cast to org.geotools.geometry.GeneralEnvelope (with version 1.0.0) and sometimes it says java.lang.ClassCastException: [Ljava.lang.String; cannot be cast to java.lang.String.

The exception that I think managed to go the furthest and which I can't seem to get past is org.geotools.coverage.grid.InvalidGridGeometryException: Unspecified image's size..

When I use the version that directly affects the link I can bypass the problem of GDAL not allowing 5 entries for the bounding box with version 1.0.0, which makes the server not cast the ArrayList to GeneralEnvelope cast exception. However, it seems that it then completely ignores the bounding box and returns all the data on the server. If I switch version to 1.1.0 then I get the unspecified image's size exception.

There should be a version 2.0.0 of WCS, but considering the problems I'm having and the fact that I have yet to see any examples of someone else using it has made me unwilling to try it right now. EDIT: Scratch that, GDAL doesn't seem to support it… At least that's what the exceptions tells me.

In short: I need help with finding a way to fetch a GeoTiff of a small area via WCS.

EDIT: This is the XML which I'm currently trying to use, as asked for by iant:

<WCS_GDAL>
    <ServiceURL>http://example.com/geoserver/wcs?</ServiceURL>
    <SERVICE>WCS</SERVICE>
    <VERSION>1.1.0</VERSION>
    <REQUEST>GetCoverage</REQUEST>
    <CoverageName>coverage_name</CoverageName>
    <CRS>urn:ogc:def:crs:EPSG:3006</CRS>
    <GetCoverageExtra>&BoundingBox=372500.0,6164500.0,373000.0,6165000.0,urn:ogc:def:crs:EPSG:3006&WIDTH=200&HEIGHT=200&RESX=200&RESY=200&SizeX=200&SizeY=200&FORMAT=image/GeoTiff</GetCoverageExtra>
</WCS_GDAL>

As you can see, I've resorted to using a lot of redundancy in hope that one of the attributes will be right. This XML isn't written like this btw. I do as I showed above where I build them in a string. Though some of the attributes and values are different since I wrote that.

EDIT 2: iant actually asked for the generated URL string. Here it is:
/geoserver/wcs?SERVICE=WCS&VERSION=1.1.0&REQUEST=GetCoverage&IDENTIFIER=coverage_name&FORMAT=image%2Ftiff&BOUNDINGBOX=200025,7699925,200075,7699975,urn:ogc:def:crs:EPSG::3006&RangeSubset=contents

It seems like GDAL has messed up the URL big time. It's supposed to include everything in the GetCoverageExtra tags straight of, from what I hear. Yet it doesn't have the WIDTH & HEIGHT… It also messed up the values in the bounding box for some reason.

EDIT 3: After URL-escaping the ampersands in the GetCoverageExtra tag the URL is more complete:
/geoserver/wcs?SERVICE=WCS&VERSION=1.1.0&REQUEST=GetCoverage&IDENTIFIER=swe%3Anh_riks_Sweref_99_TM_geotiff&FORMAT=image%2Ftiff&BOUNDINGBOX=200025,7699925,200075,7699975,urn:ogc:def:crs:EPSG::3006&RangeSubset=contents&BOUNDINGBOX=372500.0,6164500.0,373000.0,6165000.0,urn:ogc:def:crs:EPSG:3006&WIDTH=200
Something is still off though. There's two bounding boxes! Both the one which is wrong, which seems like some kind of default value, and the one which I wrote. WIDTH is included, but not HEIGHT, however I found that this was my own fault because I screwed up the ampersand before HEIGHT. After fixing it I get the following exception from the server: java.util.ArrayList cannot be cast to net.opengis.ows11.BoundingBoxType. The default value seems to have double colons between EPSG and 3006, but even after adding another colon I get the same error. I also tried to remove the ".0" from the bounding box values.

EDIT 4 (final?): By looking at the GDAL source code we seem to have found what GDAL is doing wrong. When opening a Dataset via WCS GDAL first automatically sends a DescribeCoverage to the server and uses the returned data to set a bounding box. This means that the bounding box is set to cover all the data in the coverage I'm asking for. This bounding box is then sent with the GetCoverage request regardless of what data you try to give GDAL and the server seems to be unable to handle receiving two bounding boxes in a single request. So there seems to be no way of getting a specified area in WCS via GDAL. My colleague is writing a separate class to fetch specific areas from the server without GDAL and then load that data into GDAL from memory. It seems to be the only way for now, other than compiling our own version of GDAL.

Best Answer

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

Related Solutions

[GIS] How to automate downloading “slices” of data from Web Coverage Service (WCS)

WCS is the less-used of the WMS/WFS/WCS triangle so there isn't as much out there as you've probably noticed.

However, making requests is really really easy and you don't really need a special library. Arguably its the easiest of the three from a users perspective because you get a nice simple file back.

A valid request looks something like this:

http://my.host.com/cgi-bin/mapserv?map=mywcs.map&SERVICE=wcs &VERSION=1.0.0&REQUEST=GetCoverage&coverage=bathymetry&CRS=EPSG:42304&OX=-2200000,-712631,3072800,3840000&WIDTH=3199&HEIGHT=2833&FORMAT=GTiff

It may look complex, but broken down its very simple:

http://my.host.com/cgi-bin/mapserv?map=mywcs.map
&SERVICE=wcs
&VERSION=1.0.0
&REQUEST=GetCoverage
&coverage=bathymetry
&CRS=EPSG:42304
&BBOX=-2200000,-712631,3072800,3840000
&WIDTH=3199
&HEIGHT=2833
&FORMAT=GTiff

As you can see, most of them are self explanatory. So you just need to concatenate a string that has those variables. Most of them will probably be constants.

Coverage is the name of the layer you're requesting.

BBOX is the hardest part to do as you need to get the coordinates from somewhere.

I understand the point of libraries, but when something's this simple, "re-inventing the wheel" is actually a sensible option. You do of course need some sort of HTTP library.

A search for WCS GetCoverage finds other examples. This page has a live example: http://www.mapserver.org/ogc/wcs_server.html#test-your-wcs-1-0-server

[GIS] Qgis Server and WCS service : problem to display raster layers

My problem is solved ! Not as i wanted but i have obtained the result i needed.

So, first, it's not really a bug in QGIS server.According severals sources, the reason why you can't display raster (based on a ECW file) layer in qgis server is due to a paid license : ERDAS ECW/JP2 SDK you can see here : download.intergraph.com/?ProductName=ERDAS%20ECW/JPEG2000%20SDK

You have to know that i haven't tested this solution because it's sounded too tedious for what i want (and secondarily i don't want to pay something without knowing if it will work).

But, if you want to try, here is some topic i found about this :

the problem : www.forumsig.org/showthread.php/37886-Probl%C3%A8mes-QGIS-Server
how to install the library : doc.ubuntu-fr.org/qgis
how to install the library : trac.osgeo.org/gdal/wiki/ECW

So ! This is not the solution a choosed. I think you get that the problem is about ECW files and not about the entire rasters layers. Starting from here, i simply used tif file instead of ecw as a raster and everything goes well !

Two solutions to get tif file :

-> you can convert ECW file into tif with QGIS Desktop (Menu -> Raster -> Conversion -> Convertir) : it's a very long process if you have a heavy ECW file. For example, it takes me 10mn to convert a 900Mo ECW to a 10GB tif

-> you can download tif map directly online (Example : professionnels.ign.fr/catalogue )

I realise i'm not so clear about the first point (about paid license etc). It's due to i don't really understand what is this lib...

Hope this will help some people.

Bye !

Best Answer

Related Solutions

[GIS] How to automate downloading “slices” of data from Web Coverage Service (WCS)

[GIS] Qgis Server and WCS service : problem to display raster layers

Related Question