Odd raster resampling problems

[dmorissette]

Frank,

I ran into an odd problem with a WMS layer in an application, and as I was
debugging I kept getting more and more confused.  Perhaps you can show me the
light... the same problem seems to happen with MapServer 4.2 or 4.3.

In mapresample.c's msResampleGDALToMap(), we have the following block of code:


/* -------------------------------------------------------------------- */
/*      Project desired extents out by 2 pixels, and then strip to      */
/*      available data.                                                 */
/* -------------------------------------------------------------------- */
    sSrcExtent.minx = floor(sSrcExtent.minx-1.0);
    sSrcExtent.maxx = ceil (sSrcExtent.maxx+1.0);
    sSrcExtent.miny = floor(sSrcExtent.miny-1.0);
    sSrcExtent.maxy = ceil (sSrcExtent.maxy+1.0);
        
    sSrcExtent.minx = MAX(0,sSrcExtent.minx);
    sSrcExtent.maxx = MIN(sSrcExtent.maxx, nSrcXSize );
    sSrcExtent.miny = MAX(sSrcExtent.miny, 0 );
    sSrcExtent.maxy = MIN(sSrcExtent.maxy, nSrcYSize );
        
    if( sSrcExtent.maxx <= sSrcExtent.minx 
        || sSrcExtent.maxy <= sSrcExtent.miny )
    {
        if( layer->debug )
            msDebug( "msResampleGDALToMap(): no overlap ... no result.\n" );
        return 0;
    }



If I look at sSrcExtent in the debugger before this block of code, I have the
following which seems to make sense:

(gdb) p sSrcExtent
$6 = {minx = 272501.79575150425, miny = 4909536.7036967305, 
  maxx = 281689.75219360756, maxy = 4916526.2247099709}

This corresponds to the BBOX that was used for this layer's WMS GetMap request
in the projection supported by the remote WMS layer.

Now, the calls to MAX() and MIN() above compare ground coordinates with integer
(pixel) values nSrcXSize and nSrcYSize, which results in the following value for
sSrcExtent at the end:

(gdb) p sSrcExtent
$8 = {minx = 272500, miny = 4909535, maxx = 400, maxy = 300}


... and of course the if() fails and I get a message saying "no overlap... no
result".

Am I missing something?  How come this has always worked before and I would be
the first one to run into that?


I'll attach a mapfile to reproduce. All you need to reproduce this is to copy
this mapfile somewhere, and call mapserv in the debugger using:

./mapserv
"QUERY_STRING=map=/path/to/wmsbug.map&version=1.0.0&request=GetMap&layers=usgs+topos"

Using shp2img might work as well... I'm just more used to debugging with the
mapserv CGI.

[dmorissette]

wmsbug.map - WMS client mapfile to reproduce this

[dmorissette]

Here is the WMS GetMap request that is sent to the remote server when I run this
mapfile:
http://terraservice.net/ogcmap.ashx?VERSION=1.1.1&SERVICE=wms&LAYERS=DRG&FORMAT=png&styles=&REQUEST=GetMap&HEIGHT=300&SRS=EPSG:26919&WIDTH=400&BBOX=272501%2E795752%2C4909507%2E524654%2C281689%2E752194%2C4916526%2E224710&EXCEPTIONS=application/vnd.ogc.se_inimage

[dmorissette]

okay, found it: I had forgotten to set IMAGEPATH in my mapfile... I'm not sure
why that matters since the temp images and .wld files were created in the
mapfile's directory by default, but adding the following gives me valid results:

WEB
 IMAGEPATH "/tmp/ms_tmp"
END

... and walking through the code in the debugger I see that sSrcExtents always
has a value in pixels now...

Marking as invalid.

[dmorissette]

More details: I found that the temp images were created in the debugger's CWD,
and the .wld files were created in the mapfile's directory.

So I'll reopen and make msTmpFile() more robust: I'll modify it to use the
msBuildPath() to always return an absolute path, relative to the mapfile's
directory.

[fwarmerdam]

Daniel, 

I can see why if the world file wasn't found, it would screw up the georeferencing
in the manner you saw.  So hopefully making the temporary file creation stuff
bulletproof will resolve the problem. 

The temporary file issue is is the WMS client code, right?

[dmorissette]

Done in 4.3 CVS. msTmpFile() has been modified to take an extra argument (the
mapfile path), and always return absolute paths even if the temporary path is
relative or empty.

[dmorissette]

> The temporary file issue is is the WMS client code, right?

In this case yes, but the same type of issue could have happened with other
blocks of code that used msTmpFile() if the temporary directory (i.e. imagepath)
was relative instead of being absolute (or if it was empty).  In this case the
old msTmpFile() would have still returned a relative path and if that same path
had been run through msBuildPath() later on it could have been modified to point
to a different location.  This is what was happening in the WMS client case
where the temporary images were created in the CWD (relative path), but the .wld
files were created in the mapfile directory because msBuildPath() was used only
for the .wld files.