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