1 /******************************************************************************
2 * $Id$
3 *
4 * Project: MapServer
5 * Purpose: Implementation of query operations on rasters.
6 * Author: Frank Warmerdam, warmerdam@pobox.com
7 *
8 ******************************************************************************
9 * Copyright (c) 1996-2005 Regents of the University of Minnesota.
10 *
11 * Permission is hereby granted, free of charge, to any person obtaining a
12 * copy of this software and associated documentation files (the "Software"),
13 * to deal in the Software without restriction, including without limitation
14 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 * and/or sell copies of the Software, and to permit persons to whom the
16 * Software is furnished to do so, subject to the following conditions:
17 *
18 * The above copyright notice and this permission notice shall be included in
19 * all copies of this Software or works derived from this Software.
20 *
21 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 * DEALINGS IN THE SOFTWARE.
28 ****************************************************************************/
29
30 #include <assert.h>
31 #include "mapserver.h"
32 #include "mapresample.h"
33 #include "mapthread.h"
34 #include "mapraster.h"
35
36
37 int msRASTERLayerGetShape(layerObj *layer, shapeObj *shape, resultObj *record);
38 int msRASTERLayerGetItems(layerObj *layer);
39
40 /* ==================================================================== */
41 /* For now the rasterLayerInfo lives here since it is just used */
42 /* to hold information related to queries. */
43 /* ==================================================================== */
44 typedef struct {
45
46 /* query cache results */
47 int query_results;
48 int query_alloc_max;
49 int query_request_max;
50 int query_result_hard_max;
51 int raster_query_mode;
52 int band_count;
53
54 int refcount;
55
56 rectObj which_rect;
57 int next_shape;
58
59 double *qc_x;
60 double *qc_y;
61 double *qc_x_reproj;
62 double *qc_y_reproj;
63 float *qc_values;
64 int *qc_class;
65 int *qc_red;
66 int *qc_green;
67 int *qc_blue;
68 int *qc_count;
69 int *qc_tileindex;
70
71 /* query bound in force */
72 shapeObj *searchshape;
73
74 /* Only nearest result to this point. */
75 int range_mode; /* MS_QUERY_SINGLE, MS_QUERY_MULTIPLE or -1 (skip test) */
76 double range_dist;
77 pointObj target_point;
78
79 GDALColorTableH hCT;
80
81 double shape_tolerance;
82
83 } rasterLayerInfo;
84
85 #define RQM_UNKNOWN 0
86 #define RQM_ENTRY_PER_PIXEL 1
87 #define RQM_HIST_ON_CLASS 2
88 #define RQM_HIST_ON_VALUE 3
89
90 extern int InvGeoTransform( double *gt_in, double *gt_out );
91 #define GEO_TRANS(tr,x,y) ((tr)[0]+(tr)[1]*(x)+(tr)[2]*(y))
92
93 /************************************************************************/
94 /* addResult() */
95 /* */
96 /* this is a copy of the code in mapquery.c. Should we prepare */
97 /* a small public API for managing results caches? */
98 /************************************************************************/
99
addResult(resultCacheObj * cache,int classindex,int shapeindex,int tileindex)100 static int addResult(resultCacheObj *cache, int classindex, int shapeindex, int tileindex)
101 {
102 int i;
103
104 if(cache->numresults == cache->cachesize) { /* just add it to the end */
105 if(cache->cachesize == 0)
106 cache->results = (resultObj *) malloc(sizeof(resultObj)*MS_RESULTCACHEINCREMENT);
107 else
108 cache->results = (resultObj *) realloc(cache->results, sizeof(resultObj)*(cache->cachesize+MS_RESULTCACHEINCREMENT));
109 if(!cache->results) {
110 msSetError(MS_MEMERR, "Realloc() error.", "addResult()");
111 return(MS_FAILURE);
112 }
113 cache->cachesize += MS_RESULTCACHEINCREMENT;
114 }
115
116 i = cache->numresults;
117
118 cache->results[i].classindex = classindex;
119 cache->results[i].tileindex = tileindex;
120 cache->results[i].shapeindex = shapeindex;
121 cache->results[i].resultindex = -1; /* unused */
122 cache->results[i].shape = NULL;
123 cache->numresults++;
124
125 return(MS_SUCCESS);
126 }
127
128 /************************************************************************/
129 /* msRasterLayerInfoFree() */
130 /************************************************************************/
131
msRasterLayerInfoFree(layerObj * layer)132 static void msRasterLayerInfoFree( layerObj *layer )
133
134 {
135 rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;
136
137 if( rlinfo == NULL )
138 return;
139
140 if( rlinfo->qc_x != NULL ) {
141 free( rlinfo->qc_x );
142 free( rlinfo->qc_y );
143 free( rlinfo->qc_x_reproj );
144 free( rlinfo->qc_y_reproj );
145 }
146
147 if( rlinfo->qc_values )
148 free( rlinfo->qc_values );
149
150 if( rlinfo->qc_class ) {
151 free( rlinfo->qc_class );
152 }
153
154 if( rlinfo->qc_red ) {
155 free( rlinfo->qc_red );
156 free( rlinfo->qc_green );
157 free( rlinfo->qc_blue );
158 }
159
160 if( rlinfo->qc_count != NULL )
161 free( rlinfo->qc_count );
162
163 if( rlinfo->qc_tileindex != NULL )
164 free( rlinfo->qc_tileindex );
165
166 free( rlinfo );
167
168 layer->layerinfo = NULL;
169 }
170
171 /************************************************************************/
172 /* msRasterLayerInfoInitialize() */
173 /************************************************************************/
174
msRasterLayerInfoInitialize(layerObj * layer)175 static void msRasterLayerInfoInitialize( layerObj *layer )
176
177 {
178 rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;
179
180 if( rlinfo != NULL )
181 return;
182
183 rlinfo = (rasterLayerInfo *) msSmallCalloc(1,sizeof(rasterLayerInfo));
184 layer->layerinfo = rlinfo;
185
186 rlinfo->band_count = -1;
187 rlinfo->raster_query_mode = RQM_ENTRY_PER_PIXEL;
188 rlinfo->range_mode = -1; /* inactive */
189 rlinfo->refcount = 0;
190 rlinfo->shape_tolerance = 0.0;
191
192 /* We need to do this or the layer->layerinfo will be interpreted */
193 /* as shapefile access info because the default connectiontype is */
194 /* MS_SHAPEFILE. */
195 if (layer->connectiontype != MS_WMS && layer->connectiontype != MS_KERNELDENSITY)
196 layer->connectiontype = MS_RASTER;
197
198 rlinfo->query_result_hard_max = 1000000;
199
200 if( CSLFetchNameValue( layer->processing, "RASTER_QUERY_MAX_RESULT" )
201 != NULL ) {
202 rlinfo->query_result_hard_max =
203 atoi(CSLFetchNameValue( layer->processing, "RASTER_QUERY_MAX_RESULT" ));
204 }
205 }
206
207 /************************************************************************/
208 /* msRasterQueryAddPixel() */
209 /************************************************************************/
210
msRasterQueryAddPixel(layerObj * layer,pointObj * location,pointObj * reprojectedLocation,float * values)211 static void msRasterQueryAddPixel( layerObj *layer, pointObj *location,
212 pointObj *reprojectedLocation,
213 float *values )
214
215 {
216 rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;
217 int red = 0, green = 0, blue = 0, nodata = FALSE;
218 int p_class = -1;
219
220 if( rlinfo->query_results == rlinfo->query_result_hard_max )
221 return;
222
223 /* -------------------------------------------------------------------- */
224 /* Is this our first time in? If so, do an initial allocation */
225 /* for the data arrays suitable to our purposes. */
226 /* -------------------------------------------------------------------- */
227 if( rlinfo->query_alloc_max == 0 ) {
228 rlinfo->query_alloc_max = 2;
229
230 switch( rlinfo->raster_query_mode ) {
231 case RQM_ENTRY_PER_PIXEL:
232 rlinfo->qc_x = (double *)
233 msSmallCalloc(sizeof(double),rlinfo->query_alloc_max);
234 rlinfo->qc_y = (double *)
235 msSmallCalloc(sizeof(double),rlinfo->query_alloc_max);
236 rlinfo->qc_x_reproj = (double *)
237 msSmallCalloc(sizeof(double),rlinfo->query_alloc_max);
238 rlinfo->qc_y_reproj = (double *)
239 msSmallCalloc(sizeof(double),rlinfo->query_alloc_max);
240 rlinfo->qc_values = (float *)
241 msSmallCalloc(sizeof(float),
242 rlinfo->query_alloc_max*rlinfo->band_count);
243 rlinfo->qc_red = (int *)
244 msSmallCalloc(sizeof(int),rlinfo->query_alloc_max);
245 rlinfo->qc_green = (int *)
246 msSmallCalloc(sizeof(int),rlinfo->query_alloc_max);
247 rlinfo->qc_blue = (int *)
248 msSmallCalloc(sizeof(int),rlinfo->query_alloc_max);
249 if( layer->numclasses > 0 )
250 rlinfo->qc_class = (int *)
251 msSmallCalloc(sizeof(int),rlinfo->query_alloc_max);
252 break;
253
254 case RQM_HIST_ON_CLASS:
255 break;
256
257 case RQM_HIST_ON_VALUE:
258 break;
259
260 default:
261 assert( FALSE );
262 }
263 }
264
265 /* -------------------------------------------------------------------- */
266 /* Reallocate the data arrays larger if they are near the max */
267 /* now. */
268 /* -------------------------------------------------------------------- */
269 if( rlinfo->query_results == rlinfo->query_alloc_max ) {
270 rlinfo->query_alloc_max = rlinfo->query_alloc_max * 2 + 100;
271
272 if( rlinfo->qc_x != NULL )
273 rlinfo->qc_x = msSmallRealloc(rlinfo->qc_x,
274 sizeof(double) * rlinfo->query_alloc_max);
275 if( rlinfo->qc_y != NULL )
276 rlinfo->qc_y = msSmallRealloc(rlinfo->qc_y,
277 sizeof(double) * rlinfo->query_alloc_max);
278 if( rlinfo->qc_x_reproj != NULL )
279 rlinfo->qc_x_reproj = msSmallRealloc(rlinfo->qc_x_reproj,
280 sizeof(double) * rlinfo->query_alloc_max);
281 if( rlinfo->qc_y_reproj != NULL )
282 rlinfo->qc_y_reproj = msSmallRealloc(rlinfo->qc_y_reproj,
283 sizeof(double) * rlinfo->query_alloc_max);
284 if( rlinfo->qc_values != NULL )
285 rlinfo->qc_values =
286 msSmallRealloc(rlinfo->qc_values,
287 sizeof(float) * rlinfo->query_alloc_max
288 * rlinfo->band_count );
289 if( rlinfo->qc_class != NULL )
290 rlinfo->qc_class = msSmallRealloc(rlinfo->qc_class,
291 sizeof(int) * rlinfo->query_alloc_max);
292 if( rlinfo->qc_red != NULL )
293 rlinfo->qc_red = msSmallRealloc(rlinfo->qc_red,
294 sizeof(int) * rlinfo->query_alloc_max);
295 if( rlinfo->qc_green != NULL )
296 rlinfo->qc_green = msSmallRealloc(rlinfo->qc_green,
297 sizeof(int) * rlinfo->query_alloc_max);
298 if( rlinfo->qc_blue != NULL )
299 rlinfo->qc_blue = msSmallRealloc(rlinfo->qc_blue,
300 sizeof(int) * rlinfo->query_alloc_max);
301 if( rlinfo->qc_count != NULL )
302 rlinfo->qc_count = msSmallRealloc(rlinfo->qc_count,
303 sizeof(int) * rlinfo->query_alloc_max);
304 if( rlinfo->qc_tileindex != NULL )
305 rlinfo->qc_tileindex = msSmallRealloc(rlinfo->qc_tileindex,
306 sizeof(int) * rlinfo->query_alloc_max);
307 }
308
309 /* -------------------------------------------------------------------- */
310 /* Handle colormap */
311 /* -------------------------------------------------------------------- */
312 if( rlinfo->hCT != NULL ) {
313 int pct_index = (int) floor(values[0]);
314 GDALColorEntry sEntry;
315
316 if( GDALGetColorEntryAsRGB( rlinfo->hCT, pct_index, &sEntry ) ) {
317 red = sEntry.c1;
318 green = sEntry.c2;
319 blue = sEntry.c3;
320
321 if( sEntry.c4 == 0 )
322 nodata = TRUE;
323 } else
324 nodata = TRUE;
325 }
326
327 /* -------------------------------------------------------------------- */
328 /* Color derived from greyscale value. */
329 /* -------------------------------------------------------------------- */
330 else {
331 if( rlinfo->band_count >= 3 ) {
332 red = (int) MS_MAX(0,MS_MIN(255,values[0]));
333 green = (int) MS_MAX(0,MS_MIN(255,values[1]));
334 blue = (int) MS_MAX(0,MS_MIN(255,values[2]));
335 } else {
336 red = green = blue = (int) MS_MAX(0,MS_MIN(255,values[0]));
337 }
338 }
339
340 /* -------------------------------------------------------------------- */
341 /* Handle classification. */
342 /* */
343 /* NOTE: The following is really quite inadequate to deal with */
344 /* classifications based on [red], [green] and [blue] as */
345 /* described in: */
346 /* http://mapserver.gis.umn.edu/bugs/show_bug.cgi?id=1021 */
347 /* -------------------------------------------------------------------- */
348 if( rlinfo->qc_class != NULL ) {
349 p_class = msGetClass_FloatRGB(layer, values[0],
350 red, green, blue );
351
352 if( p_class == -1 )
353 nodata = TRUE;
354 else {
355 nodata = FALSE;
356 rlinfo->qc_class[rlinfo->query_results] = p_class;
357 if( layer->class[p_class]->numstyles > 0 ) {
358 red = layer->class[p_class]->styles[0]->color.red;
359 green = layer->class[p_class]->styles[0]->color.green;
360 blue = layer->class[p_class]->styles[0]->color.blue;
361 } else {
362 red = green = blue = 0;
363 }
364 }
365 }
366
367 /* -------------------------------------------------------------------- */
368 /* Record the color. */
369 /* -------------------------------------------------------------------- */
370 rlinfo->qc_red[rlinfo->query_results] = red;
371 rlinfo->qc_green[rlinfo->query_results] = green;
372 rlinfo->qc_blue[rlinfo->query_results] = blue;
373
374 /* -------------------------------------------------------------------- */
375 /* Record spatial location. */
376 /* -------------------------------------------------------------------- */
377 if( rlinfo->qc_x != NULL ) {
378 rlinfo->qc_x[rlinfo->query_results] = location->x;
379 rlinfo->qc_y[rlinfo->query_results] = location->y;
380 rlinfo->qc_x_reproj[rlinfo->query_results] = reprojectedLocation->x;
381 rlinfo->qc_y_reproj[rlinfo->query_results] = reprojectedLocation->y;
382 }
383
384 /* -------------------------------------------------------------------- */
385 /* Record actual pixel value(s). */
386 /* -------------------------------------------------------------------- */
387 if( rlinfo->qc_values != NULL )
388 memcpy( rlinfo->qc_values + rlinfo->query_results * rlinfo->band_count,
389 values, sizeof(float) * rlinfo->band_count );
390
391 /* -------------------------------------------------------------------- */
392 /* Add to the results cache. */
393 /* -------------------------------------------------------------------- */
394 if( ! nodata ) {
395 addResult( layer->resultcache, p_class, rlinfo->query_results, 0 );
396 rlinfo->query_results++;
397 }
398 }
399
400 /************************************************************************/
401 /* msRasterQueryByRectLow() */
402 /************************************************************************/
403
404 static int
msRasterQueryByRectLow(mapObj * map,layerObj * layer,GDALDatasetH hDS,rectObj queryRect)405 msRasterQueryByRectLow(mapObj *map, layerObj *layer, GDALDatasetH hDS,
406 rectObj queryRect)
407
408 {
409 double adfGeoTransform[6], adfInvGeoTransform[6];
410 double dfXMin, dfYMin, dfXMax, dfYMax, dfX, dfY, dfAdjustedRange;
411 int nWinXOff, nWinYOff, nWinXSize, nWinYSize;
412 int nRXSize, nRYSize;
413 float *pafRaster;
414 int nBandCount, *panBandMap, iPixel, iLine;
415 CPLErr eErr;
416 rasterLayerInfo *rlinfo;
417 rectObj searchrect;
418 #if PROJ_VERSION_MAJOR < 6
419 int mayNeedLonWrapAdjustment = MS_FALSE;
420 #endif
421
422 rlinfo = (rasterLayerInfo *) layer->layerinfo;
423
424 /* -------------------------------------------------------------------- */
425 /* Reproject the search rect into the projection of this */
426 /* layer/file if needed. */
427 /* -------------------------------------------------------------------- */
428 searchrect = queryRect;
429 layer->project = msProjectionsDiffer(&(layer->projection), &(map->projection));
430 if(layer->project)
431 msProjectRect(&(map->projection), &(layer->projection), &searchrect);
432
433 /* -------------------------------------------------------------------- */
434 /* Transform the rectangle in target ground coordinates to */
435 /* pixel/line extents on the file. Process all 4 corners, to */
436 /* build extents. */
437 /* -------------------------------------------------------------------- */
438 nRXSize = GDALGetRasterXSize( hDS );
439 nRYSize = GDALGetRasterYSize( hDS );
440
441 msGetGDALGeoTransform( hDS, map, layer, adfGeoTransform );
442 InvGeoTransform( adfGeoTransform, adfInvGeoTransform );
443
444 /* top left */
445 dfXMin = dfXMax = GEO_TRANS(adfInvGeoTransform,
446 searchrect.minx, searchrect.maxy);
447 dfYMin = dfYMax = GEO_TRANS(adfInvGeoTransform+3,
448 searchrect.minx, searchrect.maxy);
449
450 /* top right */
451 dfX = GEO_TRANS(adfInvGeoTransform , searchrect.maxx, searchrect.maxy);
452 dfY = GEO_TRANS(adfInvGeoTransform+3, searchrect.maxx, searchrect.maxy);
453 dfXMin = MS_MIN(dfXMin,dfX);
454 dfXMax = MS_MAX(dfXMax,dfX);
455 dfYMin = MS_MIN(dfYMin,dfY);
456 dfYMax = MS_MAX(dfYMax,dfY);
457
458 /* bottom left */
459 dfX = GEO_TRANS(adfInvGeoTransform , searchrect.minx, searchrect.miny);
460 dfY = GEO_TRANS(adfInvGeoTransform+3, searchrect.minx, searchrect.miny);
461 dfXMin = MS_MIN(dfXMin,dfX);
462 dfXMax = MS_MAX(dfXMax,dfX);
463 dfYMin = MS_MIN(dfYMin,dfY);
464 dfYMax = MS_MAX(dfYMax,dfY);
465
466 /* bottom right */
467 dfX = GEO_TRANS(adfInvGeoTransform , searchrect.maxx, searchrect.miny);
468 dfY = GEO_TRANS(adfInvGeoTransform+3, searchrect.maxx, searchrect.miny);
469 dfXMin = MS_MIN(dfXMin,dfX);
470 dfXMax = MS_MAX(dfXMax,dfX);
471 dfYMin = MS_MIN(dfYMin,dfY);
472 dfYMax = MS_MAX(dfYMax,dfY);
473
474 /* -------------------------------------------------------------------- */
475 /* Trim the rectangle to the area of the file itself, but out */
476 /* to the edges of the touched edge pixels. */
477 /* -------------------------------------------------------------------- */
478 dfXMin = MS_MAX(0.0,MS_MIN(nRXSize,floor(dfXMin)));
479 dfYMin = MS_MAX(0.0,MS_MIN(nRYSize,floor(dfYMin)));
480 dfXMax = MS_MAX(0.0,MS_MIN(nRXSize,ceil(dfXMax)));
481 dfYMax = MS_MAX(0.0,MS_MIN(nRYSize,ceil(dfYMax)));
482
483 /* -------------------------------------------------------------------- */
484 /* Convert to integer offset/size values. */
485 /* -------------------------------------------------------------------- */
486 nWinXOff = (int) dfXMin;
487 nWinYOff = (int) dfYMin;
488 nWinXSize = (int) (dfXMax - dfXMin);
489 nWinYSize = (int) (dfYMax - dfYMin);
490
491 /* -------------------------------------------------------------------- */
492 /* What bands are we operating on? */
493 /* -------------------------------------------------------------------- */
494 panBandMap = msGetGDALBandList( layer, hDS, 0, &nBandCount );
495
496 if( rlinfo->band_count == -1 )
497 rlinfo->band_count = nBandCount;
498
499 if( nBandCount != rlinfo->band_count ) {
500 msSetError( MS_IMGERR,
501 "Got %d bands, but expected %d bands.",
502 "msRasterQueryByRectLow()",
503 nBandCount, rlinfo->band_count );
504
505 return -1;
506 }
507
508 /* -------------------------------------------------------------------- */
509 /* Try to load the raster data. For now we just load the first */
510 /* band in the file. Later we will deal with the various band */
511 /* selection criteria. */
512 /* -------------------------------------------------------------------- */
513 pafRaster = (float *)
514 calloc(sizeof(float),nWinXSize*nWinYSize*nBandCount);
515 MS_CHECK_ALLOC(pafRaster, sizeof(float)*nWinXSize*nWinYSize*nBandCount, -1);
516
517 eErr = GDALDatasetRasterIO( hDS, GF_Read,
518 nWinXOff, nWinYOff, nWinXSize, nWinYSize,
519 pafRaster, nWinXSize, nWinYSize, GDT_Float32,
520 nBandCount, panBandMap,
521 4 * nBandCount,
522 4 * nBandCount * nWinXSize,
523 4 );
524
525 if( eErr != CE_None ) {
526 msSetError( MS_IOERR, "GDALDatasetRasterIO() failed: %s",
527 "msRasterQueryByRectLow()", CPLGetLastErrorMsg() );
528
529 free( pafRaster );
530 return -1;
531 }
532
533 /* -------------------------------------------------------------------- */
534 /* Fetch color table for intepreting colors if needed. */
535 /* -------------------------------------------------------------------- */
536 rlinfo->hCT = GDALGetRasterColorTable(
537 GDALGetRasterBand( hDS, panBandMap[0] ) );
538
539 free( panBandMap );
540
541 /* -------------------------------------------------------------------- */
542 /* When computing whether pixels are within range we do it */
543 /* based on the center of the pixel to the target point but */
544 /* really it ought to be the nearest point on the pixel. It */
545 /* would be too much trouble to do this rigerously, so we just */
546 /* add a fudge factor so that a range of zero will find the */
547 /* pixel the target falls in at least. */
548 /* -------------------------------------------------------------------- */
549 dfAdjustedRange =
550 sqrt(adfGeoTransform[1] * adfGeoTransform[1]
551 + adfGeoTransform[2] * adfGeoTransform[2]
552 + adfGeoTransform[4] * adfGeoTransform[4]
553 + adfGeoTransform[5] * adfGeoTransform[5]) * 0.5 * 1.41421356237
554 + sqrt( rlinfo->range_dist );
555 dfAdjustedRange = dfAdjustedRange * dfAdjustedRange;
556
557 #if PROJ_VERSION_MAJOR < 6
558 if( layer->project &&
559 msProjIsGeographicCRS(&(layer->projection)) &&
560 msProjIsGeographicCRS(&(map->projection)) )
561 {
562 double dfLonWrap = 0;
563 mayNeedLonWrapAdjustment = msProjectHasLonWrap(&(layer->projection), &dfLonWrap);
564 }
565 #endif
566
567 reprojectionObj* reprojector = NULL;
568 if( layer->project )
569 {
570 reprojector = msProjectCreateReprojector(&(layer->projection), &(map->projection));
571 }
572
573 /* -------------------------------------------------------------------- */
574 /* Loop over all pixels determining which are "in". */
575 /* -------------------------------------------------------------------- */
576 for( iLine = 0; iLine < nWinYSize; iLine++ ) {
577 for( iPixel = 0; iPixel < nWinXSize; iPixel++ ) {
578 pointObj sPixelLocation,sReprojectedPixelLocation;
579
580 if( rlinfo->query_results == rlinfo->query_result_hard_max )
581 break;
582
583 /* transform pixel/line to georeferenced */
584 sPixelLocation.x =
585 GEO_TRANS(adfGeoTransform,
586 iPixel + nWinXOff + 0.5, iLine + nWinYOff + 0.5 );
587 sPixelLocation.y =
588 GEO_TRANS(adfGeoTransform+3,
589 iPixel + nWinXOff + 0.5, iLine + nWinYOff + 0.5 );
590
591 /* If projections differ, convert this back into the map */
592 /* projection for distance testing, and comprison to the */
593 /* search shape. Save the original pixel location coordinates */
594 /* in sPixelLocationInLayerSRS, so that we can return those */
595 /* coordinates if we have a hit */
596 sReprojectedPixelLocation = sPixelLocation;
597 if( reprojector )
598 {
599 #if PROJ_VERSION_MAJOR < 6
600 /* Works around a bug in PROJ < 6 when reprojecting from a lon_wrap */
601 /* geogCRS to a geogCRS, and the input abs(longitude) is > 180. Then */
602 /* lon_wrap was ignored and the output longitude remained as the source */
603 if( mayNeedLonWrapAdjustment )
604 {
605 if( rlinfo->target_point.x < sReprojectedPixelLocation.x - 180 )
606 sReprojectedPixelLocation.x -= 360;
607 else if( rlinfo->target_point.x > sReprojectedPixelLocation.x + 180 )
608 sReprojectedPixelLocation.x += 360;
609 }
610 #endif
611
612 msProjectPointEx( reprojector, &sReprojectedPixelLocation);
613 }
614
615 /* If we are doing QueryByShape, check against the shape now */
616 if( rlinfo->searchshape != NULL ) {
617 if( rlinfo->shape_tolerance == 0.0
618 && rlinfo->searchshape->type == MS_SHAPE_POLYGON ) {
619 if( msIntersectPointPolygon(
620 &sReprojectedPixelLocation, rlinfo->searchshape ) == MS_FALSE )
621 continue;
622 } else {
623 shapeObj tempShape;
624 lineObj tempLine;
625
626 memset( &tempShape, 0, sizeof(shapeObj) );
627 tempShape.type = MS_SHAPE_POINT;
628 tempShape.numlines = 1;
629 tempShape.line = &tempLine;
630 tempLine.numpoints = 1;
631 tempLine.point = &sReprojectedPixelLocation;
632
633 if( msDistanceShapeToShape(rlinfo->searchshape, &tempShape)
634 > rlinfo->shape_tolerance )
635 continue;
636 }
637 }
638
639 if( rlinfo->range_mode >= 0 ) {
640 double dist;
641
642 dist = (rlinfo->target_point.x - sReprojectedPixelLocation.x)
643 * (rlinfo->target_point.x - sReprojectedPixelLocation.x)
644 + (rlinfo->target_point.y - sReprojectedPixelLocation.y)
645 * (rlinfo->target_point.y - sReprojectedPixelLocation.y);
646
647 if( dist >= dfAdjustedRange )
648 continue;
649
650 /* If we can only have one feature, trim range and clear */
651 /* previous result. */
652 if( rlinfo->range_mode == MS_QUERY_SINGLE ) {
653 rlinfo->range_dist = dist;
654 rlinfo->query_results = 0;
655 }
656 }
657
658 msRasterQueryAddPixel( layer,
659 &sPixelLocation, // return coords in layer SRS
660 &sReprojectedPixelLocation,
661 pafRaster
662 + (iLine*nWinXSize + iPixel) * nBandCount );
663 }
664 }
665
666 /* -------------------------------------------------------------------- */
667 /* Cleanup. */
668 /* -------------------------------------------------------------------- */
669 free( pafRaster );
670 msProjectDestroyReprojector( reprojector );
671
672 return MS_SUCCESS;
673 }
674
675 /************************************************************************/
676 /* msRasterQueryByRect() */
677 /************************************************************************/
678
msRasterQueryByRect(mapObj * map,layerObj * layer,rectObj queryRect)679 int msRasterQueryByRect(mapObj *map, layerObj *layer, rectObj queryRect)
680
681 {
682 int status = MS_SUCCESS;
683 char *filename=NULL;
684
685 layerObj *tlp=NULL; /* pointer to the tile layer either real or temporary */
686 int tileitemindex=-1, tilelayerindex=-1, tilesrsindex=-1;
687 shapeObj tshp;
688 char tilename[MS_PATH_LENGTH], tilesrsname[1024];
689 int done, destroy_on_failure;
690
691 char szPath[MS_MAXPATHLEN];
692 rectObj searchrect;
693 rasterLayerInfo *rlinfo = NULL;
694
695 /* -------------------------------------------------------------------- */
696 /* Get the layer info. */
697 /* -------------------------------------------------------------------- */
698 if(!layer->layerinfo) {
699 msRasterLayerInfoInitialize( layer );
700 destroy_on_failure = 1;
701 } else {
702 destroy_on_failure = 0;
703 }
704 rlinfo = (rasterLayerInfo *) layer->layerinfo;
705
706 /* -------------------------------------------------------------------- */
707 /* Clear old results cache. */
708 /* -------------------------------------------------------------------- */
709 if(layer->resultcache) {
710 if(layer->resultcache->results) free(layer->resultcache->results);
711 free(layer->resultcache);
712 layer->resultcache = NULL;
713 }
714
715 /* -------------------------------------------------------------------- */
716 /* Initialize the results cache. */
717 /* -------------------------------------------------------------------- */
718 layer->resultcache = (resultCacheObj *)msSmallMalloc(sizeof(resultCacheObj));
719 layer->resultcache->results = NULL;
720 layer->resultcache->numresults = layer->resultcache->cachesize = 0;
721 layer->resultcache->bounds.minx =
722 layer->resultcache->bounds.miny =
723 layer->resultcache->bounds.maxx =
724 layer->resultcache->bounds.maxy = -1;
725
726 /* -------------------------------------------------------------------- */
727 /* Check if we should really be acting on this layer and */
728 /* provide debug info in various cases. */
729 /* -------------------------------------------------------------------- */
730 if(layer->debug > 0 || map->debug > 1)
731 msDebug( "msRasterQueryByRect(%s): entering.\n", layer->name );
732
733 if(!layer->data && !layer->tileindex) {
734 if(layer->debug > 0 || map->debug > 0 )
735 msDebug( "msRasterQueryByRect(%s): layer data and tileindex NULL ... doing nothing.", layer->name );
736 return MS_SUCCESS;
737 }
738
739 if((layer->status != MS_ON) && layer->status != MS_DEFAULT ) {
740 if(layer->debug > 0 )
741 msDebug( "msRasterQueryByRect(%s): not status ON or DEFAULT, doing nothing.", layer->name );
742 return MS_SUCCESS;
743 }
744
745 /* ==================================================================== */
746 /* Handle setting up tileindex layer. */
747 /* ==================================================================== */
748 if(layer->tileindex) { /* we have an index file */
749 msInitShape(&tshp);
750 searchrect = queryRect;
751
752 status = msDrawRasterSetupTileLayer(map, layer,
753 &searchrect,
754 MS_TRUE,
755 &tilelayerindex,
756 &tileitemindex,
757 &tilesrsindex,
758 &tlp);
759 if (status != MS_SUCCESS) {
760 goto cleanup;
761 }
762 } else {
763 /* we have to manually apply to scaletoken logic as we're not going through msLayerOpen() */
764 status = msLayerApplyScaletokens(layer,map->scaledenom);
765 if (status != MS_SUCCESS) {
766 goto cleanup;
767 }
768 }
769
770 /* -------------------------------------------------------------------- */
771 /* Iterate over all tiles (just one in untiled case). */
772 /* -------------------------------------------------------------------- */
773 done = MS_FALSE;
774 while( done == MS_FALSE && status == MS_SUCCESS ) {
775
776 GDALDatasetH hDS;
777 char *decrypted_path = NULL;
778
779 /* -------------------------------------------------------------------- */
780 /* Get filename. */
781 /* -------------------------------------------------------------------- */
782 if(layer->tileindex) {
783 status = msDrawRasterIterateTileIndex(layer, tlp, &tshp,
784 tileitemindex, tilesrsindex,
785 tilename, sizeof(tilename),
786 tilesrsname, sizeof(tilesrsname));
787 if( status == MS_FAILURE) {
788 break;
789 }
790
791 if(status == MS_DONE) break; /* no more tiles/images */
792 filename = tilename;
793 } else {
794 filename = layer->data;
795 done = MS_TRUE; /* only one image so we're done after this */
796 }
797
798 if(strlen(filename) == 0) continue;
799
800 /* -------------------------------------------------------------------- */
801 /* Open the file. */
802 /* -------------------------------------------------------------------- */
803 msGDALInitialize();
804
805 msDrawRasterBuildRasterPath(map, layer, filename, szPath);
806
807 decrypted_path = msDecryptStringTokens( map, szPath );
808 if( !decrypted_path ) {
809 status = MS_FAILURE;
810 goto cleanup;
811 }
812
813 msAcquireLock( TLOCK_GDAL );
814 if( !layer->tileindex )
815 {
816 char** connectionoptions = msGetStringListFromHashTable(&(layer->connectionoptions));
817 hDS = GDALOpenEx(decrypted_path,
818 GDAL_OF_RASTER,
819 NULL,
820 (const char* const*)connectionoptions,
821 NULL);
822 CSLDestroy(connectionoptions);
823 }
824 else
825 {
826 hDS = GDALOpen(decrypted_path, GA_ReadOnly );
827 }
828
829 if( hDS == NULL ) {
830 int ignore_missing = msMapIgnoreMissingData( map );
831 const char *cpl_error_msg = msDrawRasterGetCPLErrorMsg(decrypted_path, szPath);
832
833 msFree( decrypted_path );
834 decrypted_path = NULL;
835
836 msReleaseLock( TLOCK_GDAL );
837
838 if ( ignore_missing == MS_MISSING_DATA_FAIL ) {
839 if( layer->debug || map->debug )
840 msSetError( MS_IMGERR,
841 "Unable to open file %s for layer %s ... fatal error.\n%s",
842 "msRasterQueryByRect()",
843 szPath, layer->name, cpl_error_msg);
844 status = MS_FAILURE;
845 goto cleanup;
846 }
847 if( ignore_missing == MS_MISSING_DATA_LOG ) {
848 if( layer->debug || map->debug )
849 msDebug( "Unable to open file %s for layer %s ... ignoring this missing data.\n%s",
850 filename, layer->name, cpl_error_msg );
851 }
852 continue;
853 }
854
855 msFree( decrypted_path );
856 decrypted_path = NULL;
857
858 if( msDrawRasterLoadProjection(layer, hDS, filename, tilesrsindex, tilesrsname) != MS_SUCCESS )
859 {
860 msReleaseLock( TLOCK_GDAL );
861 status = MS_FAILURE;
862 goto cleanup;
863 }
864
865 /* -------------------------------------------------------------------- */
866 /* Perform actual query on this file. */
867 /* -------------------------------------------------------------------- */
868 if( status == MS_SUCCESS )
869 status = msRasterQueryByRectLow( map, layer, hDS, queryRect );
870
871 GDALClose( hDS );
872 msReleaseLock( TLOCK_GDAL );
873
874 } /* next tile */
875
876 /* -------------------------------------------------------------------- */
877 /* Cleanup tileindex if it is open. */
878 /* -------------------------------------------------------------------- */
879 cleanup:
880 if(layer->tileindex) { /* tiling clean-up */
881 msDrawRasterCleanupTileLayer(tlp, tilelayerindex);
882 } else {
883 msLayerRestoreFromScaletokens(layer);
884 }
885
886 /* -------------------------------------------------------------------- */
887 /* On failure, or empty result set, cleanup the rlinfo since we */
888 /* likely won't ever have it accessed or cleaned up later. */
889 /* -------------------------------------------------------------------- */
890 if( status == MS_FAILURE || rlinfo->query_results == 0 ) {
891 if(destroy_on_failure) {
892 msRasterLayerInfoFree( layer );
893 }
894 } else {
895 /* populate the items/numitems layer-level values */
896 msRASTERLayerGetItems(layer);
897 }
898
899 return status;
900 }
901
902 /************************************************************************/
903 /* msRasterQueryByShape() */
904 /************************************************************************/
905
msRasterQueryByShape(mapObj * map,layerObj * layer,shapeObj * selectshape)906 int msRasterQueryByShape(mapObj *map, layerObj *layer, shapeObj *selectshape)
907
908 {
909 rasterLayerInfo *rlinfo = NULL;
910 int status;
911 double tolerance;
912 rectObj searchrect;
913
914 msRasterLayerInfoInitialize( layer );
915
916 rlinfo = (rasterLayerInfo *) layer->layerinfo;
917
918 /* If the selection shape is point or line we use the default tolerance of
919 3, but for polygons we require an exact hit. */
920 if(layer->tolerance == -1) {
921 if(selectshape->type == MS_SHAPE_POINT ||
922 selectshape->type == MS_SHAPE_LINE)
923 tolerance = 3;
924 else
925 tolerance = 0;
926 } else
927 tolerance = layer->tolerance;
928
929 if(layer->toleranceunits == MS_PIXELS)
930 tolerance = tolerance
931 * msAdjustExtent(&(map->extent), map->width, map->height);
932 else
933 tolerance = tolerance
934 * (msInchesPerUnit(layer->toleranceunits,0)
935 / msInchesPerUnit(map->units,0));
936
937 rlinfo->searchshape = selectshape;
938 rlinfo->shape_tolerance = tolerance;
939
940 msComputeBounds( selectshape );
941 searchrect = selectshape->bounds;
942
943 searchrect.minx -= tolerance; /* expand the search box to account for layer tolerances (e.g. buffered searches) */
944 searchrect.maxx += tolerance;
945 searchrect.miny -= tolerance;
946 searchrect.maxy += tolerance;
947
948 status = msRasterQueryByRect( map, layer, searchrect );
949
950 rlinfo = (rasterLayerInfo *) layer->layerinfo;
951 if( rlinfo )
952 rlinfo->searchshape = NULL;
953
954 return status;
955 }
956
957 /************************************************************************/
958 /* msRasterQueryByPoint() */
959 /************************************************************************/
960
msRasterQueryByPoint(mapObj * map,layerObj * layer,int mode,pointObj p,double buffer,int maxresults)961 int msRasterQueryByPoint(mapObj *map, layerObj *layer, int mode, pointObj p,
962 double buffer, int maxresults)
963 {
964 int result;
965 double layer_tolerance;
966 rectObj bufferRect;
967 rasterLayerInfo *rlinfo = NULL;
968
969 msRasterLayerInfoInitialize( layer );
970
971 rlinfo = (rasterLayerInfo *) layer->layerinfo;
972
973 /* -------------------------------------------------------------------- */
974 /* If the buffer is not set, then use layer tolerances */
975 /* instead. The "buffer" distince is now in georeferenced */
976 /* units. Note that tolerances in pixels are basically map */
977 /* display pixels, not underlying raster pixels. It isn't */
978 /* clear that there is any way of requesting a buffer size in */
979 /* underlying pixels. */
980 /* -------------------------------------------------------------------- */
981 if(buffer <= 0) { /* use layer tolerance */
982 if(layer->tolerance == -1)
983 layer_tolerance = 3;
984 else
985 layer_tolerance = layer->tolerance;
986
987 if(layer->toleranceunits == MS_PIXELS)
988 buffer = layer_tolerance
989 * msAdjustExtent(&(map->extent), map->width, map->height);
990 else
991 buffer = layer_tolerance
992 * (msInchesPerUnit(layer->toleranceunits,0)
993 / msInchesPerUnit(map->units,0));
994 }
995
996 /* -------------------------------------------------------------------- */
997 /* Setup target point information, at this point they are in */
998 /* map coordinates. */
999 /* -------------------------------------------------------------------- */
1000 rlinfo->range_dist = buffer * buffer;
1001 rlinfo->target_point = p;
1002
1003 /* -------------------------------------------------------------------- */
1004 /* if we are in the MS_QUERY_SINGLE mode, first try a query with */
1005 /* zero tolerance. If this gets a raster pixel then we can be */
1006 /* reasonably assured that it is the closest to the query */
1007 /* point. This will potentially be must more efficient than */
1008 /* processing all pixels within the tolerance. */
1009 /* -------------------------------------------------------------------- */
1010 if( mode == MS_QUERY_SINGLE ) {
1011 rectObj pointRect;
1012
1013 pointRect.minx = p.x;
1014 pointRect.maxx = p.x;
1015 pointRect.miny = p.y;
1016 pointRect.maxy = p.y;
1017
1018 rlinfo->range_mode = MS_QUERY_SINGLE;
1019 result = msRasterQueryByRect( map, layer, pointRect );
1020 if( rlinfo->query_results > 0 )
1021 return result;
1022 }
1023
1024 /* -------------------------------------------------------------------- */
1025 /* Setup a rectangle that is everything within the designated */
1026 /* range and do a search against that. */
1027 /* -------------------------------------------------------------------- */
1028 bufferRect.minx = p.x - buffer;
1029 bufferRect.maxx = p.x + buffer;
1030 bufferRect.miny = p.y - buffer;
1031 bufferRect.maxy = p.y + buffer;
1032
1033 rlinfo->range_mode = mode;
1034
1035 if( maxresults != 0 ) {
1036 const int previous_maxresults = rlinfo->query_result_hard_max;
1037 rlinfo->query_result_hard_max = maxresults;
1038 result = msRasterQueryByRect( map, layer, bufferRect );
1039 rlinfo->query_result_hard_max = previous_maxresults;
1040 } else {
1041 result = msRasterQueryByRect( map, layer, bufferRect );
1042 }
1043
1044 return result;
1045 }
1046
1047 /************************************************************************/
1048 /* ==================================================================== */
1049 /* RASTER Query Layer */
1050 /* */
1051 /* The results of a raster query are treated as a set of shapes */
1052 /* that can be accessed using the normal layer semantics. */
1053 /* ==================================================================== */
1054 /************************************************************************/
1055
1056 /************************************************************************/
1057 /* msRASTERLayerOpen() */
1058 /************************************************************************/
1059
msRASTERLayerOpen(layerObj * layer)1060 int msRASTERLayerOpen(layerObj *layer)
1061 {
1062 rasterLayerInfo *rlinfo;
1063
1064 /* If we don't have info, initialize an empty one now */
1065 if( layer->layerinfo == NULL )
1066 msRasterLayerInfoInitialize( layer );
1067
1068 rlinfo = (rasterLayerInfo *) layer->layerinfo;
1069
1070 rlinfo->refcount = rlinfo->refcount + 1;
1071
1072 return MS_SUCCESS;
1073 }
1074
1075 /************************************************************************/
1076 /* msRASTERIsLayerOpen() */
1077 /************************************************************************/
1078
msRASTERLayerIsOpen(layerObj * layer)1079 int msRASTERLayerIsOpen(layerObj *layer)
1080 {
1081 if (layer->layerinfo)
1082 return MS_TRUE;
1083 return MS_FALSE;
1084 }
1085
1086
1087 /************************************************************************/
1088 /* msRASTERLayerFreeItemInfo() */
1089 /************************************************************************/
msRASTERLayerFreeItemInfo(layerObj * layer)1090 void msRASTERLayerFreeItemInfo(layerObj *layer)
1091 {}
1092
1093 /************************************************************************/
1094 /* msRASTERLayerInitItemInfo() */
1095 /* */
1096 /* Perhaps we should be validating the requested items here? */
1097 /************************************************************************/
1098
msRASTERLayerInitItemInfo(layerObj * layer)1099 int msRASTERLayerInitItemInfo(layerObj *layer)
1100 {
1101 return MS_SUCCESS;
1102 }
1103
1104 /************************************************************************/
1105 /* msRASTERLayerWhichShapes() */
1106 /************************************************************************/
msRASTERLayerWhichShapes(layerObj * layer,rectObj rect,int isQuery)1107 int msRASTERLayerWhichShapes(layerObj *layer, rectObj rect, int isQuery)
1108
1109 {
1110 rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;
1111
1112 rlinfo->which_rect = rect;
1113 rlinfo->next_shape = 0;
1114
1115 return MS_SUCCESS;
1116 }
1117
1118 /************************************************************************/
1119 /* msRASTERLayerClose() */
1120 /************************************************************************/
1121
msRASTERLayerClose(layerObj * layer)1122 int msRASTERLayerClose(layerObj *layer)
1123 {
1124 rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;
1125
1126 if( rlinfo != NULL ) {
1127 rlinfo->refcount--;
1128
1129 if( rlinfo->refcount < 0 )
1130 msRasterLayerInfoFree( layer );
1131 }
1132 return MS_SUCCESS;
1133 }
1134
1135 /************************************************************************/
1136 /* msRASTERLayerNextShape() */
1137 /************************************************************************/
1138
msRASTERLayerNextShape(layerObj * layer,shapeObj * shape)1139 int msRASTERLayerNextShape(layerObj *layer, shapeObj *shape)
1140 {
1141 rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;
1142
1143 if( rlinfo->next_shape < 0
1144 || rlinfo->next_shape >= rlinfo->query_results ) {
1145 msFreeShape(shape);
1146 shape->type = MS_SHAPE_NULL;
1147 return MS_DONE;
1148 } else {
1149 resultObj record;
1150
1151 record.shapeindex = rlinfo->next_shape++;
1152 record.tileindex = 0;
1153 record.classindex = record.resultindex = -1;
1154
1155 return msRASTERLayerGetShape( layer, shape, &record);
1156 }
1157 }
1158
1159 /************************************************************************/
1160 /* msRASTERLayerGetShape() */
1161 /************************************************************************/
1162
msRASTERLayerGetShape(layerObj * layer,shapeObj * shape,resultObj * record)1163 int msRASTERLayerGetShape(layerObj *layer, shapeObj *shape, resultObj *record)
1164 {
1165 rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;
1166 int i;
1167
1168 long shapeindex = record->shapeindex;
1169
1170 msFreeShape(shape);
1171 shape->type = MS_SHAPE_NULL;
1172
1173 /* -------------------------------------------------------------------- */
1174 /* Validate requested record id. */
1175 /* -------------------------------------------------------------------- */
1176 if( shapeindex < 0 || shapeindex >= rlinfo->query_results ) {
1177 msSetError(MS_MISCERR,
1178 "Out of range shape index requested. Requested %ld\n"
1179 "but only %d shapes available.",
1180 "msRASTERLayerGetShape()",
1181 shapeindex, rlinfo->query_results );
1182 return MS_FAILURE;
1183 }
1184
1185 /* -------------------------------------------------------------------- */
1186 /* Apply the geometry. */
1187 /* -------------------------------------------------------------------- */
1188 if( rlinfo->qc_x != NULL ) {
1189 lineObj line;
1190 pointObj point;
1191
1192 shape->type = MS_SHAPE_POINT;
1193
1194 line.numpoints = 1;
1195 line.point = &point;
1196
1197 point.x = rlinfo->qc_x[shapeindex];
1198 point.y = rlinfo->qc_y[shapeindex];
1199 #ifdef USE_POINT_Z_M
1200 point.m = 0.0;
1201 #endif
1202
1203 msAddLine( shape, &line );
1204 msComputeBounds( shape );
1205 }
1206
1207 /* -------------------------------------------------------------------- */
1208 /* Apply the requested items. */
1209 /* -------------------------------------------------------------------- */
1210 if( layer->numitems > 0 ) {
1211 shape->values = (char **) msSmallMalloc(sizeof(char *) * layer->numitems);
1212 shape->numvalues = layer->numitems;
1213
1214 for( i = 0; i < layer->numitems; i++ ) {
1215 const size_t bufferSize = 1000;
1216 char szWork[1000];
1217
1218 szWork[0] = '\0';
1219 if( EQUAL(layer->items[i],"x") && rlinfo->qc_x_reproj )
1220 snprintf( szWork, bufferSize, "%.8g", rlinfo->qc_x_reproj[shapeindex] );
1221 else if( EQUAL(layer->items[i],"y") && rlinfo->qc_y_reproj )
1222 snprintf( szWork, bufferSize, "%.8g", rlinfo->qc_y_reproj[shapeindex] );
1223
1224 else if( EQUAL(layer->items[i],"value_list") && rlinfo->qc_values ) {
1225 int iValue;
1226
1227 for( iValue = 0; iValue < rlinfo->band_count; iValue++ ) {
1228 if( iValue != 0 )
1229 strlcat( szWork, ",", bufferSize);
1230
1231 snprintf( szWork+strlen(szWork), bufferSize-strlen(szWork), "%.8g",
1232 rlinfo->qc_values[shapeindex * rlinfo->band_count
1233 + iValue] );
1234 }
1235 } else if( EQUALN(layer->items[i],"value_",6) && rlinfo->qc_values ) {
1236 int iValue = atoi(layer->items[i]+6);
1237 snprintf( szWork, bufferSize, "%.8g",
1238 rlinfo->qc_values[shapeindex*rlinfo->band_count+iValue] );
1239 } else if( EQUAL(layer->items[i],"class") && rlinfo->qc_class ) {
1240 int p_class = rlinfo->qc_class[shapeindex];
1241 if( layer->class[p_class]->name != NULL )
1242 snprintf( szWork, bufferSize, "%.999s", layer->class[p_class]->name );
1243 else
1244 snprintf( szWork, bufferSize, "%d", p_class );
1245 } else if( EQUAL(layer->items[i],"red") && rlinfo->qc_red )
1246 snprintf( szWork, bufferSize, "%d", rlinfo->qc_red[shapeindex] );
1247 else if( EQUAL(layer->items[i],"green") && rlinfo->qc_green )
1248 snprintf( szWork, bufferSize, "%d", rlinfo->qc_green[shapeindex] );
1249 else if( EQUAL(layer->items[i],"blue") && rlinfo->qc_blue )
1250 snprintf( szWork, bufferSize, "%d", rlinfo->qc_blue[shapeindex] );
1251 else if( EQUAL(layer->items[i],"count") && rlinfo->qc_count )
1252 snprintf( szWork, bufferSize, "%d", rlinfo->qc_count[shapeindex] );
1253
1254 shape->values[i] = msStrdup(szWork);
1255 }
1256 }
1257
1258 /* -------------------------------------------------------------------- */
1259 /* Eventually we should likey apply the geometry properly but */
1260 /* we don't really care about the geometry for query purposes. */
1261 /* -------------------------------------------------------------------- */
1262
1263 return MS_SUCCESS;
1264 }
1265
1266 /************************************************************************/
1267 /* msRASTERLayerGetItems() */
1268 /************************************************************************/
1269
msRASTERLayerGetItems(layerObj * layer)1270 int msRASTERLayerGetItems(layerObj *layer)
1271 {
1272 int maxnumitems = 0;
1273 rasterLayerInfo *rlinfo = (rasterLayerInfo *) layer->layerinfo;
1274
1275 if( rlinfo == NULL )
1276 return MS_FAILURE;
1277
1278 maxnumitems = 8 + (rlinfo->qc_values?rlinfo->band_count:0);
1279 layer->items = (char **) msSmallCalloc(sizeof(char *),maxnumitems);
1280
1281 layer->numitems = 0;
1282 if( rlinfo->qc_x_reproj )
1283 layer->items[layer->numitems++] = msStrdup("x");
1284 if( rlinfo->qc_y_reproj )
1285 layer->items[layer->numitems++] = msStrdup("y");
1286 if( rlinfo->qc_values ) {
1287 int i;
1288 for( i = 0; i < rlinfo->band_count; i++ ) {
1289 char szName[100];
1290 snprintf( szName, sizeof(szName), "value_%d", i );
1291 layer->items[layer->numitems++] = msStrdup(szName);
1292 }
1293 layer->items[layer->numitems++] = msStrdup("value_list");
1294 }
1295 if( rlinfo->qc_class )
1296 layer->items[layer->numitems++] = msStrdup("class");
1297 if( rlinfo->qc_red )
1298 layer->items[layer->numitems++] = msStrdup("red");
1299 if( rlinfo->qc_green )
1300 layer->items[layer->numitems++] = msStrdup("green");
1301 if( rlinfo->qc_blue )
1302 layer->items[layer->numitems++] = msStrdup("blue");
1303 if( rlinfo->qc_count )
1304 layer->items[layer->numitems++] = msStrdup("count");
1305
1306 assert(layer->numitems <= maxnumitems);
1307
1308 return msRASTERLayerInitItemInfo(layer);
1309 }
1310
1311 /************************************************************************/
1312 /* msRASTERLayerGetExtent() */
1313 /************************************************************************/
1314
msRASTERLayerGetExtent(layerObj * layer,rectObj * extent)1315 int msRASTERLayerGetExtent(layerObj *layer, rectObj *extent)
1316
1317 {
1318 char szPath[MS_MAXPATHLEN];
1319 mapObj *map = layer->map;
1320 shapefileObj *tileshpfile;
1321 int tilelayerindex = -1;
1322
1323 if( (!layer->data || strlen(layer->data) == 0)
1324 && layer->tileindex == NULL) {
1325 /* should we be issuing a specific error about not supporting
1326 extents for tileindexed raster layers? */
1327 return MS_FAILURE;
1328 }
1329
1330 if( map == NULL )
1331 return MS_FAILURE;
1332
1333 /* If the layer use a tileindex, return the extent of the tileindex shapefile/referenced layer */
1334 if (layer->tileindex) {
1335 tilelayerindex = msGetLayerIndex(map, layer->tileindex);
1336 if(tilelayerindex != -1) /* does the tileindex reference another layer */
1337 return msLayerGetExtent(GET_LAYER(map, tilelayerindex), extent);
1338 else {
1339 tileshpfile = (shapefileObj *) malloc(sizeof(shapefileObj));
1340 MS_CHECK_ALLOC(tileshpfile, sizeof(shapefileObj), MS_FAILURE);
1341
1342 if(msShapefileOpen(tileshpfile, "rb", msBuildPath3(szPath, map->mappath, map->shapepath, layer->tileindex), MS_TRUE) == -1)
1343 if(msShapefileOpen(tileshpfile, "rb", msBuildPath(szPath, map->mappath, layer->tileindex), MS_TRUE) == -1)
1344 return MS_FAILURE;
1345
1346 *extent = tileshpfile->bounds;
1347 msShapefileClose(tileshpfile);
1348 free(tileshpfile);
1349 return MS_SUCCESS;
1350 }
1351 }
1352
1353 msTryBuildPath3(szPath, map->mappath, map->shapepath, layer->data);
1354 char* decrypted_path = msDecryptStringTokens( map, szPath );
1355 if( !decrypted_path )
1356 return MS_FAILURE;
1357
1358 char** connectionoptions = msGetStringListFromHashTable(&(layer->connectionoptions));
1359 GDALDatasetH hDS = GDALOpenEx(decrypted_path,
1360 GDAL_OF_RASTER,
1361 NULL,
1362 (const char* const*)connectionoptions,
1363 NULL);
1364 CSLDestroy(connectionoptions);
1365 msFree( decrypted_path );
1366 if( hDS == NULL ) {
1367 return MS_FAILURE;
1368 }
1369
1370 const int nXSize = GDALGetRasterXSize( hDS );
1371 const int nYSize = GDALGetRasterYSize( hDS );
1372 double adfGeoTransform[6] = {0};
1373 const CPLErr eErr = GDALGetGeoTransform( hDS, adfGeoTransform );
1374 if( eErr != CE_None ) {
1375 GDALClose( hDS );
1376 return MS_FAILURE;
1377 }
1378
1379 /* If this appears to be an ungeoreferenced raster than flip it for
1380 mapservers purposes. */
1381 if( adfGeoTransform[5] == 1.0 && adfGeoTransform[3] == 0.0 ) {
1382 adfGeoTransform[5] = -1.0;
1383 adfGeoTransform[3] = nYSize;
1384 }
1385
1386 extent->minx = adfGeoTransform[0];
1387 extent->maxy = adfGeoTransform[3];
1388
1389 extent->maxx = adfGeoTransform[0] + nXSize * adfGeoTransform[1];
1390 extent->miny = adfGeoTransform[3] + nYSize * adfGeoTransform[5];
1391
1392 return MS_SUCCESS;
1393 }
1394
1395
1396 /************************************************************************/
1397 /* msRASTERLayerSetTimeFilter() */
1398 /* */
1399 /* This function is actually just used in the context of */
1400 /* setting a filter on the tileindex for time based queries. */
1401 /* For instance via WMS requests. So it isn't really related */
1402 /* to the "raster query" support at all. */
1403 /* */
1404 /* If a local shapefile tileindex is in use, we will set a */
1405 /* backtics filter (shapefile compatible). If another layer is */
1406 /* being used as the tileindex then we will forward the */
1407 /* SetTimeFilter call to it. If there is no tileindex in */
1408 /* place, we do nothing. */
1409 /************************************************************************/
1410
msRASTERLayerSetTimeFilter(layerObj * layer,const char * timestring,const char * timefield)1411 int msRASTERLayerSetTimeFilter(layerObj *layer, const char *timestring,
1412 const char *timefield)
1413 {
1414 int tilelayerindex;
1415
1416 /* -------------------------------------------------------------------- */
1417 /* If we don't have a tileindex the time filter has no effect. */
1418 /* -------------------------------------------------------------------- */
1419 if( layer->tileindex == NULL )
1420 return MS_SUCCESS;
1421
1422 /* -------------------------------------------------------------------- */
1423 /* Find the tileindex layer. */
1424 /* -------------------------------------------------------------------- */
1425 tilelayerindex = msGetLayerIndex(layer->map, layer->tileindex);
1426
1427 /* -------------------------------------------------------------------- */
1428 /* If we are using a local shapefile as our tileindex (that is */
1429 /* to say, the tileindex name is not of another layer), then we */
1430 /* just install a backtics style filter on the raster layer. */
1431 /* This is propogated to the "working layer" created for the */
1432 /* tileindex by code in mapraster.c. */
1433 /* -------------------------------------------------------------------- */
1434 if( tilelayerindex == -1 )
1435 return msLayerMakeBackticsTimeFilter( layer, timestring, timefield );
1436
1437 /* -------------------------------------------------------------------- */
1438 /* Otherwise we invoke the tileindex layers SetTimeFilter */
1439 /* method. */
1440 /* -------------------------------------------------------------------- */
1441 if ( msCheckParentPointer(layer->map,"map")==MS_FAILURE )
1442 return MS_FAILURE;
1443 return msLayerSetTimeFilter( layer->GET_LAYER(map,tilelayerindex),
1444 timestring, timefield );
1445 }
1446
1447 /************************************************************************/
1448 /* msRASTERLayerInitializeVirtualTable() */
1449 /************************************************************************/
1450
1451 int
msRASTERLayerInitializeVirtualTable(layerObj * layer)1452 msRASTERLayerInitializeVirtualTable(layerObj *layer)
1453 {
1454 assert(layer != NULL);
1455 assert(layer->vtable != NULL);
1456
1457 layer->vtable->LayerInitItemInfo = msRASTERLayerInitItemInfo;
1458 layer->vtable->LayerFreeItemInfo = msRASTERLayerFreeItemInfo;
1459 layer->vtable->LayerOpen = msRASTERLayerOpen;
1460 layer->vtable->LayerIsOpen = msRASTERLayerIsOpen;
1461 layer->vtable->LayerWhichShapes = msRASTERLayerWhichShapes;
1462 layer->vtable->LayerNextShape = msRASTERLayerNextShape;
1463 layer->vtable->LayerGetShape = msRASTERLayerGetShape;
1464 /* layer->vtable->LayerGetShapeCount, use default */
1465 layer->vtable->LayerClose = msRASTERLayerClose;
1466 layer->vtable->LayerGetItems = msRASTERLayerGetItems;
1467 layer->vtable->LayerGetExtent = msRASTERLayerGetExtent;
1468 /* layer->vtable->LayerGetAutoStyle, use default */
1469 /* layer->vtable->LayerApplyFilterToLayer, use default */
1470 /* layer->vtable->LayerCloseConnection = msRASTERLayerClose; */
1471 /* we use backtics for proper tileindex shapefile functioning */
1472 layer->vtable->LayerSetTimeFilter = msRASTERLayerSetTimeFilter;
1473 /* layer->vtable->LayerCreateItems, use default */
1474 /* layer->vtable->LayerGetNumFeatures, use default */
1475
1476 return MS_SUCCESS;
1477 }
1478
1479