[gdal-dev] GDAL WCS driver

Fri Oct 5 15:32:58 PDT 2012

I have encountered a couple of issues with the WCS driver of GDAL 1.9.0
when accessing my WCS 1.0.0 server.  The WCS driver doesn't recognize
the color interpretations of image bands.  For example (from a server
not accessible to the outside world):

[dev:~] gdalinfo gdal_wcs_terrapixel.xml 
Driver: WCS/OGC Web Coverage Service
Files: gdal_wcs_terrapixel.xml
Size is 106056, 28584
Coordinate System is:
PROJCS["NAD83 / UTM zone 15N",
...
Origin = (260489.000000000000000,3667255.000000000000000)
Pixel Size = (0.500000000000000,-0.500000000000000)
Corner Coordinates:
Upper Left  (  260489.000, 3667255.000) ( 95d34' 1.15"W, 33d 7' 3.30"N)
Lower Left  (  260489.000, 3652963.000) ( 95d33'47.69"W, 32d59'19.66"N)
Upper Right (  313517.000, 3667255.000) ( 94d59'56.38"W, 33d 7'40.78"N)
Lower Right (  313517.000, 3652963.000) ( 94d59'45.90"W, 32d59'56.96"N)
Center      (  287003.000, 3660109.000) ( 95d16'52.90"W, 33d 3'31.34"N)
Band 1 Block=1024x512 Type=Byte, ColorInterp=Undefined
  Overviews: 53028x14292, 26514x7146, 13257x3573, 6628x1786, 3314x893, 1657x446, 828x223
Band 2 Block=1024x512 Type=Byte, ColorInterp=Undefined
  Overviews: 53028x14292, 26514x7146, 13257x3573, 6628x1786, 3314x893, 1657x446, 828x223
Band 3 Block=1024x512 Type=Byte, ColorInterp=Undefined
  Overviews: 53028x14292, 26514x7146, 13257x3573, 6628x1786, 3314x893, 1657x446, 828x223
Band 4 Block=1024x512 Type=Byte, ColorInterp=Undefined
  Overviews: 53028x14292, 26514x7146, 13257x3573, 6628x1786, 3314x893, 1657x446, 828x223

They're all "Undefined".  However, if I use GDAL to examine the 2x2 GeoTIFF
"sample" image that GDAL requests from the WCS server to determine the
band types, GDAL reports the information properly, of course:

[dev:~] gdalinfo terra_gdal_sampler.tif
Driver: GTiff/GeoTIFF
Files: terra_gdal_sampler.tif
Size is 2, 2
Coordinate System is:
PROJCS["NAD83 / UTM zone 15N",
...
Origin = (260488.500000000000000,3667255.500000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_IMAGEDESCRIPTION=Geo.scene_id: terrapixel
  TIFFTAG_RESOLUTIONUNIT=1 (unitless)
  TIFFTAG_SOFTWARE=CubeWerx Suite 6.3.4
  TIFFTAG_XRESOLUTION=1
  TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  260488.500, 3667255.500) ( 95d34' 1.17"W, 33d 7' 3.32"N)
Lower Left  (  260488.500, 3667253.500) ( 95d34' 1.16"W, 33d 7' 3.25"N)
Upper Right (  260490.500, 3667255.500) ( 95d34' 1.09"W, 33d 7' 3.32"N)
Lower Right (  260490.500, 3667253.500) ( 95d34' 1.09"W, 33d 7' 3.25"N)
Center      (  260489.500, 3667254.500) ( 95d34' 1.13"W, 33d 7' 3.29"N)
Band 1 Block=2x2 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA 
Band 2 Block=2x2 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA 
Band 3 Block=2x2 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA 
Band 4 Block=2x2 Type=Byte, ColorInterp=Alpha

Shouldn't it be able to retain this information?  As it is, clients won't
be able to utilize the Alpha channel.  It looks like it's thrown away in
the inforation being cached in the <WCS_GDAL> file:

  <PreferredFormat>GeoTIFF</PreferredFormat>
  <BandCount>4</BandCount>
  <BandType>Byte</BandType>

There's no band description other than the count and data type and
it doesn't attempt to retain the 2x2 sample GeoTIFF image, which is
relatively small.

There's also an oddity with the sample image requested.  The request from
GDAL in my http log is:

base_url?datastore=cwtile&SERVICE=WCS&VERSION=1.0.0&REQUEST=GetCoverage&COVERAGE=terrapixel&FORMAT=GeoTIFF&BBOX=260489,3667254,260490,3667255&WIDTH=2&HEIGHT=2&CRS=EPSG:26915

As you can see in the gdalinfo of the WCS itself, the pixel resolution
is 0.5m, but this request is for 1.0m, since WCS BBOX= coordinates use
pixel-is-point sematics rather than pixel-is-area, this is requesting a
different resolution from the native resolution and is going outside of
the source-data area, which is advertised in the WCS DescribeCoverage as:

<gml:Envelope srsName="EPSG:26915">
  <gml:pos dimension="2">260489.25 3652963.25</gml:pos>
  <gml:pos dimension="2">313516.75 3667254.75</gml:pos>
</gml:Envelope>

The gdalinfo interprets this (and the equivalent <RectifiedGrid>)
information from the WCS correctly with the pixel-is-point semantic and
converts it properly to GDAL's pixel-is-area semantic and displays:

Upper Left  (  260489.000, 3667255.000) ( 95d34' 1.15"W, 33d 7' 3.30"N)
Lower Left  (  260489.000, 3652963.000) ( 95d33'47.69"W, 32d59'19.66"N)
Upper Right (  313517.000, 3667255.000) ( 94d59'56.38"W, 33d 7'40.78"N)
Lower Right (  313517.000, 3652963.000) ( 94d59'45.90"W, 32d59'56.96"N)
Center      (  287003.000, 3660109.000) ( 95d16'52.90"W, 33d 3'31.34"N)

The same thing happens when requesting live data, which elevates this
oddity to a bug.  If I enter the following into my own application (which
uses pixel-is-area semantics like GDAL does) for a 10m x 10m window at
the top-left corner:

[dev:~] cwplot gdal_wcs_terrapixel.xml format=gdal bbox=260489,3667245,260499,3667255 cs=26915 width=20 height=20

it results in the GDAL generating the WCS request:

base_url?datastore=cwtile&SERVICE=WCS&VERSION=1.0.0&REQUEST=GetCoverage&COVERAGE=terrapixel&FORMAT=GeoTIFF&BBOX=260489,3667245,260499,3667255&WIDTH=20&HEIGHT=20&CRS=EPSG:26915

which is wrong since it hasn't been converted to the pixel-is-point
semantic of the WCS.  It's doing an awkward, shifted resampling of the
source data.

Apparently, other people have run into this problem, too:

http://blog.minst.net/2009/12/09/perceived-wcs-inaccuracy
http://jira.codehaus.org/browse/GEOS-3702

Or are WCSes supposed report data extents as pixel-is-point in
GetCapabilities and DescribeCoverage but accept BBOXes as pixel-is-area
in GetCoverage?  That would be bizarre.  I can't find a place in the WCS
1.0.0 spec where it actually explicitly clears this up.  Does anyone know
what MapServer does about this?

--------------------------+----------------------+--------------------------
Dr. Craig S. Bruce        | Ph 819-771-8303 x205 |             CubeWerx Inc.
Senior Software Developer |   Fax 819-771-8388   |  Gatineau, Québec, Canada
csbruce at cubewerx.com      |  http://csbruce.com/ |  http://www.cubewerx.com/
--------------------------+----------------------+--------------------------
         "And now for tonight's hockey scores -- ZERO!! Ha-ha-ha!"