Image from WCS GetCoverage request contain wrong pixel size and wrong corner coordinates

[tbonfort]

Reporter: nacional@cbs.umn.edu Date: 2004/09/14 - 00:49

The bounds and resolution of the resulting tiff file don't match the
WCS GetCoverage request.  The resulting image resolution doesn't correspond to
the request resolution, returning values of almost 2500 x 2500 meters when the
requested resolution is 500 x 500 meters.  The bounding box returned is also
substantially larger than the requested BBOX.

Here's an example query:
http://hypnos.cbs.umn.edu/cgi-bin/mapserv43d?map=/data/wcs/demo.map&SERVICE=WCS&VERSION=1.0.0&REQUEST=GetCoverage&COVERAGE=modis&TIME=2002-169&CRS=EPSG:26915&FORMAT=GEOTIFF_RGB&RESX=500&RESY=500&BBOX=548691.558824,5117528.683007,703816.558824,5241628.683007&bands=2,1,4

Software info:
Gentoo Linux (2.6.8-df8727e3aa04f92baa92fc58473d49e5fc63c041 (r3) gentoo-dev kernel)
GCC 3.3.4
Apache 2.0.50

MapServer CVS:
MapServer version 4.3 OUTPUT=GIF OUTPUT=PNG OUTPUT=JPEG OUTPUT=WBMP OUTPUT=PDF
OUTPUT=SWF SUPPORTS=PROJ SUPPORTS=FREETYPE SUPPORTS=WMS_SERVER
SUPPORTS=WMS_CLIENT SUPPORTS=WFS_SERVER SUPPORTS=WFS_CLIENT SUPPORTS=WCS_SERVER
INPUT=EPPL7 INPUT=POSTGIS INPUT=OGR INPUT=GDAL INPUT=SHAPEFILE

GDAL CVS (formats):
gxf gtiff hfa aigrid aaigrid ceos ceos2 iso8211 xpm sdts raw dted mem jdem
envisat elas fit vrt usgsdem l1b nitf bmp pcidsk dods bsb hdf4 gif jpeg png grass

[tbonfort]

Author: sdlime Date: 2004/09/17 - 21:13

I'll take a look this weekend. Does the output image match the request with 
just the metadata being off or is the image wrong too?

Steve

[tbonfort]

Author: nacional@cbs.umn.edu Date: 2004/09/17 - 22:01

Steve,

Correct, as far as I know the output image has the correct dimensions as
requested.  It's the metadata that's not correct.  If I remove the internal
metadata for the GeoTIFF image and substitute it with an external world file,
the image overlays with other data properly.

-Perry

[tbonfort]

Author: sdlime Date: 2004/09/20 - 21:18

Frank: Where does the data GDAL writes to an image header come from? Is it 
someplace different than where the data comes from to actually build the image? 
Seems like a very odd error...

My guess is that the metadata object we use in the WCS code to harvest metadata 
from a coverage is getting mixed somehow with computed values that end up 
stored in the mapfile... 

Steve

[tbonfort]

Author: sdlime Date: 2004/10/04 - 05:57

This may well be a projection or datum conversion issue.  If I run MapServer with the 
following command:

~/dev/mapserver/mapserv QUERY_STRING="map=/users/sdlime/dev/projects/wcs/data/
modis/2002/demo.map&SERVICE=WCS&VERSION=1.0.0&REQUEST=GetCoverage&
COVERAGE=ndvi&TIME=2002-001&CRS=EPSG:26915&FORMAT=GEOTIFF_INT16&
RESX=500&RESY=500&BBOX=
734290.000000,4945151.000000,1044540.000000,5193351.000000"

The computed extent and cellsize during execution is exactly as requested. The resulting 
tiff image is like so, with screwed up extent and cellsize:

Driver: GTiff/GeoTIFFSize is 621, 496
Coordinate System is:
PROJCS["NAD83 / UTM zone 15N",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010002,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree (supplier to define representation)",0.01745329251994328],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-93.0000000000001],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","26915"]]
Origin = (159707.000000,5501395.000000)
Pixel Size = (2482.00000000,-2260.00000000)
Corner Coordinates:
Upper Left  (  159707.000, 5501395.000) 
Lower Left  (  159707.000, 4380435.000) 
Upper Right ( 1701029.000, 5501395.000) 
Lower Right ( 1701029.000, 4380435.000) 
Center      (  930368.000, 4940915.000) 
Band 1 Block=621x6 Type=Int16, ColorInterp=Gray

The source image gives us this (from gdal):

Driver: ENVI/ENVI .hdr Labelled
Size is 2481, 1807
Coordinate System is:
PROJCS["UTM Zone 15, Northern Hemisphere",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9108"]],
        AXIS["Lat",NORTH],
        AXIS["Long",EAST],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-93],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
Origin = (159707.000000,5501395.000000)
Pixel Size = (500.00000000,-500.00000000)
Corner Coordinates:
Upper Left  (  159707.000, 5501395.000) 
Lower Left  (  159707.000, 4597895.000) 
Upper Right ( 1400207.000, 5501395.000) 
Lower Right ( 1400207.000, 4597895.000) 
Center      (  779957.000, 5049645.000) 
Band 1 Block=2481x1 Type=Int16, ColorInterp=Undefined

The header file for that source image is like so:

ENVI
description = {MOD13A12002001ugl.hdr}
samples = 2481 
lines = 1807 
bands = 1
data type = 2
header offset = 0
file type = ENVI Standard
interleave = bip
sensor type = MODIS
byte order = 0
x start = 0
y start = 0
map info = {UTM, 1.0, 1.0, 159707.000000, 5501395.000000, 500.000000, 500.000000, 
15, North, units=Meters}

The projection code in the mapfile is EPSG:26915. Swtiching the request SRS and the 
layer projections to something that matches the source image (EPSG:32615) doesn't fix 
the problem. The header still doesn't match the request.  Frank does this make and 
sense to you?

Steve

[tbonfort]

Author: sdlime Date: 2004/10/05 - 06:48

Straight shp2img (no projections or anything defined) works ok. Man, very wierd
cause the extent, cellsize, width and height mapObj members are all correct
before handing off to the GDAL drawing code. Frank, any ideas?

Steve

[tbonfort]

Author: sdlime Date: 2004/10/07 - 05:57

The WCS code was missing a call to msMapComputeGeotransform. Apparently that was
initialized to the default size and exent values in the map file hence the funky
GDAL output. msDrawMap does call this function so that's why shp2img worked and
the WCS code (WCS uses msDrawRasterLow) didn't. Should really consider adding a
msPrepareImage (ala MapScript) function and using it rather than duplicating all
that image creation code in the WCS getCoverage function.

Steve