WCS Server does not return data at the correct resolution

[ben.tuttle@…]

If starting with a 30 arcsec (0.00833333 degrees) data set in a mapserver WCS,
when I make a WCS request using resx=0.00833333 and resy=0.00833333 the output
image returned by mapserver does not register correctly and gdalinfo returns
Pixel Size = (0.00836820,-0.00836820). The new pixel size can vary slightly
based on the bounding box of the request. The main point is that the output
image (we are using geotiff) is not the same resolution as the request. 

Example:
http://www.ngdc.noaa.gov/eog/maps/test5.php?SERVICE=WCS&REQUEST=GetCoverage&Coverage=vis_mosaic&STYLES=&CRS=EPSG:4326&BBOX=-125,30,-115,40&resx=0.00833333&resy=0.00833333&Format=geotiff
the data comes back as  0.00834028,-0.00834028 instead of
0.00833333,-0.00833333 (the native resolution of the data is
0.00833333,-0.00833333). Additionally the given bounding box actually
covers an area 1201px by 1201px in the data. However, the returned
image is 1200px by 1200px. If compared to the original data one can
see that mapserver has resampled the returned image to fit the area
which is 1201 by 1201 into an image that is 1200 by 1200, which
explains the change in the resolution. It seems to me that this may
not be the the desired behaviour for WCS. I would expect to get my data
back at the same resolution requested.

[fwarmerdam]

Added Steve to cc: list for input. 

It would be helpful to have a simple mapfile, and a data file that
setups up a service that demonstrates the problem.  If not, I'll try 
setting something up myself. 

[ben.tuttle@…]

Sorry about that, here is the mapfile from the example in the first post. Also,
the data is still quite large (410M) zipped so I put it on our ftp site
(ftp://ftp.ngdc.noaa.gov/DMSP/ms/test_data.tgz). If anyone would prefer to have
it on the bugzilla server or would prefer a subset on the bugzilla server let me
know. 
::::::::::::::::::::::::::::::::::::::::::::
NAME DMSP-Test
STATUS ON
SIZE 864 312
EXTENT -180 -65 180 65
UNITS dd
SHAPEPATH "/mnt/anchor1/eogmaps/data/"
IMAGECOLOR 255 255 255

WEB
  IMAGEPATH "/mnt/anchor1/eogmaps/www/html/tmp/"
  IMAGEURL "/tmp/"
  METADATA
    "wms_title"     "DMSP Demo Server"  ##required
    "wms_onlineresource" "http://eogmaps.ngdc.noaa.gov/cgi-bin/mapserv.exe?"  
##required
    "wms_srs"       "EPSG:4326"  ##recommended
  END
END

PROJECTION
  "init=epsg:4326"   ##required
END

OUTPUTFORMAT
  NAME GEOTIFF
  DRIVER "GDAL/GTiff"
  MIMETYPE "image/tiff"
  IMAGEMODE "BYTE"
  EXTENSION "tif"
END

LAYER
  NAME "vis_mosaic"
  METADATA
    wcs_description "F152003 11 28 mosaic"
    wcs_name "mosaic"
    wcs_label "vis_mosaic"
    ows_srs "EPSG:4326"
    ows_extent "-180 -65 180 65"
    wcs_resolution "0.00833333 -0.00833333"
    wcs_bandcount "1"
    wcs_formats "GEOTIFF"
    wcs_nativeformat "8-bit GeoTIF"
    wcs_rangeset_name "vis_mosaic"
    wcs_rangeset_label "vis"
    wcs_rangeset_axes "bands"
    wcs_bands_rangeitem "_bands"
    wcs_bands_name "bands"
    wcs_bands_semantic "http://www.ngdc.noaa.gov/dmsp/"
    wcs_bands_refSys "European Petroleum Survey Group (EPSG)"
    wcs_bands_label "F152003"
    wcs_bands_values "VIS"
  END
  STATUS off
  DUMP true
  DATA "F1520061123.nt.d.OIS.vis.mos"
  TYPE raster
  PROJECTION
    "init=epsg:4326"
  END
END

LAYER
  NAME "tir_mosaic"
  METADATA
    wcs_description "F152003 11 28 mosaic"
    wcs_name "mosaic"
    wcs_label "tir_mosaic"
    ows_srs "EPSG:4326"
    ows_extent "-180 -65 180 65"
    wcs_resolution "0.00833333 -0.00833333"
    wcs_bandcount "1"
    wcs_formats "GEOTIFF"
    wcs_nativeformat "8-bit GeoTIF"
    wcs_rangeset_name "tir_mosaic"
    wcs_rangeset_label "tir"
    wcs_rangeset_axes "bands"
    wcs_bands_rangeitem "_bands"
    wcs_bands_name "bands"
    wcs_bands_semantic "http://www.ngdc.noaa.gov/dmsp/"
    wcs_bands_refSys "European Petroleum Survey Group (EPSG)"
    wcs_bands_label "F152003"
    wcs_bands_values "TIR"
  END
  STATUS off
  DUMP true
  DATA "F1520061123.nt.d.OIS.tir.mos"
  TYPE raster
  PROJECTION
    "init=epsg:4326"
  END
END

END #Mapfile

[sdlime]

Starting to look into this. May take all night to download the example image,
but once that's here... ;-)

Steve

[ben.tuttle@…]

So here is some additional information based on work we've been doing. In order
to work around this and get back the data at full res (with no resampling) we
wrote some php that recalculates the width and height based on the bounding box
and the resolution of the data. So we wrote the following two functions:

  function round_bbox($lat0,$lon0,$lat1,$lon1, $resx, $resy)
    {
    $nlat0=round($lat0/$resy)*$resy;
    $nlat1=round($lat1/$resy)*$resy;
    $nlon0=round($lon0/$resx)*$resx;
    $nlon1=round($lon1/$resx)*$resx;
    parse_str("nlat0={$nlat0}&nlon0={$nlon0}&nlat1={$nlat1}&nlon1={$nlon1}
",$all);
    return $all;
    }

  function get_width_height($lat0,$lon0,$lat1,$lon1,$resx, $resy)
    {
    $width = ($lon1-$lon0)/$resx + 1;
    $height = ($lat1-$lat0)/$resy + 1;
    parse_str("nwidth={$width}&nheight={$height}",$all);
    return $all;
    }

We decided that if no widht/height or resx/rey were specified we would calculate
width/height based on the bounding box and return the data at full res for the
requested area. Then we decided to handle resx/resy being specified by
calculating the correct width and height for the specified resolution and return
the data for the requested area at the requested resolution. Our last decision
was that if width and height were specified we would not do any recaculating.
This is all based on a default gridsize of 30 arcseconds in the case of our
data. The block of code we are using for this is as follows:

  $gridsz=1.0/120.0;
  $resx=$request->getvaluebyname('resx');
  $resy=$request->getvaluebyname('resy');
  $width=$request->getvaluebyname('width');
  $height=$request->getvaluebyname('height');
  $service=$request->getvaluebyname('service');
  $format=$request->getvaluebyname('format');
  /* lon0 lat0 lon1 lat1 */
  $bbox=explode(',',$request->getvaluebyname('bbox'));
  if (($resx == "") && ($resy == "") && ($width == "") && ($height == "") &&
(preg_match("/^WCS$/i", $service) == 1)){
    /*Round the bbox and calculate the new width and height*/
    $nbbox=round_bbox($bbox[1],$bbox[0],$bbox[3],$bbox[2], $gridsz, $gridsz);
    $nwidth_nheight=get_width_height($nbbox["nlat0"], $nbbox["nlon0"],
$nbbox["nlat1"], $nbbox["nlon1"], $gridsz, $gridsz);
    $request->setParameter("bbox","{$nbbox["nlon0"]}, {$nbbox["nlat0"]},
{$nbbox["nlon1"]}, {$nbbox["nlat1"]}");
    $request->setParameter("width","{$nwidth_nheight["nwidth"]}");
    $request->setParameter("height","{$nwidth_nheight["nheight"]}");
  }elseif (($resx != "") && ($resy != "") && ($width == "") && ($height == "")
&& (preg_match("/^WCS$/i", $service) == 1)){
    /*Round the bbox and caculate the new width and height based on the
requested resolution*/
    $nbbox=round_bbox($bbox[1],$bbox[0],$bbox[3],$bbox[2], $resx, $resy);
    $nwidth_nheight=get_width_height($nbbox["nlat0"], $nbbox["nlon0"],
$nbbox["nlat1"], $nbbox["nlon1"], $resx, $resy);
    $request->setParameter("bbox","{$nbbox["nlon0"]}, {$nbbox["nlat0"]},
{$nbbox["nlon1"]}, {$nbbox["nlat1"]}");
    $request->setParameter("width","{$nwidth_nheight["nwidth"]}");
    $request->setParameter("height","{$nwidth_nheight["nheight"]}");
  }

Thought it might help to see a work around that seemed to solve the issue we
were having. I'd be happy to discuss this in more detail if anyone thinks it
could be helpful.

[sdlime]

Do I read this as there is not a problem with the WCS code? Rather that the 
client needs to do appropriate computational work (based on service metadata) 
before requesting imagery.

Steve

[ben.tuttle@…]

As it stands the if one made a request for data with a native resolution of 30
arcsec using resx=0.00833333 and resy=0.00833333 the returned image would in
fact not be at that resolution. Mapserver seems to resample the data due to a
minor miscalculation in the width and height (it seems to only be off by one
pixel in either direction and in our tests has been very predictable thus we
were able to get around it in php) that corespond to the bounding box. This
suggests an error in the code to me as it does not behave as I would expect. I
would expect the request to return th edata at the resolution requested. Any
thoughts?

[sdlime]

I played around with this one tonite. Seems to me like the problem is in the
msDrawRasterLow function (or lower) someplace. The requested extent, resolution
and computed height/width remain stable and correct all the way through that call.

They way the WCS code works is that the WCS params are used to configure the
mapfile used to drive the service. This ammounts to setting the EXTENT and SIZE
parameters in this case. Those values seem to be just fine. I can create a test
case outside of the WCS code with Ben's mapfile by changing:

SIZE 1200 1200
EXTENT -125 30 -115 40

and running shp2img (I'll post my version). The output tiff is the right size
but has different EXTENT and resolution values. The msAdjustExtent function
likes the values above and doesn't have to modify them (as expected). This
doesn't look to be a WCS problem .

Ideas?

Steve

[sdlime]

shp2img mapfile test case

[sdlime]

Changing component to GDAL...

Steve

[fwarmerdam]

Steve, 

How are you getting the extent of the GeoTIFF?  Map file extents are 
based on the center of the top left pixel to the center of the 
bottom right pixel.  But gdalinfo reports are the top left corner of the
top left pixel to the bottom right corner of the bottom right pixel. 

Is this accounted for?  I also *think* that the WCS extent has the same
interpretation as gdalinfo, and would need to be adjusted before used as
a mapfile extent.  Is this accounted for?  

[sdlime]

Note that I may be too quick to point to the raster code. I would have expected
msAdjustExtent to compute the resolution gdalinfo reports. That's not the case.
I've got a feeling the raster code is using a different cell model than other
spots resulting in this off-by-one error. See MS_CELLSIZE in map.h. Perusing
mapdrawgdal.c I see some coordinate computations for lr x/y that add .5 pixels.
Which seems to be compensating for a center cell pixel model. The comments in
map.h indicate we switched to a UL pixel model to match WMS.

Steve

[sdlime]

Geez, this gets confusing. As I look at msAdjustExtent it sure looks to me like
an "extent" matches the definitions described in the OGC WMS specs. That is,
those values represent the outside edges of the pixels on the corners. (section
6.5.6 in the WMS 1.1.1 spec). I went back to the 4.6 code and it's been that way
for a good while.
 
Steve

[sdlime]

Ok, more testing. The WMS support certainly assumes a center pixel model and if
I adjust the incoming BBOX for WCS by 1/2 the resolution (supplied or computed)
then I get the right output.

Frank, if a BBOX is not supplied then the code defaults to the extent of the
dataset (using msOWSGetLayerExtent). What pixel model does that extent use? Has
is already been adjusted for MapServer?

Back to the WCS component. I'm still concerned msAdjustExtent is off-by-one.

Steve

[ben.tuttle@…]

Hi all-
FYI we have had some interesting problems with tiled geotiffs that seem likely
to be related to this problem as well. I haven't been able to test it enough to
say absolutely it is the same problem, but it seems likely. Basically when using
a tiled tiff (physical tiles, not internal) the returned image (using WCS) has
systematic shifts at the tile boundaries. Thought there was a chance this
information might help track down the issue if it is indeed related. Happy to
discuss this more it it anyone thinks it might help.

[sdlime]

Ben: I committed a fix to CVS HEAD that addresses the differences in how
MapServer and OWS extents are interpreted. I also fixed another problem that
perhaps could be related to your last comment, again in CVS HEAD. Any chance you
can test things on your end and report back?

Steve

[ben.tuttle@…]

I'm at a conference all week, but I will try and get the updated components up
next week and see what happens. Thanks.

[sdlime]

Ben: Any chance you'll get a chance to try the 5.0 CVS version soon? No hurry, 
just want to make sure I get this resolved.

Steve

[ben.tuttle@…]

Steve:
Trying out the 5.0 cvs version is pretty high on my to do list, but January was
pretty hectic. This issue is pretty important to what we are doing here so I
also want to see it resolved. I will try to get this running as soon as I can
and report back. Thanks for your help on this. 

[ben.tuttle@…]

So I grabbed a copy of mapserver from cvs, here is the mapserv -v info for
reference:
MapServer version 4.99 OUTPUT=GIF OUTPUT=PNG OUTPUT=JPEG OUTPUT=WBMP OUTPUT=SVG
SUPPORTS=PROJ SUPPORTS=FREETYPE SUPPORTS=WMS_SERVER SUPPORTS=WMS_CLIENT
SUPPORTS=WFS_SERVER SUPPORTS=WFS_CLIENT SUPPORTS=WCS_SERVER SUPPORTS=GEOS
INPUT=EPPL7 INPUT=POSTGIS INPUT=OGR INPUT=GDAL INPUT=SHAPEFILE

So far I have only done a few tests, but it seems to be returning the correct
resolution now. I am planning to fire some more requests at it to try and be
sure it works in more than just the few cases I've tried thus far. So far so
good though. I've made a few requests with resx=0.00833333&resy=0.00833333 and
every response has been an image with Pixel Size = (0.00833333,-0.00833333).
I'll try to post any further results as soon as possible.

By the way, it does not seem to fix our issue with physical (not internal) tiles
being mosaicked incorrectly. I'll have to look at this some more to see if maybe
we are doing something bad when we make out tiles. Any chance something similar
could be happening somewhere else in the code that causes a similar issue for tiles?

[sdlime]

We need to look at the tiling problem more closely then. I'm not sure I 
understand the issue (haven't thought about it lately). Thanks for checking on 
this and let me know how other tests go.

Steve

[btuttle]

Hi all- Sorry for the long delay on getting back to this. The newer version does seem to fix the resolution issue as I said before. However, I finally noticed one other problem. If this should be filed as a new bug let me know. I think it may be part of the same issue though. Given a file with the following coordinates: Corner Coordinates: Upper Left (-180.0041667, 65.0041667) (180d 0'15.00"W, 65d 0'15.00"N) Lower Left (-180.0041667, -65.0041147) (180d 0'15.00"W, 65d 0'14.81"S) Upper Right ( 180.0040227, 65.0041667) (180d 0'14.48"E, 65d 0'15.00"N) Lower Right ( 180.0040227, -65.0041147) (180d 0'14.48"E, 65d 0'14.81"S)

When I make a request for the file via WCS I get back the expected resolution and file size, however, the header information is changed to be exactly what I said in the WCS bbox. THis means that the header info is not correct for the downloaded file.

For example, I would expect a request for a bbox=-180,0,-100,65 to return an upper left of -180.0041667, 65.0041667 for this image. Instead it returns -180.0000000, 65.0000000.

I would also expect a request for a bbox=-96.1234567,20,-90,31.1234567 to return an upper left something like -96.0042, -31.0042 rather than -96.1234557,31.1234557.

Unfortunately I don;t have the newer version running on a public server yet, but I would be happy to provide data, mapfiles, set up (if someone has a server I can log into), and/or answers to questions to help work this out.

Let me know if this makes sense or not and if there are questions or anything else I can do to track this down. Also, as I mentioned if this should be filed as a new bug let me know. Thanks.

[warmerdam]

Ben,

This should be addressed seperately from this ticket, so I'll close this one.

I'm not exactly sure that what you describe is a bug in MapServer. My understanding is that mapserver is expected to resample available data into a raster with the size and extents requested regardless of the extent or resolution of the source data.

But we can pursue that elsewhere.