1 /******************************************************************************
2 * $id: mapraster.c 10772 2010-11-29 18:27:02Z aboudreault $
3 *
4 * Project: MapServer
5 * Purpose: msDrawRasterLayer(): generic raster layer drawing.
6 * Author: Frank Warmerdam, warmerdam@pobox.com
7 * Pete Olson (LMIC)
8 * Steve Lime
9 *
10 ******************************************************************************
11 * Copyright (c) 1996-2010 Regents of the University of Minnesota.
12 *
13 * Permission is hereby granted, free of charge, to any person obtaining a
14 * copy of this software and associated documentation files (the "Software"),
15 * to deal in the Software without restriction, including without limitation
16 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 * and/or sell copies of the Software, and to permit persons to whom the
18 * Software is furnished to do so, subject to the following conditions:
19 *
20 * The above copyright notice and this permission notice shall be included in
21 * all copies of this Software or works derived from this Software.
22 *
23 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 * DEALINGS IN THE SOFTWARE.
30 ****************************************************************************/
31
32 #include <assert.h>
33 #include "mapserver.h"
34 #include "mapfile.h"
35 #include "mapresample.h"
36 #include "mapthread.h"
37
38
39
40 extern int msyylex_destroy(void);
41 extern int yyparse(parseObj *);
42
43 extern parseResultObj yypresult; /* result of parsing, true/false */
44
45 #include "gdal.h"
46 #include "cpl_string.h"
47 #include "mapraster.h"
48
49 #define MAXCOLORS 256
50 #define BUFLEN 1024
51 #define HDRLEN 8
52
53 #define CVT(x) ((x) >> 8) /* converts to 8-bit color value */
54
55 #define NUMGRAYS 16
56
57 /************************************************************************/
58 /* msGetClass_String() */
59 /************************************************************************/
60
msGetClass_String(layerObj * layer,colorObj * color,const char * pixel_value,int firstClassToTry)61 static int msGetClass_String( layerObj *layer, colorObj *color, const char *pixel_value, int firstClassToTry )
62
63 {
64 int i;
65 const char *tmpstr1=NULL;
66 int numitems;
67 char *item_names[4] = { "pixel", "red", "green", "blue" };
68 char *item_values[4];
69 char red_value[8], green_value[8], blue_value[8];
70
71 /* -------------------------------------------------------------------- */
72 /* No need to do a lookup in the case of one default class. */
73 /* -------------------------------------------------------------------- */
74 if((layer->numclasses == 1) && !(layer->class[0]->expression.string)) /* no need to do lookup */
75 return(0);
76
77 /* -------------------------------------------------------------------- */
78 /* Setup values list for expressions. */
79 /* -------------------------------------------------------------------- */
80 numitems = 4;
81 if( color->red == -1 && color->green == -1 && color->blue == -1 )
82 {
83 strcpy(red_value, "-1");
84 strcpy(green_value, "-1");
85 strcpy(blue_value, "-1");
86 }
87 else
88 {
89 sprintf( red_value, "%d", color->red );
90 sprintf( green_value, "%d", color->green );
91 sprintf( blue_value, "%d", color->blue );
92 }
93
94 item_values[0] = (char *)pixel_value;
95 item_values[1] = red_value;
96 item_values[2] = green_value;
97 item_values[3] = blue_value;
98
99 /* -------------------------------------------------------------------- */
100 /* Loop over classes till we find a match. */
101 /* -------------------------------------------------------------------- */
102 for(i= (firstClassToTry < 0 ) ? 0 : -1; i<layer->numclasses; i++) {
103
104 int idx = i;
105 if( i < 0 )
106 idx = firstClassToTry;
107 else if( i == firstClassToTry )
108 continue;
109
110 /* check for correct classgroup, if set */
111 if ( layer->class[idx]->group && layer->classgroup &&
112 strcasecmp(layer->class[idx]->group, layer->classgroup) != 0 )
113 continue;
114
115 /* Empty expression - always matches */
116 if (layer->class[idx]->expression.string == NULL)
117 return(i);
118
119 switch(layer->class[idx]->expression.type) {
120
121 /* -------------------------------------------------------------------- */
122 /* Simple string match */
123 /* -------------------------------------------------------------------- */
124 case(MS_STRING):
125 /* trim junk white space */
126 tmpstr1= pixel_value;
127 while( *tmpstr1 == ' ' )
128 tmpstr1++;
129
130 if(strcmp(layer->class[idx]->expression.string, tmpstr1) == 0) return(idx); /* matched */
131 break;
132
133 /* -------------------------------------------------------------------- */
134 /* Regular expression. Rarely used for raster. */
135 /* -------------------------------------------------------------------- */
136 case(MS_REGEX):
137 if(!layer->class[idx]->expression.compiled) {
138 if(ms_regcomp(&(layer->class[idx]->expression.regex), layer->class[idx]->expression.string, MS_REG_EXTENDED|MS_REG_NOSUB) != 0) { /* compile the expression */
139 msSetError(MS_REGEXERR, "Invalid regular expression.", "msGetClass()");
140 return(-1);
141 }
142 layer->class[idx]->expression.compiled = MS_TRUE;
143 }
144
145 if(ms_regexec(&(layer->class[idx]->expression.regex), pixel_value, 0, NULL, 0) == 0) return(idx); /* got a match */
146 break;
147
148 /* -------------------------------------------------------------------- */
149 /* Parsed expression. */
150 /* -------------------------------------------------------------------- */
151 case(MS_EXPRESSION): {
152 int status;
153 parseObj p;
154 shapeObj dummy_shape;
155 expressionObj *expression = &(layer->class[idx]->expression);
156
157 dummy_shape.numvalues = numitems;
158 dummy_shape.values = item_values;
159
160 if( expression->tokens == NULL )
161 msTokenizeExpression( expression, item_names, &numitems );
162
163 p.shape = &dummy_shape;
164 p.expr = expression;
165 p.expr->curtoken = p.expr->tokens; /* reset */
166 p.type = MS_PARSE_TYPE_BOOLEAN;
167
168 status = yyparse(&p);
169
170 if (status != 0) {
171 msSetError(MS_PARSEERR, "Failed to parse expression: %s", "msGetClass_FloatRGB", expression->string);
172 return -1;
173 }
174
175 if( p.result.intval )
176 return idx;
177 break;
178 }
179 }
180 }
181
182 return(-1); /* not found */
183 }
184
185 /************************************************************************/
186 /* msGetClass() */
187 /************************************************************************/
188
msGetClass(layerObj * layer,colorObj * color,int colormap_index)189 int msGetClass(layerObj *layer, colorObj *color, int colormap_index)
190 {
191 char pixel_value[12];
192
193 snprintf( pixel_value, sizeof(pixel_value), "%d", colormap_index );
194
195 return msGetClass_String( layer, color, pixel_value, -1 );
196 }
197
198 /************************************************************************/
199 /* msGetClass_FloatRGB() */
200 /* */
201 /* Returns the class based on classification of a floating */
202 /* pixel value. */
203 /************************************************************************/
204
msGetClass_FloatRGB(layerObj * layer,float fValue,int red,int green,int blue)205 int msGetClass_FloatRGB(layerObj *layer, float fValue, int red, int green, int blue )
206 {
207 return msGetClass_FloatRGB_WithFirstClassToTry(
208 layer, fValue, red, green, blue, -1);
209 }
210
211
msGetClass_FloatRGB_WithFirstClassToTry(layerObj * layer,float fValue,int red,int green,int blue,int firstClassToTry)212 int msGetClass_FloatRGB_WithFirstClassToTry(layerObj *layer, float fValue, int red, int green, int blue, int firstClassToTry )
213 {
214 char pixel_value[100];
215 colorObj color;
216
217 color.red = red;
218 color.green = green;
219 color.blue = blue;
220
221 snprintf( pixel_value, sizeof(pixel_value), "%18g", fValue );
222
223 return msGetClass_String( layer, &color, pixel_value, firstClassToTry );
224 }
225
226 /************************************************************************/
227 /* msRasterSetupTileLayer() */
228 /* */
229 /* Setup the tile layer. */
230 /************************************************************************/
231
msDrawRasterSetupTileLayer(mapObj * map,layerObj * layer,rectObj * psearchrect,int is_query,int * ptilelayerindex,int * ptileitemindex,int * ptilesrsindex,layerObj ** ptlp)232 int msDrawRasterSetupTileLayer(mapObj *map, layerObj *layer,
233 rectObj* psearchrect,
234 int is_query,
235 int* ptilelayerindex, /* output */
236 int* ptileitemindex, /* output */
237 int* ptilesrsindex, /* output */
238 layerObj **ptlp /* output */ )
239 {
240 int i;
241 char* requested_fields;
242 int status;
243 layerObj* tlp = NULL;
244
245 *ptilelayerindex = msGetLayerIndex(layer->map, layer->tileindex);
246 if(*ptilelayerindex == -1) { /* the tileindex references a file, not a layer */
247
248 /* so we create a temporary layer */
249 tlp = (layerObj *) malloc(sizeof(layerObj));
250 MS_CHECK_ALLOC(tlp, sizeof(layerObj), MS_FAILURE);
251
252 initLayer(tlp, map);
253 *ptlp = tlp;
254
255 /* set a few parameters for a very basic shapefile-based layer */
256 tlp->name = msStrdup("TILE");
257 tlp->type = MS_LAYER_TILEINDEX;
258 tlp->data = msStrdup(layer->tileindex);
259
260 if( is_query )
261 {
262 tlp->map = map; /*needed when scaletokens are applied, to extract current map scale */
263 for(i = 0; i < layer->numscaletokens; i++) {
264 if(msGrowLayerScaletokens(tlp) == NULL) {
265 return MS_FAILURE;
266 }
267 initScaleToken(&tlp->scaletokens[i]);
268 msCopyScaleToken(&layer->scaletokens[i],&tlp->scaletokens[i]);
269 tlp->numscaletokens++;
270 }
271 }
272
273 if (layer->projection.numargs > 0 &&
274 EQUAL(layer->projection.args[0], "auto"))
275 {
276 tlp->projection.numargs = 1;
277 tlp->projection.args[0] = msStrdup("auto");
278 }
279
280 if (layer->filteritem)
281 tlp->filteritem = msStrdup(layer->filteritem);
282 if (layer->filter.string) {
283 if (layer->filter.type == MS_EXPRESSION) {
284 char* pszTmp =
285 (char *)msSmallMalloc(sizeof(char)*(strlen(layer->filter.string)+3));
286 sprintf(pszTmp,"(%s)",layer->filter.string);
287 msLoadExpressionString(&tlp->filter, pszTmp);
288 free(pszTmp);
289 } else if (layer->filter.type == MS_REGEX ||
290 layer->filter.type == MS_IREGEX) {
291 char* pszTmp =
292 (char *)msSmallMalloc(sizeof(char)*(strlen(layer->filter.string)+3));
293 sprintf(pszTmp,"/%s/",layer->filter.string);
294 msLoadExpressionString(&tlp->filter, pszTmp);
295 free(pszTmp);
296 } else
297 msLoadExpressionString(&tlp->filter, layer->filter.string);
298
299 tlp->filter.type = layer->filter.type;
300 }
301
302 } else {
303 if ( msCheckParentPointer(layer->map,"map")==MS_FAILURE )
304 return MS_FAILURE;
305 tlp = (GET_LAYER(layer->map, *ptilelayerindex));
306 *ptlp = tlp;
307 }
308 status = msLayerOpen(tlp);
309 if(status != MS_SUCCESS) {
310 return status;
311 }
312
313 /* fetch tileitem and tilesrs fields */
314 requested_fields = (char*) msSmallMalloc(sizeof(char)*(strlen(layer->tileitem)+1+
315 (layer->tilesrs ? strlen(layer->tilesrs) : 0) + 1));
316 if( layer->tilesrs != NULL )
317 sprintf(requested_fields, "%s,%s", layer->tileitem, layer->tilesrs);
318 else
319 strcpy(requested_fields, layer->tileitem);
320 status = msLayerWhichItems(tlp, MS_FALSE, requested_fields);
321 msFree(requested_fields);
322 if(status != MS_SUCCESS) {
323 return status;
324 }
325
326 /* get the tileitem and tilesrs index */
327 for(i=0; i<tlp->numitems; i++) {
328 if(strcasecmp(tlp->items[i], layer->tileitem) == 0) {
329 *ptileitemindex = i;
330 }
331 if(layer->tilesrs != NULL &&
332 strcasecmp(tlp->items[i], layer->tilesrs) == 0) {
333 *ptilesrsindex = i;
334 }
335 }
336 if(*ptileitemindex < 0) { /* didn't find it */
337 msSetError(MS_MEMERR,
338 "Could not find attribute %s in tileindex.",
339 "msDrawRasterLayerLow()",
340 layer->tileitem);
341 return MS_FAILURE;
342 }
343 if(layer->tilesrs != NULL && *ptilesrsindex < 0) { /* didn't find it */
344 msSetError(MS_MEMERR,
345 "Could not find attribute %s in tileindex.",
346 "msDrawRasterLayerLow()",
347 layer->tilesrs);
348 return MS_FAILURE;
349 }
350
351 /* if necessary, project the searchrect to source coords */
352 if((map->projection.numargs > 0) && (layer->projection.numargs > 0) &&
353 !EQUAL(layer->projection.args[0], "auto")) {
354 if( msProjectRect(&map->projection, &layer->projection, psearchrect)
355 != MS_SUCCESS ) {
356 msDebug( "msDrawRasterLayerLow(%s): unable to reproject map request rectangle into layer projection, canceling.\n", layer->name );
357 return MS_FAILURE;
358 }
359 }
360 else if((map->projection.numargs > 0) && (tlp->projection.numargs > 0) &&
361 !EQUAL(tlp->projection.args[0], "auto")) {
362 if( msProjectRect(&map->projection, &tlp->projection, psearchrect)
363 != MS_SUCCESS ) {
364 msDebug( "msDrawRasterLayerLow(%s): unable to reproject map request rectangle into layer projection, canceling.\n", layer->name );
365 return MS_FAILURE;
366 }
367 }
368 return msLayerWhichShapes(tlp, *psearchrect, MS_FALSE);
369 }
370
371 /************************************************************************/
372 /* msDrawRasterCleanupTileLayer() */
373 /* */
374 /* Cleanup the tile layer. */
375 /************************************************************************/
376
msDrawRasterCleanupTileLayer(layerObj * tlp,int tilelayerindex)377 void msDrawRasterCleanupTileLayer(layerObj* tlp,
378 int tilelayerindex)
379 {
380 msLayerClose(tlp);
381 if(tilelayerindex == -1) {
382 freeLayer(tlp);
383 free(tlp);
384 }
385 }
386
387 /************************************************************************/
388 /* msDrawRasterIterateTileIndex() */
389 /* */
390 /* Iterate over the tile layer. */
391 /************************************************************************/
392
msDrawRasterIterateTileIndex(layerObj * layer,layerObj * tlp,shapeObj * ptshp,int tileitemindex,int tilesrsindex,char * tilename,size_t sizeof_tilename,char * tilesrsname,size_t sizeof_tilesrsname)393 int msDrawRasterIterateTileIndex(layerObj *layer,
394 layerObj* tlp,
395 shapeObj* ptshp, /* input-output */
396 int tileitemindex,
397 int tilesrsindex,
398 char* tilename, /* output */
399 size_t sizeof_tilename,
400 char* tilesrsname, /* output */
401 size_t sizeof_tilesrsname)
402 {
403 int status;
404
405 status = msLayerNextShape(tlp, ptshp);
406 if( status == MS_FAILURE || status == MS_DONE ) {
407 return status;
408 }
409
410 if(layer->data == NULL || strlen(layer->data) == 0 ) { /* assume whole filename is in attribute field */
411 strlcpy( tilename, ptshp->values[tileitemindex], sizeof_tilename);
412 } else
413 snprintf(tilename, sizeof_tilename, "%s/%s", ptshp->values[tileitemindex], layer->data);
414
415 tilesrsname[0] = '\0';
416
417 if( tilesrsindex >= 0 )
418 {
419 if(ptshp->values[tilesrsindex] != NULL )
420 strlcpy( tilesrsname, ptshp->values[tilesrsindex], sizeof_tilesrsname );
421 }
422
423 msFreeShape(ptshp); /* done with the shape */
424
425 return status;
426 }
427
428 /************************************************************************/
429 /* msDrawRasterBuildRasterPath() */
430 /* */
431 /* Build the path of the raster to open. */
432 /************************************************************************/
433
msDrawRasterBuildRasterPath(mapObj * map,layerObj * layer,const char * filename,char szPath[MS_MAXPATHLEN])434 int msDrawRasterBuildRasterPath(mapObj *map, layerObj *layer,
435 const char* filename,
436 char szPath[MS_MAXPATHLEN] /* output */)
437 {
438 /*
439 ** If using a tileindex then build the path relative to that file if SHAPEPATH is not set.
440 */
441 if(layer->tileindex && !map->shapepath) {
442 char tiAbsFilePath[MS_MAXPATHLEN];
443 char *tiAbsDirPath = NULL;
444
445 msTryBuildPath(tiAbsFilePath, map->mappath, layer->tileindex); /* absolute path to tileindex file */
446 tiAbsDirPath = msGetPath(tiAbsFilePath); /* tileindex file's directory */
447 msBuildPath(szPath, tiAbsDirPath, filename);
448 free(tiAbsDirPath);
449 } else {
450 msTryBuildPath3(szPath, map->mappath, map->shapepath, filename);
451 }
452
453 return MS_SUCCESS;
454 }
455
456 /************************************************************************/
457 /* msDrawRasterGetCPLErrorMsg() */
458 /* */
459 /* Return the CPL error message, and filter out sensitive info. */
460 /************************************************************************/
461
msDrawRasterGetCPLErrorMsg(const char * decrypted_path,const char * szPath)462 const char* msDrawRasterGetCPLErrorMsg(const char* decrypted_path,
463 const char* szPath)
464 {
465 const char *cpl_error_msg = CPLGetLastErrorMsg();
466
467 /* we wish to avoid reporting decrypted paths */
468 if( cpl_error_msg != NULL
469 && strstr(cpl_error_msg,decrypted_path) != NULL
470 && strcmp(decrypted_path,szPath) != 0 )
471 cpl_error_msg = NULL;
472
473 /* we wish to avoid reporting the stock GDALOpen error messages */
474 if( cpl_error_msg != NULL
475 && (strstr(cpl_error_msg,"not recognised as a supported") != NULL
476 || strstr(cpl_error_msg,"does not exist") != NULL) )
477 cpl_error_msg = NULL;
478
479 if( cpl_error_msg == NULL )
480 cpl_error_msg = "";
481
482 return cpl_error_msg;
483 }
484
485 /************************************************************************/
486 /* msDrawRasterLoadProjection() */
487 /* */
488 /* Handle TILESRS or PROJECTION AUTO for each tile. */
489 /************************************************************************/
490
msDrawRasterLoadProjection(layerObj * layer,GDALDatasetH hDS,const char * filename,int tilesrsindex,const char * tilesrsname)491 int msDrawRasterLoadProjection(layerObj *layer,
492 GDALDatasetH hDS,
493 const char* filename,
494 int tilesrsindex,
495 const char* tilesrsname)
496 {
497 /*
498 ** Generate the projection information if using TILESRS.
499 */
500 if( tilesrsindex >= 0 )
501 {
502 const char *pszWKT;
503 if( tilesrsname[0] != '\0' )
504 pszWKT = tilesrsname;
505 else
506 pszWKT = GDALGetProjectionRef( hDS );
507 if( pszWKT != NULL && strlen(pszWKT) > 0 ) {
508 if( msOGCWKT2ProjectionObj(pszWKT, &(layer->projection),
509 layer->debug ) != MS_SUCCESS ) {
510 char szLongMsg[MESSAGELENGTH*2];
511 errorObj *ms_error = msGetErrorObj();
512
513 snprintf( szLongMsg, sizeof(szLongMsg),
514 "%s\n"
515 "PROJECTION '%s' cannot be used for this "
516 "GDAL raster (`%s').",
517 ms_error->message, pszWKT, filename);
518 szLongMsg[MESSAGELENGTH-1] = '\0';
519
520 msSetError(MS_OGRERR, "%s","msDrawRasterLayer()",
521 szLongMsg);
522
523 return MS_FAILURE;
524 }
525 }
526 }
527 /*
528 ** Generate the projection information if using AUTO.
529 */
530 else if (layer->projection.numargs > 0 &&
531 EQUAL(layer->projection.args[0], "auto")) {
532 const char *pszWKT;
533
534 pszWKT = GDALGetProjectionRef( hDS );
535
536 if( pszWKT != NULL && strlen(pszWKT) > 0 ) {
537 if( msOGCWKT2ProjectionObj(pszWKT, &(layer->projection),
538 layer->debug ) != MS_SUCCESS ) {
539 char szLongMsg[MESSAGELENGTH*2];
540 errorObj *ms_error = msGetErrorObj();
541
542 snprintf( szLongMsg, sizeof(szLongMsg),
543 "%s\n"
544 "PROJECTION AUTO cannot be used for this "
545 "GDAL raster (`%s').",
546 ms_error->message, filename);
547 szLongMsg[MESSAGELENGTH-1] = '\0';
548
549 msSetError(MS_OGRERR, "%s","msDrawRasterLayer()",
550 szLongMsg);
551
552 return MS_FAILURE;
553 }
554 }
555 }
556
557 return MS_SUCCESS;
558 }
559
560 typedef enum
561 {
562 CDRT_OK,
563 CDRT_RETURN_MS_FAILURE,
564 CDRT_CONTINUE_NEXT_TILE
565 } CheckDatasetReturnType;
566
567 /************************************************************************/
568 /* msDrawRasterLayerLowCheckDataset() */
569 /************************************************************************/
570
571 static
msDrawRasterLayerLowCheckDataset(mapObj * map,layerObj * layer,GDALDatasetH hDS,const char * decrypted_path,const char * szPath)572 CheckDatasetReturnType msDrawRasterLayerLowCheckDataset(mapObj *map,
573 layerObj *layer,
574 GDALDatasetH hDS,
575 const char* decrypted_path,
576 const char* szPath)
577 {
578 /*
579 ** If GDAL doesn't recognise it, and it wasn't successfully opened
580 ** Generate an error.
581 */
582 if(hDS == NULL) {
583 int ignore_missing = msMapIgnoreMissingData(map);
584 const char *cpl_error_msg = msDrawRasterGetCPLErrorMsg(decrypted_path, szPath);
585
586 if(ignore_missing == MS_MISSING_DATA_FAIL) {
587 msSetError(MS_IOERR, "Corrupt, empty or missing file '%s' for layer '%s'. %s", "msDrawRasterLayerLow()", szPath, layer->name, cpl_error_msg );
588 return(CDRT_RETURN_MS_FAILURE);
589 } else if( ignore_missing == MS_MISSING_DATA_LOG ) {
590 if( layer->debug || layer->map->debug ) {
591 msDebug( "Corrupt, empty or missing file '%s' for layer '%s' ... ignoring this missing data. %s\n", szPath, layer->name, cpl_error_msg );
592 }
593 return(CDRT_CONTINUE_NEXT_TILE);
594 } else if( ignore_missing == MS_MISSING_DATA_IGNORE ) {
595 return(CDRT_CONTINUE_NEXT_TILE);
596 } else {
597 /* never get here */
598 msSetError(MS_IOERR, "msIgnoreMissingData returned unexpected value.", "msDrawRasterLayerLow()");
599 return(CDRT_RETURN_MS_FAILURE);
600 }
601 }
602
603 return CDRT_OK;
604 }
605
606 /************************************************************************/
607 /* msDrawRasterLayerLowOpenDataset() */
608 /************************************************************************/
609
msDrawRasterLayerLowOpenDataset(mapObj * map,layerObj * layer,const char * filename,char szPath[MS_MAXPATHLEN],char ** p_decrypted_path)610 void* msDrawRasterLayerLowOpenDataset(mapObj *map, layerObj *layer,
611 const char* filename,
612 char szPath[MS_MAXPATHLEN],
613 char** p_decrypted_path)
614 {
615 const char* pszPath;
616
617 msGDALInitialize();
618
619 if(layer->debug == MS_TRUE)
620 msDebug( "msDrawRasterLayerLow(%s): Filename is: %s\n", layer->name, filename);
621
622 if( strncmp(filename, "<VRTDataset", strlen("<VRTDataset")) == 0 )
623 {
624 pszPath = filename;
625 }
626 else
627 {
628 msDrawRasterBuildRasterPath(map, layer, filename, szPath);
629 pszPath = szPath;
630 }
631 if(layer->debug == MS_TRUE)
632 msDebug("msDrawRasterLayerLow(%s): Path is: %s\n", layer->name, pszPath);
633
634 /*
635 ** Note: because we do decryption after the above path expansion
636 ** which depends on actually finding a file, it essentially means that
637 ** fancy path manipulation is essentially disabled when using encrypted
638 ** components. But that is mostly ok, since stuff like sde,postgres and
639 ** oracle georaster do not use real paths.
640 */
641 *p_decrypted_path = msDecryptStringTokens( map, pszPath );
642 if( *p_decrypted_path == NULL )
643 return NULL;
644
645 msAcquireLock( TLOCK_GDAL );
646 if( !layer->tileindex )
647 {
648 char** connectionoptions = msGetStringListFromHashTable(&(layer->connectionoptions));
649 GDALDatasetH hDS = GDALOpenEx( *p_decrypted_path,
650 GDAL_OF_RASTER | GDAL_OF_SHARED,
651 NULL,
652 (const char* const*)connectionoptions,
653 NULL);
654 CSLDestroy(connectionoptions);
655 return hDS;
656 }
657 else
658 {
659 return GDALOpenShared( *p_decrypted_path, GA_ReadOnly );
660 }
661 }
662
663 /************************************************************************/
664 /* msDrawRasterLayerLowCloseDataset() */
665 /************************************************************************/
666
msDrawRasterLayerLowCloseDataset(layerObj * layer,void * hDS)667 void msDrawRasterLayerLowCloseDataset(layerObj *layer, void* hDS)
668 {
669 if( hDS )
670 {
671 const char *close_connection;
672 close_connection = msLayerGetProcessingKey( layer,
673 "CLOSE_CONNECTION" );
674
675 if( close_connection == NULL && layer->tileindex == NULL )
676 close_connection = "DEFER";
677
678 {
679 /* Due to how GDAL processes OVERVIEW_LEVEL, datasets returned are */
680 /* not shared, despite being asked to, so close them for real */
681 char** connectionoptions = msGetStringListFromHashTable(&(layer->connectionoptions));
682 if( CSLFetchNameValue(connectionoptions, "OVERVIEW_LEVEL") )
683 close_connection = NULL;
684 CSLDestroy(connectionoptions);
685 }
686
687 if( close_connection != NULL
688 && strcasecmp(close_connection,"DEFER") == 0 ) {
689 GDALDereferenceDataset( (GDALDatasetH)hDS );
690 } else {
691 GDALClose( (GDALDatasetH)hDS );
692 }
693 msReleaseLock( TLOCK_GDAL );
694 }
695 }
696
697
698 /************************************************************************/
699 /* msDrawRasterLayerLowCheckIfMustDraw() */
700 /* */
701 /* Return 1 if the layer should be drawn. */
702 /************************************************************************/
703
msDrawRasterLayerLowCheckIfMustDraw(mapObj * map,layerObj * layer)704 int msDrawRasterLayerLowCheckIfMustDraw(mapObj *map, layerObj *layer)
705 {
706 if(!layer->data && !layer->tileindex && !(layer->connectiontype==MS_KERNELDENSITY)) {
707 if(layer->debug == MS_TRUE)
708 msDebug( "msDrawRasterLayerLow(%s): layer data and tileindex NULL ... doing nothing.", layer->name );
709 return(0);
710 }
711
712 if((layer->status != MS_ON) && (layer->status != MS_DEFAULT)) {
713 if(layer->debug == MS_TRUE)
714 msDebug( "msDrawRasterLayerLow(%s): not status ON or DEFAULT, doing nothing.", layer->name );
715 return(0);
716 }
717
718 if(map->scaledenom > 0) {
719 if((layer->maxscaledenom > 0) && (map->scaledenom > layer->maxscaledenom)) {
720 if(layer->debug == MS_TRUE)
721 msDebug( "msDrawRasterLayerLow(%s): skipping, map scale %.2g > MAXSCALEDENOM=%g\n",
722 layer->name, map->scaledenom, layer->maxscaledenom );
723 return(0);
724 }
725 if((layer->minscaledenom > 0) && (map->scaledenom <= layer->minscaledenom)) {
726 if(layer->debug == MS_TRUE)
727 msDebug( "msDrawRasterLayerLow(%s): skipping, map scale %.2g < MINSCALEDENOM=%g\n",
728 layer->name, map->scaledenom, layer->minscaledenom );
729 return(0);
730 }
731 }
732
733 if(layer->maxscaledenom <= 0 && layer->minscaledenom <= 0) {
734 if((layer->maxgeowidth > 0) && ((map->extent.maxx - map->extent.minx) > layer->maxgeowidth)) {
735 if(layer->debug == MS_TRUE)
736 msDebug( "msDrawRasterLayerLow(%s): skipping, map width %.2g > MAXSCALEDENOM=%g\n", layer->name,
737 (map->extent.maxx - map->extent.minx), layer->maxgeowidth );
738 return(0);
739 }
740 if((layer->mingeowidth > 0) && ((map->extent.maxx - map->extent.minx) < layer->mingeowidth)) {
741 if(layer->debug == MS_TRUE)
742 msDebug( "msDrawRasterLayerLow(%s): skipping, map width %.2g < MINSCALEDENOM=%g\n", layer->name,
743 (map->extent.maxx - map->extent.minx), layer->mingeowidth );
744 return(0);
745 }
746 }
747
748 return 1;
749 }
750
751 /************************************************************************/
752 /* msDrawRasterLayerLow() */
753 /* */
754 /* Check for various file types and act appropriately. Handle */
755 /* tile indexing. */
756 /************************************************************************/
757
msDrawRasterLayerLow(mapObj * map,layerObj * layer,imageObj * image,rasterBufferObj * rb)758 int msDrawRasterLayerLow(mapObj *map, layerObj *layer, imageObj *image,
759 rasterBufferObj *rb )
760 {
761 return msDrawRasterLayerLowWithDataset(map, layer, image, rb, NULL);
762 }
763
msDrawRasterLayerLowWithDataset(mapObj * map,layerObj * layer,imageObj * image,rasterBufferObj * rb,void * hDatasetIn)764 int msDrawRasterLayerLowWithDataset(mapObj *map, layerObj *layer, imageObj *image,
765 rasterBufferObj *rb, void* hDatasetIn )
766 {
767 /* -------------------------------------------------------------------- */
768 /* As of MapServer 6.0 GDAL is required for rendering raster */
769 /* imagery. */
770 /* -------------------------------------------------------------------- */
771 int status, done;
772 char *filename=NULL, tilename[MS_MAXPATHLEN], tilesrsname[1024];
773
774 layerObj *tlp=NULL; /* pointer to the tile layer either real or temporary */
775 int tileitemindex=-1, tilelayerindex=-1, tilesrsindex=-1;
776 shapeObj tshp;
777
778 char szPath[MS_MAXPATHLEN] = { 0 };
779 char *decrypted_path = NULL;
780 int final_status = MS_SUCCESS;
781
782 rectObj searchrect;
783 GDALDatasetH hDS;
784 double adfGeoTransform[6];
785 void *kernel_density_cleanup_ptr = NULL;
786
787 if(layer->debug > 0 || map->debug > 1)
788 msDebug( "msDrawRasterLayerLow(%s): entering.\n", layer->name );
789
790 if( hDatasetIn == NULL &&
791 !msDrawRasterLayerLowCheckIfMustDraw(map, layer) ) {
792 return MS_SUCCESS;
793 }
794
795 msGDALInitialize();
796
797 if(layer->tileindex) { /* we have an index file */
798 msInitShape(&tshp);
799 searchrect = map->extent;
800
801 status = msDrawRasterSetupTileLayer(map, layer,
802 &searchrect,
803 MS_FALSE,
804 &tilelayerindex,
805 &tileitemindex,
806 &tilesrsindex,
807 &tlp);
808 if(status != MS_SUCCESS) {
809 if (status != MS_DONE)
810 final_status = status;
811 goto cleanup;
812 }
813 }
814
815 done = MS_FALSE;
816 while(done != MS_TRUE) {
817
818 if(layer->tileindex) {
819 status = msDrawRasterIterateTileIndex(layer, tlp, &tshp,
820 tileitemindex, tilesrsindex,
821 tilename, sizeof(tilename),
822 tilesrsname, sizeof(tilesrsname));
823 if( status == MS_FAILURE) {
824 final_status = MS_FAILURE;
825 break;
826 }
827
828 if(status == MS_DONE) break; /* no more tiles/images */
829 filename = tilename;
830 } else {
831 filename = layer->data;
832 done = MS_TRUE; /* only one image so we're done after this */
833 }
834
835 if(layer->connectiontype == MS_KERNELDENSITY) {
836 msAcquireLock( TLOCK_GDAL );
837 status = msComputeKernelDensityDataset(map, image, layer, &hDS, &kernel_density_cleanup_ptr);
838 if(status != MS_SUCCESS) {
839 msReleaseLock( TLOCK_GDAL );
840 final_status = status;
841 goto cleanup;
842 }
843 done = MS_TRUE;
844 if(msProjectionsDiffer(&map->projection,&layer->projection)) {
845 char *mapProjStr = msGetProjectionString(&(map->projection));
846
847 /* Set the projection to the map file projection */
848 if (msLoadProjectionString(&(layer->projection), mapProjStr) != 0) {
849 GDALClose( hDS );
850 msReleaseLock( TLOCK_GDAL );
851 msSetError(MS_CGIERR, "Unable to set projection on interpolation layer.", "msDrawRasterLayerLow()");
852 return(MS_FAILURE);
853 }
854 free(mapProjStr);
855 }
856 } else {
857 if(strlen(filename) == 0) continue;
858 if( hDatasetIn )
859 {
860 hDS = (GDALDatasetH)hDatasetIn;
861 }
862 else
863 {
864 hDS = (GDALDatasetH)msDrawRasterLayerLowOpenDataset(
865 map, layer, filename, szPath, &decrypted_path);
866 }
867 }
868
869 if( hDatasetIn == NULL )
870 {
871 CheckDatasetReturnType eRet =
872 msDrawRasterLayerLowCheckDataset(map,layer,hDS,decrypted_path,szPath);
873
874 msFree( decrypted_path );
875 decrypted_path = NULL;
876
877 if( eRet == CDRT_CONTINUE_NEXT_TILE )
878 {
879 msReleaseLock( TLOCK_GDAL );
880 continue;
881 }
882 if( eRet == CDRT_RETURN_MS_FAILURE )
883 {
884 msReleaseLock( TLOCK_GDAL );
885 return MS_FAILURE;
886 }
887 }
888
889 if( msDrawRasterLoadProjection(layer, hDS, filename, tilesrsindex, tilesrsname) != MS_SUCCESS )
890 {
891 if( hDatasetIn == NULL )
892 {
893 GDALClose( hDS );
894 msReleaseLock( TLOCK_GDAL );
895 }
896 final_status = MS_FAILURE;
897 break;
898 }
899
900 msGetGDALGeoTransform( hDS, map, layer, adfGeoTransform );
901
902 /*
903 ** We want to resample if the source image is rotated, if
904 ** the projections differ or if resampling has been explicitly
905 ** requested, or if the image has north-down instead of north-up.
906 */
907
908 if( ((adfGeoTransform[2] != 0.0 || adfGeoTransform[4] != 0.0
909 || adfGeoTransform[5] > 0.0 || adfGeoTransform[1] < 0.0 )
910 && layer->transform )
911 || msProjectionsDiffer( &(map->projection),
912 &(layer->projection) )
913 || CSLFetchNameValue( layer->processing, "RESAMPLE" ) != NULL ) {
914 status = msResampleGDALToMap( map, layer, image, rb, hDS );
915 }
916 else
917 {
918 if( adfGeoTransform[2] != 0.0 || adfGeoTransform[4] != 0.0 ) {
919 if( layer->debug || map->debug )
920 msDebug(
921 "Layer %s has rotational coefficients but we\n"
922 "are unable to use them, projections support\n"
923 "needs to be built in.",
924 layer->name );
925
926 }
927 status = msDrawRasterLayerGDAL(map, layer, image, rb, hDS );
928 }
929
930 if( status == -1 ) {
931 if( hDatasetIn == NULL )
932 {
933 GDALClose( hDS );
934 msReleaseLock( TLOCK_GDAL );
935 }
936 final_status = MS_FAILURE;
937 break;
938 }
939
940 /*
941 ** Should we keep this file open for future use?
942 ** default to keeping open for single data files, and
943 ** to closing for tile indexes
944 */
945 if(layer->connectiontype == MS_KERNELDENSITY) {
946 /*
947 ** Fix issue #5330
948 ** The in-memory kernel density heatmap gdal dataset handle (hDS) gets re-used
949 ** but the associated rasterband cache doesn't get flushed, which causes old data
950 ** to be rendered instead of the newly generated imagery. To fix, simply close the
951 ** the handle and prevent further re-use.
952 ** Note that instead of this workaround, we could explicitely set
953 ** CLOSE_CONNECTION=ALWAYS on the kerneldensity layer.
954 */
955 GDALClose( hDS );
956 msReleaseLock( TLOCK_GDAL );
957 }
958 else {
959 if( hDatasetIn == NULL)
960 msDrawRasterLayerLowCloseDataset(layer, hDS);
961 }
962
963 } /* next tile */
964
965 cleanup:
966 if(layer->tileindex) { /* tiling clean-up */
967 msDrawRasterCleanupTileLayer(tlp, tilelayerindex);
968 }
969 if(layer->connectiontype == MS_KERNELDENSITY && kernel_density_cleanup_ptr) {
970 msCleanupKernelDensityDataset(map, image, layer, kernel_density_cleanup_ptr);
971 }
972
973 return final_status;
974 }
975
976 /************************************************************************/
977 /* msDrawReferenceMap() */
978 /************************************************************************/
979
msDrawReferenceMap(mapObj * map)980 imageObj *msDrawReferenceMap(mapObj *map)
981 {
982 double cellsize;
983 int x1,y1,x2,y2;
984 char szPath[MS_MAXPATHLEN];
985 int status = MS_SUCCESS;
986
987 imageObj *image = NULL;
988 styleObj style;
989
990 /* check to see if we have enough information to actually proceed */
991 if(!map->reference.image || map->reference.height == 0 || map->reference.width == 0) {
992 msSetError(MS_MISCERR, "Reference map configuration error.", "msDrawReferenceMap()");
993 return NULL;
994 }
995
996 rendererVTableObj *renderer = MS_MAP_RENDERER(map);
997 rasterBufferObj *refImage = (rasterBufferObj*)calloc(1,sizeof(rasterBufferObj));
998 MS_CHECK_ALLOC(refImage, sizeof(rasterBufferObj), NULL);
999
1000 if(MS_SUCCESS != renderer->loadImageFromFile(msBuildPath(szPath, map->mappath, map->reference.image),refImage)) {
1001 msSetError(MS_MISCERR,"Error loading reference image %s.","msDrawReferenceMap()",szPath);
1002 free(refImage);
1003 return NULL;
1004 }
1005
1006 image = msImageCreate(refImage->width, refImage->height, map->outputformat,
1007 map->web.imagepath, map->web.imageurl, map->resolution, map->defresolution, &(map->reference.color));
1008 if(!image) return NULL;
1009
1010 status = renderer->mergeRasterBuffer(image,refImage,1.0,0,0,0,0,refImage->width, refImage->height);
1011 msFreeRasterBuffer(refImage);
1012 free(refImage);
1013 if(UNLIKELY(status == MS_FAILURE))
1014 return NULL;
1015
1016 /* make sure the extent given in mapfile fits the image */
1017 cellsize = msAdjustExtent(&(map->reference.extent),
1018 image->width, image->height);
1019
1020 /* convert map extent to reference image coordinates */
1021 x1 = MS_MAP2IMAGE_X(map->extent.minx, map->reference.extent.minx, cellsize);
1022 x2 = MS_MAP2IMAGE_X(map->extent.maxx, map->reference.extent.minx, cellsize);
1023 y1 = MS_MAP2IMAGE_Y(map->extent.maxy, map->reference.extent.maxy, cellsize);
1024 y2 = MS_MAP2IMAGE_Y(map->extent.miny, map->reference.extent.maxy, cellsize);
1025
1026 initStyle(&style);
1027 style.color = map->reference.color;
1028 style.outlinecolor = map->reference.outlinecolor;
1029
1030 /* if extent are smaller than minbox size */
1031 /* draw that extent on the reference image */
1032 if( (abs(x2 - x1) > map->reference.minboxsize) ||
1033 (abs(y2 - y1) > map->reference.minboxsize) ) {
1034 shapeObj rect;
1035 lineObj line;
1036 pointObj points[5];
1037 msInitShape(&rect);
1038
1039 line.point = points;
1040 line.numpoints = 5;
1041 rect.line = &line;
1042 rect.numlines = 1;
1043 rect.type = MS_SHAPE_POLYGON;
1044
1045 line.point[0].x = x1;
1046 line.point[0].y = y1;
1047 line.point[1].x = x1;
1048 line.point[1].y = y2;
1049 line.point[2].x = x2;
1050 line.point[2].y = y2;
1051 line.point[3].x = x2;
1052 line.point[3].y = y1;
1053 line.point[4].x = line.point[0].x;
1054 line.point[4].y = line.point[0].y;
1055
1056 line.numpoints = 5;
1057
1058 if( map->reference.maxboxsize == 0 ||
1059 ((abs(x2 - x1) < map->reference.maxboxsize) &&
1060 (abs(y2 - y1) < map->reference.maxboxsize)) ) {
1061 if(UNLIKELY(MS_FAILURE == msDrawShadeSymbol(map, image, &rect, &style, 1.0))) {
1062 msFreeImage(image);
1063 return NULL;
1064 }
1065 }
1066
1067 } else { /* else draw the marker symbol */
1068 if( map->reference.maxboxsize == 0 ||
1069 ((abs(x2 - x1) < map->reference.maxboxsize) &&
1070 (abs(y2 - y1) < map->reference.maxboxsize)) ) {
1071 style.size = map->reference.markersize;
1072
1073 /* if the marker symbol is specify draw this symbol else draw a cross */
1074 if(map->reference.marker || map->reference.markername) {
1075 pointObj point;
1076 point.x = (double)(x1 + x2)/2;
1077 point.y = (double)(y1 + y2)/2;
1078
1079 if(map->reference.marker) {
1080 style.symbol = map->reference.marker;
1081 } else {
1082 style.symbol = msGetSymbolIndex(&map->symbolset, map->reference.markername, MS_TRUE);
1083 }
1084
1085 if(UNLIKELY(MS_FAILURE == msDrawMarkerSymbol(map, image, &point, &style, 1.0))) {
1086 msFreeImage(image);
1087 return NULL;
1088 }
1089 } else {
1090 int x21, y21;
1091 shapeObj cross;
1092 lineObj lines[4];
1093 pointObj point[8];
1094 int i;
1095 /* determine the center point */
1096 x21 = MS_NINT((x1 + x2)/2);
1097 y21 = MS_NINT((y1 + y2)/2);
1098
1099 msInitShape(&cross);
1100 cross.numlines = 4;
1101 cross.line = lines;
1102 for(i=0; i<4; i++) {
1103 cross.line[i].numpoints = 2;
1104 cross.line[i].point = &(point[2*i]);
1105 }
1106 /* draw a cross */
1107 cross.line[0].point[0].x = x21-8;
1108 cross.line[0].point[0].y = y21;
1109 cross.line[0].point[1].x = x21-3;
1110 cross.line[0].point[1].y = y21;
1111 cross.line[1].point[0].x = x21;
1112 cross.line[1].point[0].y = y21-8;
1113 cross.line[1].point[1].x = x21;
1114 cross.line[1].point[1].y = y21-3;
1115 cross.line[2].point[0].x = x21;
1116 cross.line[2].point[0].y = y21+3;
1117 cross.line[2].point[1].x = x21;
1118 cross.line[2].point[1].y = y21+8;
1119 cross.line[3].point[0].x = x21+3;
1120 cross.line[3].point[0].y = y21;
1121 cross.line[3].point[1].x = x21+8;
1122 cross.line[3].point[1].y = y21;
1123
1124 if(UNLIKELY(MS_FAILURE == msDrawLineSymbol(map,image,&cross,&style,1.0))) {
1125 msFreeImage(image);
1126 return NULL;
1127 }
1128 }
1129 }
1130 }
1131
1132 return(image);
1133 }
1134