maprasterquery.c 56 KB
Newer Older
1
/******************************************************************************
2
 * $Id$
3 4
 *
 * Project:  MapServer
5
 * Purpose:  Implementation of query operations on rasters.
6 7 8 9 10 11 12 13 14 15 16 17
 * Author:   Frank Warmerdam, warmerdam@pobox.com
 *
 ******************************************************************************
 * Copyright (c) 1996-2005 Regents of the University of Minnesota.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
18
 * The above copyright notice and this permission notice shall be included in
19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
 * all copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 ****************************************************************************/

#include <assert.h>
#include "mapserver.h"
#include "mapresample.h"
#include "mapthread.h"
34
#include "mapraster.h"
35

36

37
int msRASTERLayerGetShape(layerObj *layer, shapeObj *shape, resultObj *record);
38 39 40 41 42 43 44 45 46 47
int msRASTERLayerGetItems(layerObj *layer);

#ifdef USE_GDAL

/* ==================================================================== */
/*      For now the rasterLayerInfo lives here since it is just used    */
/*      to hold information related to queries.                         */
/* ==================================================================== */
typedef struct {

48 49 50 51 52 53 54
  /* query cache results */
  int query_results;
  int query_alloc_max;
  int query_request_max;
  int query_result_hard_max;
  int raster_query_mode;
  int band_count;
55

56
  int refcount;
57

58 59
  rectObj which_rect;
  int     next_shape;
60

61 62 63 64 65 66 67 68 69 70 71
  double *qc_x;
  double *qc_y;
  double *qc_x_reproj;
  double *qc_y_reproj;
  float *qc_values;
  int    *qc_class;
  int    *qc_red;
  int    *qc_green;
  int    *qc_blue;
  int    *qc_count;
  int    *qc_tileindex;
72

73 74
  /* query bound in force */
  shapeObj *searchshape;
75

76 77 78 79
  /* Only nearest result to this point. */
  int      range_mode; /* MS_QUERY_SINGLE, MS_QUERY_MULTIPLE or -1 (skip test) */
  double   range_dist;
  pointObj target_point;
80

81
  GDALColorTableH hCT;
82

83
  double   shape_tolerance;
84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107

} rasterLayerInfo;

#define RQM_UNKNOWN               0
#define RQM_ENTRY_PER_PIXEL       1
#define RQM_HIST_ON_CLASS         2
#define RQM_HIST_ON_VALUE         3

extern int InvGeoTransform( double *gt_in, double *gt_out );
#define GEO_TRANS(tr,x,y)  ((tr)[0]+(tr)[1]*(x)+(tr)[2]*(y))

/************************************************************************/
/*                             addResult()                              */
/*                                                                      */
/*      this is a copy of the code in mapquery.c.  Should we prepare    */
/*      a small public API for managing results caches?                 */
/************************************************************************/

static int addResult(resultCacheObj *cache, int classindex, int shapeindex, int tileindex)
{
  int i;

  if(cache->numresults == cache->cachesize) { /* just add it to the end */
    if(cache->cachesize == 0)
108
      cache->results = (resultObj *) malloc(sizeof(resultObj)*MS_RESULTCACHEINCREMENT);
109
    else
110
      cache->results = (resultObj *) realloc(cache->results, sizeof(resultObj)*(cache->cachesize+MS_RESULTCACHEINCREMENT));
111 112 113
    if(!cache->results) {
      msSetError(MS_MEMERR, "Realloc() error.", "addResult()");
      return(MS_FAILURE);
114
    }
115 116 117 118 119 120 121 122
    cache->cachesize += MS_RESULTCACHEINCREMENT;
  }

  i = cache->numresults;

  cache->results[i].classindex = classindex;
  cache->results[i].tileindex = tileindex;
  cache->results[i].shapeindex = shapeindex;
123
  cache->results[i].resultindex = -1; /* unused */
124 125 126 127 128 129 130 131 132 133 134 135
  cache->numresults++;

  return(MS_SUCCESS);
}

/************************************************************************/
/*                       msRasterLayerInfoFree()                        */
/************************************************************************/

static void msRasterLayerInfoFree( layerObj *layer )

{
136
  rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;
137

138 139
  if( rlinfo == NULL )
    return;
140

141 142 143 144 145 146
  if( rlinfo->qc_x != NULL ) {
    free( rlinfo->qc_x );
    free( rlinfo->qc_y );
    free( rlinfo->qc_x_reproj );
    free( rlinfo->qc_y_reproj );
  }
147

148 149
  if( rlinfo->qc_values )
    free( rlinfo->qc_values );
150

151 152 153
  if( rlinfo->qc_class ) {
    free( rlinfo->qc_class );
  }
154

155 156 157 158 159
  if( rlinfo->qc_red ) {
    free( rlinfo->qc_red );
    free( rlinfo->qc_green );
    free( rlinfo->qc_blue );
  }
160

161 162 163 164 165
  if( rlinfo->qc_count != NULL )
    free( rlinfo->qc_count );

  if( rlinfo->qc_tileindex != NULL )
    free( rlinfo->qc_tileindex );
166

167
  free( rlinfo );
168

169
  layer->layerinfo = NULL;
170 171 172 173 174 175 176 177 178
}

/************************************************************************/
/*                    msRasterLayerInfoInitialize()                     */
/************************************************************************/

static void msRasterLayerInfoInitialize( layerObj *layer )

{
179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205
  rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;

  if( rlinfo != NULL )
    return;

  rlinfo = (rasterLayerInfo *) msSmallCalloc(1,sizeof(rasterLayerInfo));
  layer->layerinfo = rlinfo;

  rlinfo->band_count = -1;
  rlinfo->raster_query_mode = RQM_ENTRY_PER_PIXEL;
  rlinfo->range_mode = -1; /* inactive */
  rlinfo->refcount = 0;
  rlinfo->shape_tolerance = 0.0;

  /* We need to do this or the layer->layerinfo will be interpreted */
  /* as shapefile access info because the default connectiontype is */
  /* MS_SHAPEFILE. */
  if (layer->connectiontype != MS_WMS)
    layer->connectiontype = MS_RASTER;

  rlinfo->query_result_hard_max = 1000000;

  if( CSLFetchNameValue( layer->processing, "RASTER_QUERY_MAX_RESULT" )
      != NULL ) {
    rlinfo->query_result_hard_max =
      atoi(CSLFetchNameValue( layer->processing, "RASTER_QUERY_MAX_RESULT" ));
  }
206 207 208 209 210 211
}

/************************************************************************/
/*                       msRasterQueryAddPixel()                        */
/************************************************************************/

212 213
static void msRasterQueryAddPixel( layerObj *layer, pointObj *location,
                                   pointObj *reprojectedLocation,
214 215 216
                                   float *values )

{
217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262
  rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;
  int red = 0, green = 0, blue = 0, nodata = FALSE;
  int p_class = -1;

  if( rlinfo->query_results == rlinfo->query_result_hard_max )
    return;

  /* -------------------------------------------------------------------- */
  /*      Is this our first time in?  If so, do an initial allocation     */
  /*      for the data arrays suitable to our purposes.                   */
  /* -------------------------------------------------------------------- */
  if( rlinfo->query_alloc_max == 0 ) {
    rlinfo->query_alloc_max = 2;

    switch( rlinfo->raster_query_mode ) {
      case RQM_ENTRY_PER_PIXEL:
        rlinfo->qc_x = (double *)
                       msSmallCalloc(sizeof(double),rlinfo->query_alloc_max);
        rlinfo->qc_y = (double *)
                       msSmallCalloc(sizeof(double),rlinfo->query_alloc_max);
        rlinfo->qc_x_reproj = (double *)
                       msSmallCalloc(sizeof(double),rlinfo->query_alloc_max);
        rlinfo->qc_y_reproj = (double *)
                       msSmallCalloc(sizeof(double),rlinfo->query_alloc_max);
        rlinfo->qc_values = (float *)
                            msSmallCalloc(sizeof(float),
                                          rlinfo->query_alloc_max*rlinfo->band_count);
        rlinfo->qc_red = (int *)
                         msSmallCalloc(sizeof(int),rlinfo->query_alloc_max);
        rlinfo->qc_green = (int *)
                           msSmallCalloc(sizeof(int),rlinfo->query_alloc_max);
        rlinfo->qc_blue = (int *)
                          msSmallCalloc(sizeof(int),rlinfo->query_alloc_max);
        if( layer->numclasses > 0 )
          rlinfo->qc_class = (int *)
                             msSmallCalloc(sizeof(int),rlinfo->query_alloc_max);
        break;

      case RQM_HIST_ON_CLASS:
        break;

      case RQM_HIST_ON_VALUE:
        break;

      default:
        assert( FALSE );
263
    }
264
  }
265

266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300
  /* -------------------------------------------------------------------- */
  /*      Reallocate the data arrays larger if they are near the max      */
  /*      now.                                                            */
  /* -------------------------------------------------------------------- */
  if( rlinfo->query_results == rlinfo->query_alloc_max ) {
    rlinfo->query_alloc_max = rlinfo->query_alloc_max * 2 + 100;

    if( rlinfo->qc_x != NULL )
      rlinfo->qc_x = msSmallRealloc(rlinfo->qc_x,
                                    sizeof(double) * rlinfo->query_alloc_max);
    if( rlinfo->qc_y != NULL )
      rlinfo->qc_y = msSmallRealloc(rlinfo->qc_y,
                                    sizeof(double) * rlinfo->query_alloc_max);
    if( rlinfo->qc_x_reproj != NULL )
      rlinfo->qc_x_reproj = msSmallRealloc(rlinfo->qc_x_reproj,
                                    sizeof(double) * rlinfo->query_alloc_max);
    if( rlinfo->qc_y_reproj != NULL )
      rlinfo->qc_y_reproj = msSmallRealloc(rlinfo->qc_y_reproj,
                                    sizeof(double) * rlinfo->query_alloc_max);
    if( rlinfo->qc_values != NULL )
      rlinfo->qc_values =
        msSmallRealloc(rlinfo->qc_values,
                       sizeof(float) * rlinfo->query_alloc_max
                       * rlinfo->band_count );
    if( rlinfo->qc_class != NULL )
      rlinfo->qc_class = msSmallRealloc(rlinfo->qc_class,
                                        sizeof(int) * rlinfo->query_alloc_max);
    if( rlinfo->qc_red != NULL )
      rlinfo->qc_red = msSmallRealloc(rlinfo->qc_red,
                                      sizeof(int) * rlinfo->query_alloc_max);
    if( rlinfo->qc_green != NULL )
      rlinfo->qc_green = msSmallRealloc(rlinfo->qc_green,
                                        sizeof(int) * rlinfo->query_alloc_max);
    if( rlinfo->qc_blue != NULL )
      rlinfo->qc_blue = msSmallRealloc(rlinfo->qc_blue,
301
                                       sizeof(int) * rlinfo->query_alloc_max);
302 303 304 305 306 307 308
    if( rlinfo->qc_count != NULL )
      rlinfo->qc_count = msSmallRealloc(rlinfo->qc_count,
                                        sizeof(int) * rlinfo->query_alloc_max);
    if( rlinfo->qc_tileindex != NULL )
      rlinfo->qc_tileindex = msSmallRealloc(rlinfo->qc_tileindex,
                                            sizeof(int) * rlinfo->query_alloc_max);
  }
309

310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326
  /* -------------------------------------------------------------------- */
  /*      Handle colormap                                                 */
  /* -------------------------------------------------------------------- */
  if( rlinfo->hCT != NULL ) {
    int pct_index = (int) floor(values[0]);
    GDALColorEntry sEntry;

    if( GDALGetColorEntryAsRGB( rlinfo->hCT, pct_index, &sEntry ) ) {
      red = sEntry.c1;
      green = sEntry.c2;
      blue = sEntry.c3;

      if( sEntry.c4 == 0 )
        nodata = TRUE;
    } else
      nodata = TRUE;
  }
327

328 329 330 331 332 333 334 335 336 337
  /* -------------------------------------------------------------------- */
  /*      Color derived from greyscale value.                             */
  /* -------------------------------------------------------------------- */
  else {
    if( rlinfo->band_count >= 3 ) {
      red = (int) MAX(0,MIN(255,values[0]));
      green = (int) MAX(0,MIN(255,values[1]));
      blue = (int) MAX(0,MIN(255,values[2]));
    } else {
      red = green = blue = (int) MAX(0,MIN(255,values[0]));
338
    }
339
  }
340

341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364
  /* -------------------------------------------------------------------- */
  /*      Handle classification.                                          */
  /*                                                                      */
  /*      NOTE: The following is really quite inadequate to deal with     */
  /*      classifications based on [red], [green] and [blue] as           */
  /*      described in:                                                   */
  /*       http://mapserver.gis.umn.edu/bugs/show_bug.cgi?id=1021         */
  /* -------------------------------------------------------------------- */
  if( rlinfo->qc_class != NULL ) {
    p_class = msGetClass_FloatRGB(layer, values[0],
                                  red, green, blue );

    if( p_class == -1 )
      nodata = TRUE;
    else {
      nodata = FALSE;
      rlinfo->qc_class[rlinfo->query_results] = p_class;
      if( layer->class[p_class]->numstyles > 0 ) {
        red   = layer->class[p_class]->styles[0]->color.red;
        green = layer->class[p_class]->styles[0]->color.green;
        blue  = layer->class[p_class]->styles[0]->color.blue;
      } else {
        red = green = blue = 0;
      }
365
    }
366
  }
367

368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383
  /* -------------------------------------------------------------------- */
  /*      Record the color.                                               */
  /* -------------------------------------------------------------------- */
  rlinfo->qc_red[rlinfo->query_results] = red;
  rlinfo->qc_green[rlinfo->query_results] = green;
  rlinfo->qc_blue[rlinfo->query_results] = blue;

  /* -------------------------------------------------------------------- */
  /*      Record spatial location.                                        */
  /* -------------------------------------------------------------------- */
  if( rlinfo->qc_x != NULL ) {
    rlinfo->qc_x[rlinfo->query_results] = location->x;
    rlinfo->qc_y[rlinfo->query_results] = location->y;
    rlinfo->qc_x_reproj[rlinfo->query_results] = reprojectedLocation->x;
    rlinfo->qc_y_reproj[rlinfo->query_results] = reprojectedLocation->y;
  }
384

385 386 387 388 389 390 391 392 393 394 395 396 397 398
  /* -------------------------------------------------------------------- */
  /*      Record actual pixel value(s).                                   */
  /* -------------------------------------------------------------------- */
  if( rlinfo->qc_values != NULL )
    memcpy( rlinfo->qc_values + rlinfo->query_results * rlinfo->band_count,
            values, sizeof(float) * rlinfo->band_count );

  /* -------------------------------------------------------------------- */
  /*      Add to the results cache.                                       */
  /* -------------------------------------------------------------------- */
  if( ! nodata ) {
    addResult( layer->resultcache, p_class, rlinfo->query_results, 0 );
    rlinfo->query_results++;
  }
399 400 401 402 403 404
}

/************************************************************************/
/*                       msRasterQueryByRectLow()                       */
/************************************************************************/

405
static int
406
msRasterQueryByRectLow(mapObj *map, layerObj *layer, GDALDatasetH hDS,
407
                       rectObj queryRect)
408 409

{
410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427
  double  adfGeoTransform[6], adfInvGeoTransform[6];
  double      dfXMin, dfYMin, dfXMax, dfYMax, dfX, dfY, dfAdjustedRange;
  int         nWinXOff, nWinYOff, nWinXSize, nWinYSize;
  int         nRXSize, nRYSize;
  float       *pafRaster;
  int         nBandCount, *panBandMap, iPixel, iLine;
  CPLErr      eErr;
  rasterLayerInfo *rlinfo;
  rectObj     searchrect;
  int         needReproject = MS_FALSE;

  rlinfo = (rasterLayerInfo *) layer->layerinfo;

  /* -------------------------------------------------------------------- */
  /*      Reproject the search rect into the projection of this           */
  /*      layer/file if needed.                                           */
  /* -------------------------------------------------------------------- */
  searchrect = queryRect;
428
#ifdef USE_PROJ
429 430 431 432 433 434
  if(layer->project
      && msProjectionsDiffer(&(layer->projection), &(map->projection))) {
    msProjectRect(&(map->projection), &(layer->projection), &searchrect);
    needReproject = MS_TRUE;
  } else
    layer->project = MS_FALSE;
435 436
#endif

437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531
  /* -------------------------------------------------------------------- */
  /*      Transform the rectangle in target ground coordinates to         */
  /*      pixel/line extents on the file.  Process all 4 corners, to      */
  /*      build extents.                                                  */
  /* -------------------------------------------------------------------- */
  nRXSize = GDALGetRasterXSize( hDS );
  nRYSize = GDALGetRasterYSize( hDS );

  msGetGDALGeoTransform( hDS, map, layer, adfGeoTransform );
  InvGeoTransform( adfGeoTransform, adfInvGeoTransform );

  /* top left */
  dfXMin = dfXMax = GEO_TRANS(adfInvGeoTransform,
                              searchrect.minx, searchrect.maxy);
  dfYMin = dfYMax = GEO_TRANS(adfInvGeoTransform+3,
                              searchrect.minx, searchrect.maxy);

  /* top right */
  dfX = GEO_TRANS(adfInvGeoTransform  , searchrect.maxx, searchrect.maxy);
  dfY = GEO_TRANS(adfInvGeoTransform+3, searchrect.maxx, searchrect.maxy);
  dfXMin = MIN(dfXMin,dfX);
  dfXMax = MAX(dfXMax,dfX);
  dfYMin = MIN(dfYMin,dfY);
  dfYMax = MAX(dfYMax,dfY);

  /* bottom left */
  dfX = GEO_TRANS(adfInvGeoTransform  , searchrect.minx, searchrect.miny);
  dfY = GEO_TRANS(adfInvGeoTransform+3, searchrect.minx, searchrect.miny);
  dfXMin = MIN(dfXMin,dfX);
  dfXMax = MAX(dfXMax,dfX);
  dfYMin = MIN(dfYMin,dfY);
  dfYMax = MAX(dfYMax,dfY);

  /* bottom right */
  dfX = GEO_TRANS(adfInvGeoTransform  , searchrect.maxx, searchrect.miny);
  dfY = GEO_TRANS(adfInvGeoTransform+3, searchrect.maxx, searchrect.miny);
  dfXMin = MIN(dfXMin,dfX);
  dfXMax = MAX(dfXMax,dfX);
  dfYMin = MIN(dfYMin,dfY);
  dfYMax = MAX(dfYMax,dfY);

  /* -------------------------------------------------------------------- */
  /*      Trim the rectangle to the area of the file itself, but out      */
  /*      to the edges of the touched edge pixels.                        */
  /* -------------------------------------------------------------------- */
  dfXMin = MAX(0.0,MIN(nRXSize,floor(dfXMin)));
  dfYMin = MAX(0.0,MIN(nRYSize,floor(dfYMin)));
  dfXMax = MAX(0.0,MIN(nRXSize,ceil(dfXMax)));
  dfYMax = MAX(0.0,MIN(nRYSize,ceil(dfYMax)));

  /* -------------------------------------------------------------------- */
  /*      Convert to integer offset/size values.                          */
  /* -------------------------------------------------------------------- */
  nWinXOff = (int) dfXMin;
  nWinYOff = (int) dfYMin;
  nWinXSize = (int) (dfXMax - dfXMin);
  nWinYSize = (int) (dfYMax - dfYMin);

  /* -------------------------------------------------------------------- */
  /*      What bands are we operating on?                                 */
  /* -------------------------------------------------------------------- */
  panBandMap = msGetGDALBandList( layer, hDS, 0, &nBandCount );

  if( rlinfo->band_count == -1 )
    rlinfo->band_count = nBandCount;

  if( nBandCount != rlinfo->band_count ) {
    msSetError( MS_IMGERR,
                "Got %d bands, but expected %d bands.",
                "msRasterQueryByRectLow()",
                nBandCount, rlinfo->band_count );

    return -1;
  }

  /* -------------------------------------------------------------------- */
  /*      Try to load the raster data.  For now we just load the first    */
  /*      band in the file.  Later we will deal with the various band     */
  /*      selection criteria.                                             */
  /* -------------------------------------------------------------------- */
  pafRaster = (float *)
              calloc(sizeof(float),nWinXSize*nWinYSize*nBandCount);
  MS_CHECK_ALLOC(pafRaster, sizeof(float)*nWinXSize*nWinYSize*nBandCount, -1);

  eErr = GDALDatasetRasterIO( hDS, GF_Read,
                              nWinXOff, nWinYOff, nWinXSize, nWinYSize,
                              pafRaster, nWinXSize, nWinYSize, GDT_Float32,
                              nBandCount, panBandMap,
                              4 * nBandCount,
                              4 * nBandCount * nWinXSize,
                              4 );

  if( eErr != CE_None ) {
    msSetError( MS_IOERR, "GDALDatasetRasterIO() failed: %s",
                "msRasterQueryByRectLow()", CPLGetLastErrorMsg() );
532 533

    free( pafRaster );
534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611
    return -1;
  }

  /* -------------------------------------------------------------------- */
  /*      Fetch color table for intepreting colors if needed.             */
  /* -------------------------------------------------------------------- */
  rlinfo->hCT = GDALGetRasterColorTable(
                  GDALGetRasterBand( hDS, panBandMap[0] ) );

  free( panBandMap );

  /* -------------------------------------------------------------------- */
  /*      When computing whether pixels are within range we do it         */
  /*      based on the center of the pixel to the target point but        */
  /*      really it ought to be the nearest point on the pixel.  It       */
  /*      would be too much trouble to do this rigerously, so we just     */
  /*      add a fudge factor so that a range of zero will find the        */
  /*      pixel the target falls in at least.                             */
  /* -------------------------------------------------------------------- */
  dfAdjustedRange =
    sqrt(adfGeoTransform[1] * adfGeoTransform[1]
         + adfGeoTransform[2] * adfGeoTransform[2]
         + adfGeoTransform[4] * adfGeoTransform[4]
         + adfGeoTransform[5] * adfGeoTransform[5]) * 0.5 * 1.41421356237
    + sqrt( rlinfo->range_dist );
  dfAdjustedRange = dfAdjustedRange * dfAdjustedRange;

  /* -------------------------------------------------------------------- */
  /*      Loop over all pixels determining which are "in".                */
  /* -------------------------------------------------------------------- */
  for( iLine = 0; iLine < nWinYSize; iLine++ ) {
    for( iPixel = 0; iPixel < nWinXSize; iPixel++ ) {
      pointObj  sPixelLocation,sReprojectedPixelLocation;

      if( rlinfo->query_results == rlinfo->query_result_hard_max )
        break;

      /* transform pixel/line to georeferenced */
      sPixelLocation.x =
        GEO_TRANS(adfGeoTransform,
                  iPixel + nWinXOff + 0.5, iLine + nWinYOff + 0.5 );
      sPixelLocation.y =
        GEO_TRANS(adfGeoTransform+3,
                  iPixel + nWinXOff + 0.5, iLine + nWinYOff + 0.5 );

      /* If projections differ, convert this back into the map  */
      /* projection for distance testing, and comprison to the  */
      /* search shape.  Save the original pixel location coordinates */
      /* in sPixelLocationInLayerSRS, so that we can return those */
      /* coordinates if we have a hit */
      sReprojectedPixelLocation = sPixelLocation;
      if( needReproject )
        msProjectPoint( &(layer->projection), &(map->projection),
                        &sReprojectedPixelLocation);

      /* If we are doing QueryByShape, check against the shape now */
      if( rlinfo->searchshape != NULL ) {
        if( rlinfo->shape_tolerance == 0.0
            && rlinfo->searchshape->type == MS_SHAPE_POLYGON ) {
          if( msIntersectPointPolygon(
                &sReprojectedPixelLocation, rlinfo->searchshape ) == MS_FALSE )
            continue;
        } else {
          shapeObj  tempShape;
          lineObj   tempLine;

          memset( &tempShape, 0, sizeof(shapeObj) );
          tempShape.type = MS_SHAPE_POINT;
          tempShape.numlines = 1;
          tempShape.line = &tempLine;
          tempLine.numpoints = 1;
          tempLine.point = &sReprojectedPixelLocation;

          if( msDistanceShapeToShape(rlinfo->searchshape, &tempShape)
              > rlinfo->shape_tolerance )
            continue;
        }
      }
612

613 614
      if( rlinfo->range_mode >= 0 ) {
        double dist;
615

616 617 618 619
        dist = (rlinfo->target_point.x - sReprojectedPixelLocation.x)
               * (rlinfo->target_point.x - sReprojectedPixelLocation.x)
               + (rlinfo->target_point.y - sReprojectedPixelLocation.y)
               * (rlinfo->target_point.y - sReprojectedPixelLocation.y);
620

621 622 623 624 625 626 627 628
        if( dist >= dfAdjustedRange )
          continue;

        /* If we can only have one feature, trim range and clear */
        /* previous result.  */
        if( rlinfo->range_mode == MS_QUERY_SINGLE ) {
          rlinfo->range_dist = dist;
          rlinfo->query_results = 0;
629
        }
630 631 632 633 634 635 636
      }

      msRasterQueryAddPixel( layer,
			                       &sPixelLocation, // return coords in layer SRS
                             &sReprojectedPixelLocation,
                             pafRaster
                             + (iLine*nWinXSize + iPixel) * nBandCount );
637
    }
638
  }
639

640 641 642 643
  /* -------------------------------------------------------------------- */
  /*      Cleanup.                                                        */
  /* -------------------------------------------------------------------- */
  free( pafRaster );
644

645
  return MS_SUCCESS;
646 647 648 649 650 651 652
}
#endif /* def USE_GDAL */

/************************************************************************/
/*                        msRasterQueryByRect()                         */
/************************************************************************/

653
int msRasterQueryByRect(mapObj *map, layerObj *layer, rectObj queryRect)
654 655 656

{
#ifndef USE_GDAL
657 658 659 660
  msSetError( MS_IMGERR,
              "Rasters queries only supported with GDAL support enabled.",
              "msRasterQueryByRect()" );
  return MS_FAILURE;
661
#else
662 663 664 665
  int status = MS_SUCCESS;
  char *filename=NULL;

  layerObj *tlp=NULL; /* pointer to the tile layer either real or temporary */
666
  int tileitemindex=-1, tilelayerindex=-1, tilesrsindex=-1;
667
  shapeObj tshp;
668
  char tilename[MS_PATH_LENGTH], tilesrsname[1024];
669
  int  done, destroy_on_failure;
670

671 672 673 674 675 676 677
  char szPath[MS_MAXPATHLEN];
  rectObj searchrect;
  rasterLayerInfo *rlinfo = NULL;

  /* -------------------------------------------------------------------- */
  /*      Get the layer info.                                             */
  /* -------------------------------------------------------------------- */
678 679 680 681 682 683
  if(!layer->layerinfo) {
    msRasterLayerInfoInitialize( layer );
    destroy_on_failure = 1;
  } else {
    destroy_on_failure = 0;
  }
684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703
  rlinfo = (rasterLayerInfo *) layer->layerinfo;

  /* -------------------------------------------------------------------- */
  /*      Clear old results cache.                                        */
  /* -------------------------------------------------------------------- */
  if(layer->resultcache) {
    if(layer->resultcache->results) free(layer->resultcache->results);
    free(layer->resultcache);
    layer->resultcache = NULL;
  }

  /* -------------------------------------------------------------------- */
  /*      Initialize the results cache.                                   */
  /* -------------------------------------------------------------------- */
  layer->resultcache = (resultCacheObj *)msSmallMalloc(sizeof(resultCacheObj));
  layer->resultcache->results = NULL;
  layer->resultcache->numresults = layer->resultcache->cachesize = 0;
  layer->resultcache->bounds.minx =
    layer->resultcache->bounds.miny =
      layer->resultcache->bounds.maxx =
704 705
        layer->resultcache->bounds.maxy = -1;

706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729
  /* -------------------------------------------------------------------- */
  /*      Check if we should really be acting on this layer and           */
  /*      provide debug info in various cases.                            */
  /* -------------------------------------------------------------------- */
  if(layer->debug > 0 || map->debug > 1)
    msDebug( "msRasterQueryByRect(%s): entering.\n", layer->name );

  if(!layer->data && !layer->tileindex) {
    if(layer->debug > 0 || map->debug > 0 )
      msDebug( "msRasterQueryByRect(%s): layer data and tileindex NULL ... doing nothing.", layer->name );
    return MS_SUCCESS;
  }

  if((layer->status != MS_ON) && layer->status != MS_DEFAULT ) {
    if(layer->debug > 0 )
      msDebug( "msRasterQueryByRect(%s): not status ON or DEFAULT, doing nothing.", layer->name );
    return MS_SUCCESS;
  }

  /* ==================================================================== */
  /*      Handle setting up tileindex layer.                              */
  /* ==================================================================== */
  if(layer->tileindex) { /* we have an index file */
    msInitShape(&tshp);
730
    searchrect = queryRect;
731

732 733 734 735 736 737 738 739
    status = msDrawRasterSetupTileLayer(map, layer,
                           &searchrect,
                           MS_TRUE,
                           &tilelayerindex,
                           &tileitemindex,
                           &tilesrsindex,
                           &tlp);
    if (status != MS_SUCCESS) {
740 741
      goto cleanup;
    }
742 743 744
  } else {
    /* we have to manually apply to scaletoken logic as we're not going through msLayerOpen() */
    status = msLayerApplyScaletokens(layer,map->scaledenom);
745 746
    if (status != MS_SUCCESS) {
      goto cleanup;
747
    }
748
  }
749

750 751 752 753 754 755 756 757 758 759 760 761 762
  /* -------------------------------------------------------------------- */
  /*      Iterate over all tiles (just one in untiled case).              */
  /* -------------------------------------------------------------------- */
  done = MS_FALSE;
  while( done == MS_FALSE && status == MS_SUCCESS ) {

    GDALDatasetH  hDS;
    char *decrypted_path = NULL;

    /* -------------------------------------------------------------------- */
    /*      Get filename.                                                   */
    /* -------------------------------------------------------------------- */
    if(layer->tileindex) {
763 764 765 766 767
      status = msDrawRasterIterateTileIndex(layer, tlp, &tshp,
                                            tileitemindex, tilesrsindex,
                                            tilename, sizeof(tilename),
                                            tilesrsname, sizeof(tilesrsname));
      if( status == MS_FAILURE) {
768
        break;
769
      }
770 771 772 773 774 775 776

      if(status == MS_DONE) break; /* no more tiles/images */
      filename = tilename;
    } else {
      filename = layer->data;
      done = MS_TRUE; /* only one image so we're done after this */
    }
777

778
    if(strlen(filename) == 0) continue;
779

780 781 782 783
    /* -------------------------------------------------------------------- */
    /*      Open the file.                                                  */
    /* -------------------------------------------------------------------- */
    msGDALInitialize();
784

785
    msDrawRasterBuildRasterPath(map, layer, filename, szPath);
786

787
    decrypted_path = msDecryptStringTokens( map, szPath );
788 789 790 791
    if( !decrypted_path ) {
      status = MS_FAILURE;
      goto cleanup;
    }
792 793 794 795 796 797

    msAcquireLock( TLOCK_GDAL );
    hDS = GDALOpen(decrypted_path, GA_ReadOnly );

    if( hDS == NULL ) {
      int ignore_missing = msMapIgnoreMissingData( map );
798
      const char *cpl_error_msg = msDrawRasterGetCPLErrorMsg(decrypted_path, szPath);
799 800 801 802 803 804 805 806 807 808 809 810

      msFree( decrypted_path );
      decrypted_path = NULL;

      msReleaseLock( TLOCK_GDAL );

      if ( ignore_missing == MS_MISSING_DATA_FAIL ) {
        if( layer->debug || map->debug )
          msSetError( MS_IMGERR,
                      "Unable to open file %s for layer %s ... fatal error.\n%s",
                      "msRasterQueryByRect()",
                      szPath, layer->name, cpl_error_msg);
811 812
        status = MS_FAILURE;
        goto cleanup;
813 814 815 816 817 818 819 820 821 822 823 824
      }
      if( ignore_missing == MS_MISSING_DATA_LOG ) {
        if( layer->debug || map->debug )
          msDebug( "Unable to open file %s for layer %s ... ignoring this missing data.\n%s",
                   filename, layer->name, cpl_error_msg );
      }
      continue;
    }

    msFree( decrypted_path );
    decrypted_path = NULL;

825 826 827 828 829
    if( msDrawRasterLoadProjection(layer, hDS, filename, tilesrsindex, tilesrsname) != MS_SUCCESS )
    {
        msReleaseLock( TLOCK_GDAL );
        status = MS_FAILURE;
        goto cleanup;
830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847
    }

    /* -------------------------------------------------------------------- */
    /*      Perform actual query on this file.                              */
    /* -------------------------------------------------------------------- */
    if( status == MS_SUCCESS )
      status = msRasterQueryByRectLow( map, layer, hDS, queryRect );

    GDALClose( hDS );
    msReleaseLock( TLOCK_GDAL );

  } /* next tile */

  /* -------------------------------------------------------------------- */
  /*      Cleanup tileindex if it is open.                                */
  /* -------------------------------------------------------------------- */
cleanup:
  if(layer->tileindex) { /* tiling clean-up */
848 849 850
    msDrawRasterCleanupTileLayer(tlp, tilelayerindex);
  } else {
    msLayerRestoreFromScaletokens(layer);
851
  }
852

853 854 855 856
  /* -------------------------------------------------------------------- */
  /*      On failure, or empty result set, cleanup the rlinfo since we    */
  /*      likely won't ever have it accessed or cleaned up later.         */
  /* -------------------------------------------------------------------- */
857 858 859 860 861 862
  if( status == MS_FAILURE || rlinfo->query_results == 0 ) {
    if(destroy_on_failure) {
      msRasterLayerInfoFree( layer );
    }
  } else {
    /* populate the items/numitems layer-level values */
863
    msRASTERLayerGetItems(layer);
864
  }
865

866
  return status;
867 868 869 870 871 872 873 874 875 876 877
#endif /* def USE_GDAL */
}

/************************************************************************/
/*                        msRasterQueryByShape()                        */
/************************************************************************/

int msRasterQueryByShape(mapObj *map, layerObj *layer, shapeObj *selectshape)

{
#ifndef USE_GDAL
878 879 880 881
  msSetError( MS_IMGERR,
              "Rasters queries only supported with GDAL support enabled.",
              "msRasterQueryByRect()" );
  return MS_FAILURE;
882
#else
883 884 885 886
  rasterLayerInfo *rlinfo = NULL;
  int status;
  double tolerance;
  rectObj searchrect;
887

888
  msRasterLayerInfoInitialize( layer );
889

890
  rlinfo = (rasterLayerInfo *) layer->layerinfo;
891

892 893 894 895 896 897
  /* If the selection shape is point or line we use the default tolerance of
     3, but for polygons we require an exact hit. */
  if(layer->tolerance == -1) {
    if(selectshape->type == MS_SHAPE_POINT ||
        selectshape->type == MS_SHAPE_LINE)
      tolerance = 3;
898
    else
899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928
      tolerance = 0;
  } else
    tolerance = layer->tolerance;

  if(layer->toleranceunits == MS_PIXELS)
    tolerance = tolerance
                * msAdjustExtent(&(map->extent), map->width, map->height);
  else
    tolerance = tolerance
                * (msInchesPerUnit(layer->toleranceunits,0)
                   / msInchesPerUnit(map->units,0));

  rlinfo->searchshape = selectshape;
  rlinfo->shape_tolerance = tolerance;

  msComputeBounds( selectshape );
  searchrect = selectshape->bounds;

  searchrect.minx -= tolerance; /* expand the search box to account for layer tolerances (e.g. buffered searches) */
  searchrect.maxx += tolerance;
  searchrect.miny -= tolerance;
  searchrect.maxy += tolerance;

  status = msRasterQueryByRect( map, layer, searchrect );

  rlinfo = (rasterLayerInfo *) layer->layerinfo;
  if( rlinfo )
    rlinfo->searchshape = NULL;

  return status;
929 930 931 932 933 934 935
#endif /* USE_GDAL  */
}

/************************************************************************/
/*                        msRasterQueryByPoint()                        */
/************************************************************************/

936
int msRasterQueryByPoint(mapObj *map, layerObj *layer, int mode, pointObj p,
937 938 939
                         double buffer, int maxresults)
{
#ifndef USE_GDAL
940 941 942 943
  msSetError( MS_IMGERR,
              "Rasters queries only supported with GDAL support enabled.",
              "msRasterQueryByRect()" );
  return MS_FAILURE;
944
#else
945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967
  int result;
  int previous_maxresults;
  double layer_tolerance;
  rectObj bufferRect;
  rasterLayerInfo *rlinfo = NULL;

  msRasterLayerInfoInitialize( layer );

  rlinfo = (rasterLayerInfo *) layer->layerinfo;

  /* -------------------------------------------------------------------- */
  /*      If the buffer is not set, then use layer tolerances             */
  /*      instead.   The "buffer" distince is now in georeferenced        */
  /*      units.  Note that tolerances in pixels are basically map        */
  /*      display pixels, not underlying raster pixels.  It isn't         */
  /*      clear that there is any way of requesting a buffer size in      */
  /*      underlying pixels.                                              */
  /* -------------------------------------------------------------------- */
  if(buffer <= 0) { /* use layer tolerance */
    if(layer->tolerance == -1)
      layer_tolerance = 3;
    else
      layer_tolerance = layer->tolerance;
968

969 970 971 972 973 974 975 976
    if(layer->toleranceunits == MS_PIXELS)
      buffer = layer_tolerance
               * msAdjustExtent(&(map->extent), map->width, map->height);
    else
      buffer = layer_tolerance
               * (msInchesPerUnit(layer->toleranceunits,0)
                  / msInchesPerUnit(map->units,0));
  }
977

978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004
  /* -------------------------------------------------------------------- */
  /*      Setup target point information, at this point they are in       */
  /*      map coordinates.                                                */
  /* -------------------------------------------------------------------- */
  rlinfo->range_dist = buffer * buffer;
  rlinfo->target_point = p;

  /* -------------------------------------------------------------------- */
  /*      if we are in the MS_QUERY_SINGLE mode, first try a query with   */
  /*      zero tolerance.  If this gets a raster pixel then we can be     */
  /*      reasonably assured that it is the closest to the query          */
  /*      point.  This will potentially be must more efficient than       */
  /*      processing all pixels within the tolerance.                     */
  /* -------------------------------------------------------------------- */
  if( mode == MS_QUERY_SINGLE ) {
    rectObj pointRect;

    pointRect.minx = p.x;
    pointRect.maxx = p.x;
    pointRect.miny = p.y;
    pointRect.maxy = p.y;

    rlinfo->range_mode = MS_QUERY_SINGLE;
    result = msRasterQueryByRect( map, layer, pointRect );
    if( rlinfo->query_results > 0 )
      return result;
  }
1005

1006 1007 1008 1009 1010 1011 1012 1013
  /* -------------------------------------------------------------------- */
  /*      Setup a rectangle that is everything within the designated      */
  /*      range and do a search against that.                             */
  /* -------------------------------------------------------------------- */
  bufferRect.minx = p.x - buffer;
  bufferRect.maxx = p.x + buffer;
  bufferRect.miny = p.y - buffer;
  bufferRect.maxy = p.y + buffer;
1014

1015
  rlinfo->range_mode = mode;
1016

1017 1018 1019 1020
  if( maxresults != 0 ) {
    previous_maxresults = rlinfo->query_result_hard_max;
    rlinfo->query_result_hard_max = maxresults;
  }
1021

1022
  result = msRasterQueryByRect( map, layer, bufferRect );
1023

1024 1025
  if( rlinfo != NULL && maxresults != 0 )
    rlinfo->query_result_hard_max = previous_maxresults;
1026

1027
  return result;
1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043
#endif /* USE_GDAL  */
}

/************************************************************************/
/* ==================================================================== */
/*                          RASTER Query Layer                          */
/*                                                                      */
/*      The results of a raster query are treated as a set of shapes    */
/*      that can be accessed using the normal layer semantics.          */
/* ==================================================================== */
/************************************************************************/

/************************************************************************/
/*                         msRASTERLayerOpen()                          */
/************************************************************************/

1044
int msRASTERLayerOpen(layerObj *layer)
1045 1046
{
#ifndef USE_GDAL
1047 1048 1049 1050
  msSetError( MS_IMGERR,
              "Rasters queries only supported with GDAL support enabled.",
              "msRASTERLayerOpen()" );
  return MS_FAILURE;
1051
#else
1052
  rasterLayerInfo *rlinfo;
1053

1054 1055 1056
  /* If we don't have info, initialize an empty one now */
  if( layer->layerinfo == NULL )
    msRasterLayerInfoInitialize( layer );
1057

1058
  rlinfo = (rasterLayerInfo *) layer->layerinfo;
1059

1060 1061 1062
  rlinfo->refcount = rlinfo->refcount + 1;

  return MS_SUCCESS;
1063 1064 1065 1066 1067 1068 1069
#endif /* def USE_GDAL */
}

/************************************************************************/
/*                         msRASTERIsLayerOpen()                        */
/************************************************************************/

1070
int msRASTERLayerIsOpen(layerObj *layer)
1071 1072
{
#ifndef USE_GDAL
1073 1074 1075 1076
  msSetError( MS_IMGERR,
              "Rasters queries only supported with GDAL support enabled.",
              "msRasterLayerIsOpen()" );
  return MS_FALSE;
1077
#else
1078 1079 1080
  if (layer->layerinfo)
    return MS_TRUE;
  return MS_FALSE;
1081 1082 1083 1084 1085 1086 1087
#endif
}


/************************************************************************/
/*                     msRASTERLayerFreeItemInfo()                      */
/************************************************************************/
1088 1089
void msRASTERLayerFreeItemInfo(layerObj *layer)
{}
1090 1091 1092 1093 1094 1095 1096

/************************************************************************/
/*                     msRASTERLayerInitItemInfo()                      */
/*                                                                      */
/*      Perhaps we should be validating the requested items here?       */
/************************************************************************/

1097 1098 1099 1100
int msRASTERLayerInitItemInfo(layerObj *layer)
{
  return MS_SUCCESS;
}
1101 1102 1103 1104

/************************************************************************/
/*                      msRASTERLayerWhichShapes()                      */
/************************************************************************/
1105
int msRASTERLayerWhichShapes(layerObj *layer, rectObj rect, int isQuery)
1106 1107 1108

{
#ifndef USE_GDAL
1109
  return MS_FAILURE;
1110
#else
1111
  rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;
1112

1113 1114
  rlinfo->which_rect = rect;
  rlinfo->next_shape = 0;
1115

1116
  return MS_SUCCESS;
1117 1118 1119 1120 1121 1122 1123
#endif /* def USE_GDAL */
}

/************************************************************************/
/*                         msRASTERLayerClose()                         */
/************************************************************************/

1124
int msRASTERLayerClose(layerObj *layer)
1125 1126
{
#ifndef USE_GDAL
1127
  return MS_FAILURE;
1128
#else
1129
  rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;
1130

1131 1132
  if( rlinfo != NULL ) {
    rlinfo->refcount--;
1133

1134 1135 1136 1137
    if( rlinfo->refcount < 0 )
      msRasterLayerInfoFree( layer );
  }
  return MS_SUCCESS;
1138 1139 1140 1141 1142 1143 1144 1145 1146 1147
#endif /* def USE_GDAL */
}

/************************************************************************/
/*                       msRASTERLayerNextShape()                       */
/************************************************************************/

int msRASTERLayerNextShape(layerObj *layer, shapeObj *shape)
{
#ifndef USE_GDAL
1148
  return MS_FAILURE;
1149
#else
1150 1151 1152 1153 1154 1155 1156 1157 1158
  rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;

  if( rlinfo->next_shape < 0
      || rlinfo->next_shape >= rlinfo->query_results ) {
    msFreeShape(shape);
    shape->type = MS_SHAPE_NULL;
    return MS_DONE;
  } else {
    resultObj record;
1159

1160 1161 1162
    record.shapeindex = rlinfo->next_shape++;
    record.tileindex = 0;
    record.classindex = record.resultindex = -1;
1163

1164 1165
    return msRASTERLayerGetShape( layer, shape, &record);
  }
1166 1167 1168 1169 1170 1171 1172
#endif /* def USE_GDAL */
}

/************************************************************************/
/*                       msRASTERLayerGetShape()                        */
/************************************************************************/

1173
int msRASTERLayerGetShape(layerObj *layer, shapeObj *shape, resultObj *record)
1174 1175
{
#ifndef USE_GDAL
1176
  return MS_FAILURE;
1177
#else
1178 1179
  rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;
  int i;
1180

1181
  long shapeindex = record->shapeindex;
1182

1183 1184
  msFreeShape(shape);
  shape->type = MS_SHAPE_NULL;
1185

1186 1187 1188 1189 1190
  /* -------------------------------------------------------------------- */
  /*      Validate requested record id.                                   */
  /* -------------------------------------------------------------------- */
  if( shapeindex < 0 || shapeindex >= rlinfo->query_results ) {
    msSetError(MS_MISCERR,
1191
               "Out of range shape index requested.  Requested %ld\n"
1192 1193 1194 1195 1196
               "but only %d shapes available.",
               "msRASTERLayerGetShape()",
               shapeindex, rlinfo->query_results );
    return MS_FAILURE;
  }
1197

1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211
  /* -------------------------------------------------------------------- */
  /*      Apply the geometry.                                             */
  /* -------------------------------------------------------------------- */
  if( rlinfo->qc_x != NULL ) {
    lineObj line;
    pointObj point;

    shape->type = MS_SHAPE_POINT;

    line.numpoints = 1;
    line.point = &point;

    point.x = rlinfo->qc_x[shapeindex];
    point.y = rlinfo->qc_y[shapeindex];
1212
#ifdef USE_POINT_Z_M
1213
    point.m = 0.0;
1214 1215
#endif

1216 1217 1218
    msAddLine( shape, &line );
    msComputeBounds( shape );
  }
1219

1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246
  /* -------------------------------------------------------------------- */
  /*      Apply the requested items.                                      */
  /* -------------------------------------------------------------------- */
  if( layer->numitems > 0 ) {
    shape->values = (char **) msSmallMalloc(sizeof(char *) * layer->numitems);
    shape->numvalues = layer->numitems;

    for( i = 0; i < layer->numitems; i++ ) {
      const size_t bufferSize = 1000;
      char szWork[1000];

      szWork[0] = '\0';
      if( EQUAL(layer->items[i],"x") && rlinfo->qc_x_reproj )
        snprintf( szWork, bufferSize, "%.8g", rlinfo->qc_x_reproj[shapeindex] );
      else if( EQUAL(layer->items[i],"y") && rlinfo->qc_y_reproj )
        snprintf( szWork, bufferSize, "%.8g", rlinfo->qc_y_reproj[shapeindex] );

      else if( EQUAL(layer->items[i],"value_list") && rlinfo->qc_values ) {
        int iValue;

        for( iValue = 0; iValue < rlinfo->band_count; iValue++ ) {
          if( iValue != 0 )
            strlcat( szWork, ",", bufferSize);

          snprintf( szWork+strlen(szWork), bufferSize-strlen(szWork), "%.8g",
                    rlinfo->qc_values[shapeindex * rlinfo->band_count
                                      + iValue] );
1247
        }
1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267
      } else if( EQUALN(layer->items[i],"value_",6) && rlinfo->qc_values ) {
        int iValue = atoi(layer->items[i]+6);
        snprintf( szWork, bufferSize, "%.8g",
                  rlinfo->qc_values[shapeindex*rlinfo->band_count+iValue] );
      } else if( EQUAL(layer->items[i],"class") && rlinfo->qc_class ) {
        int p_class = rlinfo->qc_class[shapeindex];
        if( layer->class[p_class]->name != NULL )
          snprintf( szWork, bufferSize, "%.999s", layer->class[p_class]->name );
        else
          snprintf( szWork, bufferSize, "%d", p_class );
      } else if( EQUAL(layer->items[i],"red") && rlinfo->qc_red )
        snprintf( szWork, bufferSize, "%d", rlinfo->qc_red[shapeindex] );
      else if( EQUAL(layer->items[i],"green") && rlinfo->qc_green )
        snprintf( szWork, bufferSize, "%d", rlinfo->qc_green[shapeindex] );
      else if( EQUAL(layer->items[i],"blue") && rlinfo->qc_blue )
        snprintf( szWork, bufferSize, "%d", rlinfo->qc_blue[shapeindex] );
      else if( EQUAL(layer->items[i],"count") && rlinfo->qc_count )
        snprintf( szWork, bufferSize, "%d", rlinfo->qc_count[shapeindex] );

      shape->values[i] = msStrdup(szWork);
1268
    }
1269
  }
1270

1271 1272 1273 1274
  /* -------------------------------------------------------------------- */
  /*      Eventually we should likey apply the geometry properly but      */
  /*      we don't really care about the geometry for query purposes.     */
  /* -------------------------------------------------------------------- */
1275

1276
  return MS_SUCCESS;
1277 1278 1279 1280 1281 1282 1283 1284 1285 1286
#endif /* def USE_GDAL */
}

/************************************************************************/
/*                       msRASTERLayerGetItems()                        */
/************************************************************************/

int msRASTERLayerGetItems(layerObj *layer)
{
#ifndef USE_GDAL
1287
  return MS_FAILURE;
1288
#else
1289
  rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;
1290

1291 1292
  if( rlinfo == NULL )
    return MS_FAILURE;
1293

1294 1295 1296 1297 1298 1299 1300 1301 1302 1303 1304 1305 1306
  layer->items = (char **) msSmallCalloc(sizeof(char *),10);

  layer->numitems = 0;
  if( rlinfo->qc_x_reproj )
    layer->items[layer->numitems++] = msStrdup("x");
  if( rlinfo->qc_y_reproj )
    layer->items[layer->numitems++] = msStrdup("y");
  if( rlinfo->qc_values ) {
    int i;
    for( i = 0; i < rlinfo->band_count; i++ ) {
      char szName[100];
      snprintf( szName, sizeof(szName), "value_%d", i );
      layer->items[layer->numitems++] = msStrdup(szName);
1307
    }
1308 1309 1310 1311 1312 1313 1314 1315 1316 1317 1318 1319 1320 1321
    layer->items[layer->numitems++] = msStrdup("value_list");
  }
  if( rlinfo->qc_class )
    layer->items[layer->numitems++] = msStrdup("class");
  if( rlinfo->qc_red )
    layer->items[layer->numitems++] = msStrdup("red");
  if( rlinfo->qc_green )
    layer->items[layer->numitems++] = msStrdup("green");
  if( rlinfo->qc_blue )
    layer->items[layer->numitems++] = msStrdup("blue");
  if( rlinfo->qc_count )
    layer->items[layer->numitems++] = msStrdup("count");

  return msRASTERLayerInitItemInfo(layer);
1322 1323 1324 1325 1326 1327 1328 1329 1330
#endif /* def USE_GDAL */
}

/************************************************************************/
/*                       msRASTERLayerGetExtent()                       */
/************************************************************************/

int msRASTERLayerGetExtent(layerObj *layer, rectObj *extent)

1331
{
1332
#ifndef USE_GDAL
1333
  return MS_FAILURE;
1334 1335 1336 1337 1338 1339
#else
  char szPath[MS_MAXPATHLEN];
  mapObj *map = layer->map;
  double adfGeoTransform[6];
  int nXSize, nYSize;
  GDALDatasetH hDS;
1340
  shapefileObj *tileshpfile;
1341
  int tilelayerindex = -1;
1342
  CPLErr eErr = CE_Failure;
1343
  char *decrypted_path;
1344

1345 1346 1347
  if( (!layer->data || strlen(layer->data) == 0)
      && layer->tileindex == NULL) {
    /* should we be issuing a specific error about not supporting
1348 1349 1350 1351 1352 1353 1354
       extents for tileindexed raster layers? */
    return MS_FAILURE;
  }

  if( map == NULL )
    return MS_FAILURE;

1355
  /* If the layer use a tileindex, return the extent of the tileindex shapefile/referenced layer */
1356
  if (layer->tileindex) {
1357 1358 1359
    tilelayerindex = msGetLayerIndex(map, layer->tileindex);
    if(tilelayerindex != -1) /* does the tileindex reference another layer */
      return msLayerGetExtent(GET_LAYER(map, tilelayerindex), extent);
1360
    else {
1361 1362 1363 1364 1365 1366
      tileshpfile = (shapefileObj *) malloc(sizeof(shapefileObj));
      MS_CHECK_ALLOC(tileshpfile, sizeof(shapefileObj), MS_FAILURE);

      if(msShapefileOpen(tileshpfile, "rb", msBuildPath3(szPath, map->mappath, map->shapepath, layer->tileindex), MS_TRUE) == -1)
        if(msShapefileOpen(tileshpfile, "rb", msBuildPath(szPath, map->mappath, layer->tileindex), MS_TRUE) == -1)
          return MS_FAILURE;
1367

1368 1369 1370 1371 1372 1373 1374
      *extent = tileshpfile->bounds;
      msShapefileClose(tileshpfile);
      free(tileshpfile);
      return MS_SUCCESS;
    }
  }

1375
  msTryBuildPath3(szPath, map->mappath, map->shapepath, layer->data);
1376
  decrypted_path = msDecryptStringTokens( map, szPath );
1377 1378

  msAcquireLock( TLOCK_GDAL );
1379 1380 1381 1382 1383 1384 1385
  if( decrypted_path ) {
    hDS = GDALOpen(decrypted_path, GA_ReadOnly );
    msFree( decrypted_path );
  } else
    hDS = NULL;

  if( hDS != NULL ) {
1386 1387 1388
    nXSize = GDALGetRasterXSize( hDS );
    nYSize = GDALGetRasterYSize( hDS );
    eErr = GDALGetGeoTransform( hDS, adfGeoTransform );
1389

1390 1391
    GDALClose( hDS );
  }
1392

1393 1394
  msReleaseLock( TLOCK_GDAL );

1395 1396
  if( hDS == NULL || eErr != CE_None ) {
    return MS_FAILURE;
1397 1398 1399 1400
  }

  /* If this appears to be an ungeoreferenced raster than flip it for
     mapservers purposes. */
1401 1402 1403
  if( adfGeoTransform[5] == 1.0 && adfGeoTransform[3] == 0.0 ) {
    adfGeoTransform[5] = -1.0;
    adfGeoTransform[3] = nYSize;
1404
  }
1405

1406 1407
  extent->minx = adfGeoTransform[0];
  extent->maxy = adfGeoTransform[3];
1408

1409 1410
  extent->maxx = adfGeoTransform[0] + nXSize * adfGeoTransform[1];
  extent->miny = adfGeoTransform[3] + nYSize * adfGeoTransform[5];
1411

1412 1413 1414 1415 1416 1417 1418 1419 1420 1421 1422 1423 1424 1425 1426 1427 1428 1429 1430 1431
  return MS_SUCCESS;
#endif /* def USE_GDAL */
}


/************************************************************************/
/*                     msRASTERLayerSetTimeFilter()                     */
/*                                                                      */
/*      This function is actually just used in the context of           */
/*      setting a filter on the tileindex for time based queries.       */
/*      For instance via WMS requests.  So it isn't really related      */
/*      to the "raster query" support at all.                           */
/*                                                                      */
/*      If a local shapefile tileindex is in use, we will set a         */
/*      backtics filter (shapefile compatible).  If another layer is    */
/*      being used as the tileindex then we will forward the            */
/*      SetTimeFilter call to it.  If there is no tileindex in          */
/*      place, we do nothing.                                           */
/************************************************************************/

1432
int msRASTERLayerSetTimeFilter(layerObj *layer, const char *timestring,
1433 1434
                               const char *timefield)
{
1435
  int tilelayerindex;
1436

1437 1438 1439 1440 1441
  /* -------------------------------------------------------------------- */
  /*      If we don't have a tileindex the time filter has no effect.     */
  /* -------------------------------------------------------------------- */
  if( layer->tileindex == NULL )
    return MS_SUCCESS;
1442

1443 1444 1445 1446 1447 1448 1449 1450 1451 1452 1453 1454 1455 1456 1457 1458 1459 1460 1461 1462 1463 1464 1465
  /* -------------------------------------------------------------------- */
  /*      Find the tileindex layer.                                       */
  /* -------------------------------------------------------------------- */
  tilelayerindex = msGetLayerIndex(layer->map, layer->tileindex);

  /* -------------------------------------------------------------------- */
  /*      If we are using a local shapefile as our tileindex (that is     */
  /*      to say, the tileindex name is not of another layer), then we    */
  /*      just install a backtics style filter on the raster layer.       */
  /*      This is propogated to the "working layer" created for the       */
  /*      tileindex by code in mapraster.c.                               */
  /* -------------------------------------------------------------------- */
  if( tilelayerindex == -1 )
    return msLayerMakeBackticsTimeFilter( layer, timestring, timefield );

  /* -------------------------------------------------------------------- */
  /*      Otherwise we invoke the tileindex layers SetTimeFilter          */
  /*      method.                                                         */
  /* -------------------------------------------------------------------- */
  if ( msCheckParentPointer(layer->map,"map")==MS_FAILURE )
    return MS_FAILURE;
  return msLayerSetTimeFilter( layer->GET_LAYER(map,tilelayerindex),
                               timestring, timefield );
1466 1467 1468 1469 1470 1471 1472 1473 1474
}

/************************************************************************/
/*                msRASTERLayerInitializeVirtualTable()                 */
/************************************************************************/

int
msRASTERLayerInitializeVirtualTable(layerObj *layer)
{
1475 1476 1477 1478 1479 1480 1481 1482 1483 1484 1485 1486 1487 1488 1489 1490 1491 1492 1493 1494
  assert(layer != NULL);
  assert(layer->vtable != NULL);

  layer->vtable->LayerInitItemInfo = msRASTERLayerInitItemInfo;
  layer->vtable->LayerFreeItemInfo = msRASTERLayerFreeItemInfo;
  layer->vtable->LayerOpen = msRASTERLayerOpen;
  layer->vtable->LayerIsOpen = msRASTERLayerIsOpen;
  layer->vtable->LayerWhichShapes = msRASTERLayerWhichShapes;
  layer->vtable->LayerNextShape = msRASTERLayerNextShape;
  layer->vtable->LayerGetShape = msRASTERLayerGetShape;
  layer->vtable->LayerClose = msRASTERLayerClose;
  layer->vtable->LayerGetItems = msRASTERLayerGetItems;
  layer->vtable->LayerGetExtent = msRASTERLayerGetExtent;
  /* layer->vtable->LayerGetAutoStyle, use default */
  /* layer->vtable->LayerApplyFilterToLayer, use default */
  /* layer->vtable->LayerCloseConnection = msRASTERLayerClose; */
  /* we use backtics for proper tileindex shapefile functioning */
  layer->vtable->LayerSetTimeFilter = msRASTERLayerSetTimeFilter;
  /* layer->vtable->LayerCreateItems, use default */
  /* layer->vtable->LayerGetNumFeatures, use default */
1495

1496
  return MS_SUCCESS;
1497 1498
}