1 /**********************************************************************
2 * $Id: mapuv.c 12629 2011-10-06 18:06:34Z aboudreault $
3 *
4 * Project: MapServer
5 * Purpose: UV Layer
6 * Author: Alan Boudreault (aboudreault@mapgears.com)
7 *
8 **********************************************************************
9 * Copyright (c) 2011, Alan Boudreault, MapGears
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 OR
22 * 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 "mapserver.h"
31
32 #include <assert.h>
33 #include <math.h>
34 #include "mapows.h"
35 #include "mapresample.h"
36 #include "mapthread.h"
37
38 #define MSUVRASTER_NUMITEMS 6
39 #define MSUVRASTER_ANGLE "uv_angle"
40 #define MSUVRASTER_ANGLEINDEX -100
41 #define MSUVRASTER_MINUS_ANGLE "uv_minus_angle"
42 #define MSUVRASTER_MINUSANGLEINDEX -101
43 #define MSUVRASTER_LENGTH "uv_length"
44 #define MSUVRASTER_LENGTHINDEX -102
45 #define MSUVRASTER_LENGTH_2 "uv_length_2"
46 #define MSUVRASTER_LENGTH2INDEX -103
47 #define MSUVRASTER_U "u"
48 #define MSUVRASTER_UINDEX -104
49 #define MSUVRASTER_V "v"
50 #define MSUVRASTER_VINDEX -105
51
52 #define RQM_UNKNOWN 0
53 #define RQM_ENTRY_PER_PIXEL 1
54 #define RQM_HIST_ON_CLASS 2
55 #define RQM_HIST_ON_VALUE 3
56
57 typedef struct {
58
59 /* query cache results */
60 int query_results;
61 /* int query_alloc_max;
62 int query_request_max;
63 int query_result_hard_max;
64 int raster_query_mode; */
65 int band_count;
66
67 int refcount;
68
69 /* query bound in force
70 shapeObj *searchshape;*/
71
72 /* Only nearest result to this point.
73 int range_mode; MS_QUERY_SINGLE, MS_QUERY_MULTIPLE or -1 (skip test)
74 double range_dist;
75 pointObj target_point;*/
76
77 /* double shape_tolerance; */
78
79 float **u; /* u values */
80 float **v; /* v values */
81 int width;
82 int height;
83 rectObj extent;
84 int next_shape;
85 int x, y; /* used internally in msUVRasterLayerNextShape() */
86
87 mapObj* mapToUseForWhichShapes; /* set if the map->extent and map->projection are valid in msUVRASTERLayerWhichShapes() */
88
89 } uvRasterLayerInfo;
90
msUVRASTERLayerUseMapExtentAndProjectionForNextWhichShapes(layerObj * layer,mapObj * map)91 void msUVRASTERLayerUseMapExtentAndProjectionForNextWhichShapes(layerObj* layer,
92 mapObj* map)
93 {
94 uvRasterLayerInfo *uvlinfo = (uvRasterLayerInfo *) layer->layerinfo;
95 uvlinfo->mapToUseForWhichShapes = map;
96 }
97
msUVRASTERLayerInitItemInfo(layerObj * layer)98 static int msUVRASTERLayerInitItemInfo(layerObj *layer)
99 {
100 uvRasterLayerInfo *uvlinfo = (uvRasterLayerInfo *) layer->layerinfo;
101 int i;
102 int *itemindexes;
103 int failed=0;
104
105 if (layer->numitems == 0)
106 return MS_SUCCESS;
107
108 if( uvlinfo == NULL ) {
109 msSetError(MS_MISCERR, "Assertion failed: GDAL layer not opened!!!",
110 "msUVRASTERLayerInitItemInfo()");
111 return(MS_FAILURE);
112 }
113
114 if (layer->iteminfo)
115 free(layer->iteminfo);
116
117 if((layer->iteminfo = (int *)malloc(sizeof(int)*layer->numitems))== NULL) {
118 msSetError(MS_MEMERR, NULL, "msUVRASTERLayerInitItemInfo()");
119 return(MS_FAILURE);
120 }
121
122 itemindexes = (int*)layer->iteminfo;
123 for(i=0; i<layer->numitems; i++) {
124 /* Special case for handling text string and angle coming from */
125 /* OGR style strings. We use special attribute snames. */
126 if (EQUAL(layer->items[i], MSUVRASTER_ANGLE))
127 itemindexes[i] = MSUVRASTER_ANGLEINDEX;
128 else if (EQUAL(layer->items[i], MSUVRASTER_MINUS_ANGLE))
129 itemindexes[i] = MSUVRASTER_MINUSANGLEINDEX;
130 else if (EQUAL(layer->items[i], MSUVRASTER_LENGTH))
131 itemindexes[i] = MSUVRASTER_LENGTHINDEX;
132 else if (EQUAL(layer->items[i], MSUVRASTER_LENGTH_2))
133 itemindexes[i] = MSUVRASTER_LENGTH2INDEX;
134 else if (EQUAL(layer->items[i], MSUVRASTER_U))
135 itemindexes[i] = MSUVRASTER_UINDEX;
136 else if (EQUAL(layer->items[i], MSUVRASTER_V))
137 itemindexes[i] = MSUVRASTER_VINDEX;
138 else {
139 itemindexes[i] = -1;
140 msSetError(MS_OGRERR,
141 "Invalid Field name: %s",
142 "msUVRASTERLayerInitItemInfo()",
143 layer->items[i]);
144 failed=1;
145 }
146 }
147
148 return failed ? (MS_FAILURE) : (MS_SUCCESS);
149 }
150
151
msUVRASTERLayerFreeItemInfo(layerObj * layer)152 void msUVRASTERLayerFreeItemInfo(layerObj *layer)
153 {
154 if (layer->iteminfo)
155 free(layer->iteminfo);
156 layer->iteminfo = NULL;
157 }
158
msUVRasterLayerInfoInitialize(layerObj * layer)159 static void msUVRasterLayerInfoInitialize(layerObj *layer)
160 {
161 uvRasterLayerInfo *uvlinfo = (uvRasterLayerInfo *) layer->layerinfo;
162
163 if( uvlinfo != NULL )
164 return;
165
166 uvlinfo = (uvRasterLayerInfo *) msSmallCalloc(1,sizeof(uvRasterLayerInfo));
167 layer->layerinfo = uvlinfo;
168
169 uvlinfo->band_count = -1;
170 /* uvlinfo->raster_query_mode = RQM_ENTRY_PER_PIXEL; */
171 /* uvlinfo->range_mode = -1; /\* inactive *\/ */
172 /* uvlinfo->refcount = 0; */
173 /* uvlinfo->shape_tolerance = 0.0; */
174 uvlinfo->u = NULL;
175 uvlinfo->v = NULL;
176 uvlinfo->width = 0;
177 uvlinfo->height = 0;
178
179 /* Set attribute type to Real, unless the user has explicitly set */
180 /* something else. */
181 {
182 const char* const items[] = {
183 MSUVRASTER_ANGLE,
184 MSUVRASTER_MINUS_ANGLE,
185 MSUVRASTER_LENGTH,
186 MSUVRASTER_LENGTH_2,
187 MSUVRASTER_U,
188 MSUVRASTER_V,
189 };
190 size_t i;
191 for( i = 0; i < sizeof(items)/sizeof(items[0]); ++i ) {
192 char szTmp[100];
193 snprintf(szTmp, sizeof(szTmp), "%s_type", items[i]);
194 if (msOWSLookupMetadata(&(layer->metadata), "OFG", szTmp) == NULL) {
195 snprintf(szTmp, sizeof(szTmp), "gml_%s_type", items[i]);
196 msInsertHashTable(&(layer->metadata), szTmp, "Real");
197 }
198 }
199 }
200
201 /* uvlinfo->query_result_hard_max = 1000000; */
202
203 /* if( CSLFetchNameValue( layer->processing, "RASTER_QUERY_MAX_RESULT" ) */
204 /* != NULL ) */
205 /* { */
206 /* uvlinfo->query_result_hard_max = */
207 /* atoi(CSLFetchNameValue( layer->processing, "RASTER_QUERY_MAX_RESULT" )); */
208 /* } */
209 }
210
msUVRasterLayerInfoFree(layerObj * layer)211 static void msUVRasterLayerInfoFree( layerObj *layer )
212
213 {
214 uvRasterLayerInfo *uvlinfo = (uvRasterLayerInfo *) layer->layerinfo;
215 int i;
216
217 if( uvlinfo == NULL )
218 return;
219
220 if (uvlinfo->u) {
221 for (i=0; i<uvlinfo->width; ++i) {
222 free(uvlinfo->u[i]);
223 }
224 free(uvlinfo->u);
225 }
226
227 if (uvlinfo->v) {
228 for (i=0; i<uvlinfo->width; ++i) {
229 free(uvlinfo->v[i]);
230 }
231 free(uvlinfo->v);
232 }
233
234 free( uvlinfo );
235
236 layer->layerinfo = NULL;
237 }
238
msUVRASTERLayerOpen(layerObj * layer)239 int msUVRASTERLayerOpen(layerObj *layer)
240 {
241 uvRasterLayerInfo *uvlinfo;
242
243 /* If we don't have info, initialize an empty one now */
244 if( layer->layerinfo == NULL )
245 msUVRasterLayerInfoInitialize( layer );
246
247 uvlinfo = (uvRasterLayerInfo *) layer->layerinfo;
248
249 uvlinfo->refcount = uvlinfo->refcount + 1;
250
251 return MS_SUCCESS;
252 }
253
msUVRASTERLayerIsOpen(layerObj * layer)254 int msUVRASTERLayerIsOpen(layerObj *layer)
255 {
256 if (layer->layerinfo)
257 return MS_TRUE;
258 return MS_FALSE;
259 }
260
261
msUVRASTERLayerClose(layerObj * layer)262 int msUVRASTERLayerClose(layerObj *layer)
263 {
264 uvRasterLayerInfo *uvlinfo = (uvRasterLayerInfo *) layer->layerinfo;
265
266 if( uvlinfo != NULL ) {
267 uvlinfo->refcount--;
268
269 if( uvlinfo->refcount < 1 )
270 msUVRasterLayerInfoFree( layer );
271 }
272 return MS_SUCCESS;
273 }
274
msUVRASTERLayerGetItems(layerObj * layer)275 int msUVRASTERLayerGetItems(layerObj *layer)
276 {
277 uvRasterLayerInfo *uvlinfo = (uvRasterLayerInfo *) layer->layerinfo;
278
279 if( uvlinfo == NULL )
280 return MS_FAILURE;
281
282 layer->numitems = 0;
283 layer->items = (char **) msSmallCalloc(sizeof(char *),10);;
284
285 layer->items[layer->numitems++] = msStrdup(MSUVRASTER_ANGLE);
286 layer->items[layer->numitems++] = msStrdup(MSUVRASTER_MINUS_ANGLE);
287 layer->items[layer->numitems++] = msStrdup(MSUVRASTER_LENGTH);
288 layer->items[layer->numitems++] = msStrdup(MSUVRASTER_LENGTH_2);
289 layer->items[layer->numitems++] = msStrdup(MSUVRASTER_U);
290 layer->items[layer->numitems++] = msStrdup(MSUVRASTER_V);
291 layer->items[layer->numitems] = NULL;
292
293 return msUVRASTERLayerInitItemInfo(layer);
294 }
295
296 /**********************************************************************
297 * msUVRASTERGetValues()
298 *
299 * Special attribute names are used to return some UV params: uv_angle,
300 * uv_length, u and v.
301 **********************************************************************/
msUVRASTERGetValues(layerObj * layer,float * u,float * v)302 static char **msUVRASTERGetValues(layerObj *layer, float *u, float *v)
303 {
304 char **values;
305 int i = 0;
306 char tmp[100];
307 float size_scale;
308 int *itemindexes = (int*)layer->iteminfo;
309
310 if(layer->numitems == 0)
311 return(NULL);
312
313 if(!layer->iteminfo) { /* Should not happen... but just in case! */
314 if (msUVRASTERLayerInitItemInfo(layer) != MS_SUCCESS)
315 return NULL;
316 itemindexes = (int*)layer->iteminfo; /* reassign after malloc */
317 }
318
319 if((values = (char **)malloc(sizeof(char *)*layer->numitems)) == NULL) {
320 msSetError(MS_MEMERR, NULL, "msUVRASTERGetValues()");
321 return(NULL);
322 }
323
324 /* -------------------------------------------------------------------- */
325 /* Determine desired size_scale. Default to 1 if not otherwise set */
326 /* -------------------------------------------------------------------- */
327 size_scale = 1;
328 if( CSLFetchNameValue( layer->processing, "UV_SIZE_SCALE" ) != NULL ) {
329 size_scale =
330 atof(CSLFetchNameValue( layer->processing, "UV_SIZE_SCALE" ));
331 }
332
333 for(i=0; i<layer->numitems; i++) {
334 if (itemindexes[i] == MSUVRASTER_ANGLEINDEX) {
335 snprintf(tmp, 100, "%f", (atan2((double)*v, (double)*u) * 180 / MS_PI));
336 values[i] = msStrdup(tmp);
337 } else if (itemindexes[i] == MSUVRASTER_MINUSANGLEINDEX) {
338 double minus_angle;
339 minus_angle = (atan2((double)*v, (double)*u) * 180 / MS_PI)+180;
340 if (minus_angle >= 360)
341 minus_angle -= 360;
342 snprintf(tmp, 100, "%f", minus_angle);
343 values[i] = msStrdup(tmp);
344 } else if ( (itemindexes[i] == MSUVRASTER_LENGTHINDEX) ||
345 (itemindexes[i] == MSUVRASTER_LENGTH2INDEX) ) {
346 float length = sqrt((*u**u)+(*v**v))*size_scale;
347
348 if (itemindexes[i] == MSUVRASTER_LENGTHINDEX)
349 snprintf(tmp, 100, "%f", length);
350 else
351 snprintf(tmp, 100, "%f", length/2);
352
353 values[i] = msStrdup(tmp);
354 } else if (itemindexes[i] == MSUVRASTER_UINDEX) {
355 snprintf(tmp, 100, "%f",*u);
356 values[i] = msStrdup(tmp);
357 } else if (itemindexes[i] == MSUVRASTER_VINDEX) {
358 snprintf(tmp, 100, "%f",*v);
359 values[i] = msStrdup(tmp);
360 } else {
361 values[i] = NULL;
362 }
363 }
364
365 return values;
366 }
367
msUVRASTERGetSearchRect(layerObj * layer,mapObj * map)368 rectObj msUVRASTERGetSearchRect( layerObj* layer, mapObj* map )
369 {
370 rectObj searchrect = map->extent;
371 int bDone = MS_FALSE;
372
373 /* For UVRaster, it is important that the searchrect is not too large */
374 /* to avoid insufficient intermediate raster resolution, which could */
375 /* happen if we use the default code path, given potential reprojection */
376 /* issues when using a map extent that is not in the validity area of */
377 /* the layer projection. */
378 if( !layer->projection.gt.need_geotransform &&
379 !(msProjIsGeographicCRS(&(map->projection)) &&
380 msProjIsGeographicCRS(&(layer->projection))) ) {
381 rectObj layer_ori_extent;
382
383 if( msLayerGetExtent(layer, &layer_ori_extent) == MS_SUCCESS ) {
384 projectionObj map_proj;
385
386 double map_extent_minx = map->extent.minx;
387 double map_extent_miny = map->extent.miny;
388 double map_extent_maxx = map->extent.maxx;
389 double map_extent_maxy = map->extent.maxy;
390 rectObj layer_extent = layer_ori_extent;
391
392 /* Create a variant of map->projection without geotransform for */
393 /* conveniency */
394 msInitProjection(&map_proj);
395 msCopyProjection(&map_proj, &map->projection);
396 map_proj.gt.need_geotransform = MS_FALSE;
397 if( map->projection.gt.need_geotransform ) {
398 map_extent_minx = map->projection.gt.geotransform[0]
399 + map->projection.gt.geotransform[1] * map->extent.minx
400 + map->projection.gt.geotransform[2] * map->extent.miny;
401 map_extent_miny = map->projection.gt.geotransform[3]
402 + map->projection.gt.geotransform[4] * map->extent.minx
403 + map->projection.gt.geotransform[5] * map->extent.miny;
404 map_extent_maxx = map->projection.gt.geotransform[0]
405 + map->projection.gt.geotransform[1] * map->extent.maxx
406 + map->projection.gt.geotransform[2] * map->extent.maxy;
407 map_extent_maxy = map->projection.gt.geotransform[3]
408 + map->projection.gt.geotransform[4] * map->extent.maxx
409 + map->projection.gt.geotransform[5] * map->extent.maxy;
410 }
411
412 /* Reproject layer extent to map projection */
413 msProjectRect(&layer->projection, &map_proj, &layer_extent);
414
415 if( layer_extent.minx <= map_extent_minx &&
416 layer_extent.miny <= map_extent_miny &&
417 layer_extent.maxx >= map_extent_maxx &&
418 layer_extent.maxy >= map_extent_maxy ) {
419 /* do nothing special if area to map is inside layer extent */
420 }
421 else {
422 if( layer_extent.minx >= map_extent_minx &&
423 layer_extent.maxx <= map_extent_maxx &&
424 layer_extent.miny >= map_extent_miny &&
425 layer_extent.maxy <= map_extent_maxy ) {
426 /* if the area to map is larger than the layer extent, then */
427 /* use full layer extent and add some margin to reflect the */
428 /* proportion of the useful area over the requested bbox */
429 double extra_x =
430 (map_extent_maxx - map_extent_minx) /
431 (layer_extent.maxx - layer_extent.minx) *
432 (layer_ori_extent.maxx - layer_ori_extent.minx);
433 double extra_y =
434 (map_extent_maxy - map_extent_miny) /
435 (layer_extent.maxy - layer_extent.miny) *
436 (layer_ori_extent.maxy - layer_ori_extent.miny);
437 searchrect.minx = layer_ori_extent.minx - extra_x / 2;
438 searchrect.maxx = layer_ori_extent.maxx + extra_x / 2;
439 searchrect.miny = layer_ori_extent.miny - extra_y / 2;
440 searchrect.maxy = layer_ori_extent.maxy + extra_y / 2;
441 }
442 else
443 {
444 /* otherwise clip the map extent with the reprojected layer */
445 /* extent */
446 searchrect.minx = MS_MAX( map_extent_minx, layer_extent.minx );
447 searchrect.maxx = MS_MIN( map_extent_maxx, layer_extent.maxx );
448 searchrect.miny = MS_MAX( map_extent_miny, layer_extent.miny );
449 searchrect.maxy = MS_MIN( map_extent_maxy, layer_extent.maxy );
450 /* and reproject into the layer projection */
451 msProjectRect(&map_proj, &layer->projection, &searchrect);
452 }
453 bDone = MS_TRUE;
454 }
455
456 msFreeProjection(&map_proj);
457 }
458 }
459
460 if( !bDone )
461 msProjectRect(&map->projection, &layer->projection, &searchrect); /* project the searchrect to source coords */
462
463 return searchrect;
464 }
465
msUVRASTERLayerWhichShapes(layerObj * layer,rectObj rect,int isQuery)466 int msUVRASTERLayerWhichShapes(layerObj *layer, rectObj rect, int isQuery)
467 {
468 uvRasterLayerInfo *uvlinfo = (uvRasterLayerInfo *) layer->layerinfo;
469 imageObj *image_tmp;
470 outputFormatObj *outputformat = NULL;
471 mapObj *map_tmp;
472 double map_cellsize;
473 unsigned int spacing;
474 int width, height, u_src_off, v_src_off, i, x, y;
475 char **alteredProcessing = NULL, *saved_layer_mask;
476 char **savedProcessing = NULL;
477 int bHasLonWrap = MS_FALSE;
478 double dfLonWrap = 0.0;
479 rectObj oldLayerExtent = {0};
480 char* oldLayerData = NULL;
481 projectionObj oldLayerProjection={0};
482 int ret;
483
484 if (layer->debug)
485 msDebug("Entering msUVRASTERLayerWhichShapes().\n");
486
487 if( uvlinfo == NULL )
488 return MS_FAILURE;
489
490 if( CSLFetchNameValue( layer->processing, "BANDS" ) == NULL ) {
491 msSetError( MS_MISCERR, "BANDS processing option is required for UV layer. You have to specified 2 bands.",
492 "msUVRASTERLayerWhichShapes()" );
493 return MS_FAILURE;
494 }
495
496 /*
497 ** Allocate mapObj structure
498 */
499 map_tmp = (mapObj *)msSmallCalloc(sizeof(mapObj),1);
500 if(initMap(map_tmp) == -1) { /* initialize this map */
501 msFree(map_tmp);
502 return(MS_FAILURE);
503 }
504
505 /* -------------------------------------------------------------------- */
506 /* Determine desired spacing. Default to 32 if not otherwise set */
507 /* -------------------------------------------------------------------- */
508 spacing = 32;
509 if( CSLFetchNameValue( layer->processing, "UV_SPACING" ) != NULL ) {
510 spacing =
511 atoi(CSLFetchNameValue( layer->processing, "UV_SPACING" ));
512 }
513
514 width = (int)ceil(layer->map->width/spacing);
515 height = (int)ceil(layer->map->height/spacing);
516
517 /* Initialize our dummy map */
518 MS_INIT_COLOR(map_tmp->imagecolor, 255,255,255,255);
519 map_tmp->resolution = layer->map->resolution;
520 map_tmp->defresolution = layer->map->defresolution;
521
522 outputformat = (outputFormatObj *) msSmallCalloc(1,sizeof(outputFormatObj));
523 outputformat->bands = uvlinfo->band_count = 2;
524 outputformat->name = NULL;
525 outputformat->driver = NULL;
526 outputformat->refcount = 0;
527 outputformat->vtable = NULL;
528 outputformat->device = NULL;
529 outputformat->renderer = MS_RENDER_WITH_RAWDATA;
530 outputformat->imagemode = MS_IMAGEMODE_FLOAT32;
531 msAppendOutputFormat(map_tmp, outputformat);
532
533 msCopyHashTable(&map_tmp->configoptions, &layer->map->configoptions);
534 map_tmp->mappath = msStrdup(layer->map->mappath);
535 map_tmp->shapepath = msStrdup(layer->map->shapepath);
536 map_tmp->gt.rotation_angle = 0.0;
537
538 /* Custom msCopyProjection() that removes lon_wrap parameter */
539 {
540 int i;
541
542 map_tmp->projection.numargs = 0;
543 map_tmp->projection.gt = layer->projection.gt;
544 map_tmp->projection.automatic = layer->projection.automatic;
545
546 for (i = 0; i < layer->projection.numargs; i++) {
547 if( strncmp(layer->projection.args[i], "lon_wrap=",
548 strlen("lon_wrap=")) == 0 ) {
549 bHasLonWrap = MS_TRUE;
550 dfLonWrap = atof( layer->projection.args[i] + strlen("lon_wrap=") );
551 }
552 else {
553 map_tmp->projection.args[map_tmp->projection.numargs ++] =
554 msStrdup(layer->projection.args[i]);
555 }
556 }
557 if (map_tmp->projection.numargs != 0) {
558 msProcessProjection(&(map_tmp->projection));
559 }
560
561 map_tmp->projection.wellknownprojection = layer->projection.wellknownprojection;
562 }
563
564 /* Very special case to improve quality for rasters referenced from lon=0 to 360 */
565 /* We create a temporary VRT that swiches the 2 hemispheres, and then we */
566 /* modify the georeferncing to be in the more standard [-180, 180] range */
567 /* and we adjust the layer->data, extent and projection accordingly */
568 if( layer->tileindex == NULL &&
569 uvlinfo->mapToUseForWhichShapes && bHasLonWrap && dfLonWrap == 180.0 )
570 {
571 rectObj layerExtent;
572 msLayerGetExtent(layer, &layerExtent);
573 if( layerExtent.minx == 0 && layerExtent.maxx == 360 )
574 {
575 GDALDatasetH hDS = NULL;
576 char* decrypted_path;
577
578 if( strncmp(layer->data, "<VRTDataset", strlen("<VRTDataset")) == 0 )
579 {
580 decrypted_path = msStrdup(layer->data);
581 }
582 else
583 {
584 char szPath[MS_MAXPATHLEN];
585 msTryBuildPath3(szPath, layer->map->mappath,
586 layer->map->shapepath, layer->data);
587 decrypted_path = msDecryptStringTokens( layer->map, szPath );
588 }
589
590 if( decrypted_path )
591 {
592 char** connectionoptions;
593 GDALAllRegister();
594 connectionoptions = msGetStringListFromHashTable(&(layer->connectionoptions));
595 hDS = GDALOpenEx(decrypted_path,
596 GDAL_OF_RASTER,
597 NULL,
598 (const char* const*)connectionoptions,
599 NULL);
600 CSLDestroy(connectionoptions);
601 }
602 if( hDS != NULL )
603 {
604 int iBand;
605 int nXSize = GDALGetRasterXSize( hDS );
606 int nYSize = GDALGetRasterYSize( hDS );
607 int nBands = GDALGetRasterCount( hDS );
608 int nMaxLen = 100 + nBands * (800 + 2 * strlen(decrypted_path));
609 int nOffset = 0;
610 char* pszInlineVRT = msSmallMalloc( nMaxLen );
611
612 snprintf(pszInlineVRT, nMaxLen,
613 "<VRTDataset rasterXSize=\"%d\" rasterYSize=\"%d\">",
614 nXSize, nYSize);
615 nOffset = strlen(pszInlineVRT);
616 for( iBand = 1; iBand <= nBands; iBand++ )
617 {
618 const char* pszDataType = "Byte";
619 switch( GDALGetRasterDataType(GDALGetRasterBand(hDS, iBand)) )
620 {
621 case GDT_Byte: pszDataType = "Byte"; break;
622 case GDT_Int16: pszDataType = "Int16"; break;
623 case GDT_UInt16: pszDataType = "UInt16"; break;
624 case GDT_Int32: pszDataType = "Int32"; break;
625 case GDT_UInt32: pszDataType = "UInt32"; break;
626 case GDT_Float32: pszDataType = "Float32"; break;
627 case GDT_Float64: pszDataType = "Float64"; break;
628 default: break;
629 }
630
631 snprintf( pszInlineVRT + nOffset, nMaxLen - nOffset,
632 " <VRTRasterBand dataType=\"%s\" band=\"%d\">"
633 " <SimpleSource>"
634 " <SourceFilename relativeToVrt=\"1\"><![CDATA[%s]]></SourceFilename>"
635 " <SourceBand>%d</SourceBand>"
636 " <SrcRect xOff=\"%d\" yOff=\"%d\" xSize=\"%d\" ySize=\"%d\"/>"
637 " <DstRect xOff=\"%d\" yOff=\"%d\" xSize=\"%d\" ySize=\"%d\"/>"
638 " </SimpleSource>"
639 " <SimpleSource>"
640 " <SourceFilename relativeToVrt=\"1\"><![CDATA[%s]]></SourceFilename>"
641 " <SourceBand>%d</SourceBand>"
642 " <SrcRect xOff=\"%d\" yOff=\"%d\" xSize=\"%d\" ySize=\"%d\"/>"
643 " <DstRect xOff=\"%d\" yOff=\"%d\" xSize=\"%d\" ySize=\"%d\"/>"
644 " </SimpleSource>"
645 " </VRTRasterBand>",
646 pszDataType, iBand,
647 decrypted_path, iBand,
648 nXSize / 2, 0, nXSize - nXSize / 2, nYSize,
649 0, 0, nXSize - nXSize / 2, nYSize,
650 decrypted_path, iBand,
651 0, 0, nXSize / 2, nYSize,
652 nXSize - nXSize / 2, 0, nXSize / 2, nYSize );
653
654 nOffset += strlen(pszInlineVRT + nOffset);
655 }
656 snprintf(pszInlineVRT + nOffset, nMaxLen - nOffset,
657 "</VRTDataset>");
658
659 oldLayerExtent = layer->extent;
660 oldLayerData = layer->data;
661 oldLayerProjection = layer->projection;
662 layer->extent.minx = -180;
663 layer->extent.maxx = 180;
664 layer->data = pszInlineVRT;
665 layer->projection = map_tmp->projection;
666
667
668 /* map_tmp->projection is actually layer->projection without lon_wrap */
669 rect = uvlinfo->mapToUseForWhichShapes->extent;
670 msProjectRect(&uvlinfo->mapToUseForWhichShapes->projection,
671 &map_tmp->projection, &rect);
672 bHasLonWrap = MS_FALSE;
673
674 GDALClose(hDS);
675 }
676 msFree( decrypted_path );
677 }
678 }
679
680 if( isQuery )
681 {
682 /* For query mode, use layer->map->extent reprojected rather than */
683 /* the provided rect. Generic query code will filter returned features. */
684 rect = msUVRASTERGetSearchRect(layer, layer->map);
685 }
686
687 map_cellsize = MS_MAX(MS_CELLSIZE(rect.minx, rect.maxx,layer->map->width),
688 MS_CELLSIZE(rect.miny,rect.maxy,layer->map->height));
689 map_tmp->cellsize = map_cellsize*spacing;
690 map_tmp->extent.minx = rect.minx-(0.5*map_cellsize)+(0.5*map_tmp->cellsize);
691 map_tmp->extent.miny = rect.miny-(0.5*map_cellsize)+(0.5*map_tmp->cellsize);
692 map_tmp->extent.maxx = map_tmp->extent.minx+((width-1)*map_tmp->cellsize);
693 map_tmp->extent.maxy = map_tmp->extent.miny+((height-1)*map_tmp->cellsize);
694
695 if( bHasLonWrap && dfLonWrap == 180.0) {
696 if( map_tmp->extent.minx >= 180 ) {
697 /* Request on the right half of the shifted raster (= western hemisphere) */
698 map_tmp->extent.minx -= 360;
699 map_tmp->extent.maxx -= 360;
700 }
701 else if( map_tmp->extent.maxx >= 180.0 ) {
702 /* Request spanning on the 2 hemispheres => drawing whole planet */
703 /* Take only into account vertical resolution, as horizontal one */
704 /* will be unreliable (assuming square pixels...) */
705 map_cellsize = MS_CELLSIZE(rect.miny,rect.maxy,layer->map->height);
706 map_tmp->cellsize = map_cellsize*spacing;
707
708 width = 360.0 / map_tmp->cellsize;
709 map_tmp->extent.minx = -180.0+(0.5*map_tmp->cellsize);
710 map_tmp->extent.maxx = 180.0-(0.5*map_tmp->cellsize);
711 }
712 }
713
714 if (layer->debug)
715 msDebug("msUVRASTERLayerWhichShapes(): width: %d, height: %d, cellsize: %g\n",
716 width, height, map_tmp->cellsize);
717
718 if (layer->debug == 5)
719 msDebug("msUVRASTERLayerWhichShapes(): extent: %g %g %g %g\n",
720 map_tmp->extent.minx, map_tmp->extent.miny,
721 map_tmp->extent.maxx, map_tmp->extent.maxy);
722
723 /* important to use that function, to compute map
724 geotransform, used by the resampling*/
725 msMapSetSize(map_tmp, width, height);
726
727 if (layer->debug == 5)
728 msDebug("msUVRASTERLayerWhichShapes(): geotransform: %g %g %g %g %g %g\n",
729 map_tmp->gt.geotransform[0], map_tmp->gt.geotransform[1],
730 map_tmp->gt.geotransform[2], map_tmp->gt.geotransform[3],
731 map_tmp->gt.geotransform[4], map_tmp->gt.geotransform[5]);
732
733 uvlinfo->extent = map_tmp->extent;
734
735 image_tmp = msImageCreate(width, height, map_tmp->outputformatlist[0],
736 NULL, NULL, map_tmp->resolution, map_tmp->defresolution,
737 &(map_tmp->imagecolor));
738
739 /* Default set to AVERAGE resampling */
740 if( CSLFetchNameValue( layer->processing, "RESAMPLE" ) == NULL ) {
741 alteredProcessing = CSLDuplicate( layer->processing );
742 alteredProcessing =
743 CSLSetNameValue( alteredProcessing, "RESAMPLE",
744 "AVERAGE");
745 savedProcessing = layer->processing;
746 layer->processing = alteredProcessing;
747 }
748
749 /* disable masking at this level: we don't want to apply the mask at the raster level,
750 * it will be applied with the correct cellsize and image size in the vector rendering
751 * phase.
752 */
753 saved_layer_mask = layer->mask;
754 layer->mask = NULL;
755 ret = msDrawRasterLayerLow(map_tmp, layer, image_tmp, NULL );
756
757 /* restore layer attributes if we went through the above on-the-fly VRT */
758 if( oldLayerData )
759 {
760 msFree(layer->data);
761 layer->data = oldLayerData;
762 layer->extent = oldLayerExtent;
763 layer->projection = oldLayerProjection;
764 }
765
766 /* restore layer mask */
767 layer->mask = saved_layer_mask;
768
769 /* restore the saved processing */
770 if (alteredProcessing != NULL) {
771 layer->processing = savedProcessing;
772 CSLDestroy(alteredProcessing);
773 }
774
775 if( ret == MS_FAILURE) {
776 msSetError(MS_MISCERR, "Unable to draw raster data.", "msUVRASTERLayerWhichShapes()");
777
778 msFreeMap(map_tmp);
779 msFreeImage(image_tmp);
780
781 return MS_FAILURE;
782 }
783
784 /* free old query arrays */
785 if (uvlinfo->u) {
786 for (i=0; i<uvlinfo->width; ++i) {
787 free(uvlinfo->u[i]);
788 }
789 free(uvlinfo->u);
790 }
791
792 if (uvlinfo->v) {
793 for (i=0; i<uvlinfo->width; ++i) {
794 free(uvlinfo->v[i]);
795 }
796 free(uvlinfo->v);
797 }
798
799 /* Update our uv layer structure */
800 uvlinfo->width = width;
801 uvlinfo->height = height;
802 uvlinfo->query_results = width*height;
803
804 uvlinfo->u = (float **)msSmallMalloc(sizeof(float *)*width);
805 uvlinfo->v = (float **)msSmallMalloc(sizeof(float *)*width);
806
807 for (x = 0; x < width; ++x) {
808 uvlinfo->u[x] = (float *)msSmallMalloc(height * sizeof(float));
809 uvlinfo->v[x] = (float *)msSmallMalloc(height * sizeof(float));
810
811 for (y = 0; y < height; ++y) {
812 u_src_off = v_src_off = x + y * width;
813 v_src_off += width*height;
814
815 uvlinfo->u[x][y] = image_tmp->img.raw_float[u_src_off];
816 uvlinfo->v[x][y] = image_tmp->img.raw_float[v_src_off];
817
818 /* null vector? update the number of results */
819 if (uvlinfo->u[x][y] == 0 && uvlinfo->v[x][y] == 0)
820 --uvlinfo->query_results;
821 }
822 }
823
824 msFreeImage(image_tmp); /* we do not need the imageObj anymore */
825 msFreeMap(map_tmp);
826
827 uvlinfo->next_shape = 0;
828
829 return MS_SUCCESS;
830 }
831
msUVRASTERLayerGetShape(layerObj * layer,shapeObj * shape,resultObj * record)832 int msUVRASTERLayerGetShape(layerObj *layer, shapeObj *shape, resultObj *record)
833 {
834 uvRasterLayerInfo *uvlinfo = (uvRasterLayerInfo *) layer->layerinfo;
835 lineObj line ;
836 pointObj point;
837 int i, j, k, x=0, y=0;
838 long shapeindex = record->shapeindex;
839
840 msFreeShape(shape);
841 shape->type = MS_SHAPE_NULL;
842
843 if( shapeindex < 0 || shapeindex >= uvlinfo->query_results ) {
844 msSetError(MS_MISCERR,
845 "Out of range shape index requested. Requested %ld\n"
846 "but only %d shapes available.",
847 "msUVRASTERLayerGetShape()",
848 shapeindex, uvlinfo->query_results );
849 return MS_FAILURE;
850 }
851
852 /* loop to the next non null vector */
853 k = 0;
854 for (i=0, x=-1; i<uvlinfo->width && k<=shapeindex; ++i, ++x) {
855 for (j=0, y=-1; j<uvlinfo->height && k<=shapeindex; ++j, ++k, ++y) {
856 if (uvlinfo->u[i][j] == 0 && uvlinfo->v[i][j] == 0)
857 --k;
858 }
859 }
860
861 point.x = Pix2Georef(x, 0, uvlinfo->width-1,
862 uvlinfo->extent.minx, uvlinfo->extent.maxx, MS_FALSE);
863 point.y = Pix2Georef(y, 0, uvlinfo->height-1,
864 uvlinfo->extent.miny, uvlinfo->extent.maxy, MS_TRUE);
865 if (layer->debug == 5)
866 msDebug("msUVRASTERLayerWhichShapes(): shapeindex: %ld, x: %g, y: %g\n",
867 shapeindex, point.x, point.y);
868
869 #ifdef USE_POINT_Z_M
870 point.m = 0.0;
871 #endif
872
873 shape->type = MS_SHAPE_POINT;
874 line.numpoints = 1;
875 line.point = &point;
876 msAddLine( shape, &line );
877 msComputeBounds( shape );
878
879 shape->numvalues = layer->numitems;
880 shape->values = msUVRASTERGetValues(layer, &uvlinfo->u[x][y], &uvlinfo->v[x][y]);
881 shape->index = shapeindex;
882 shape->resultindex = shapeindex;
883
884 return MS_SUCCESS;
885
886 }
887
msUVRASTERLayerNextShape(layerObj * layer,shapeObj * shape)888 int msUVRASTERLayerNextShape(layerObj *layer, shapeObj *shape)
889 {
890 uvRasterLayerInfo *uvlinfo = (uvRasterLayerInfo *) layer->layerinfo;
891
892 if( uvlinfo->next_shape < 0
893 || uvlinfo->next_shape >= uvlinfo->query_results ) {
894 msFreeShape(shape);
895 shape->type = MS_SHAPE_NULL;
896 return MS_DONE;
897 } else {
898 resultObj record;
899
900 record.shapeindex = uvlinfo->next_shape++;
901 record.tileindex = 0;
902 record.classindex = record.resultindex = -1;
903
904 return msUVRASTERLayerGetShape( layer, shape, &record);
905 }
906 }
907
908 /************************************************************************/
909 /* msUVRASTERLayerGetExtent() */
910 /* Simple copy of the maprasterquery.c file. might change in the future */
911 /************************************************************************/
912
msUVRASTERLayerGetExtent(layerObj * layer,rectObj * extent)913 int msUVRASTERLayerGetExtent(layerObj *layer, rectObj *extent)
914
915 {
916 char szPath[MS_MAXPATHLEN];
917 mapObj *map = layer->map;
918 shapefileObj *tileshpfile;
919 int tilelayerindex = -1;
920
921 if( (!layer->data || strlen(layer->data) == 0)
922 && layer->tileindex == NULL) {
923 /* should we be issuing a specific error about not supporting
924 extents for tileindexed raster layers? */
925 return MS_FAILURE;
926 }
927
928 if( map == NULL )
929 return MS_FAILURE;
930
931 /* If the layer use a tileindex, return the extent of the tileindex shapefile/referenced layer */
932 if (layer->tileindex) {
933 tilelayerindex = msGetLayerIndex(map, layer->tileindex);
934 if(tilelayerindex != -1) /* does the tileindex reference another layer */
935 return msLayerGetExtent(GET_LAYER(map, tilelayerindex), extent);
936 else {
937 tileshpfile = (shapefileObj *) malloc(sizeof(shapefileObj));
938 MS_CHECK_ALLOC(tileshpfile, sizeof(shapefileObj), MS_FAILURE);
939
940 if(msShapefileOpen(tileshpfile, "rb", msBuildPath3(szPath, map->mappath, map->shapepath, layer->tileindex), MS_TRUE) == -1)
941 if(msShapefileOpen(tileshpfile, "rb", msBuildPath(szPath, map->mappath, layer->tileindex), MS_TRUE) == -1)
942 return MS_FAILURE;
943
944 *extent = tileshpfile->bounds;
945 msShapefileClose(tileshpfile);
946 free(tileshpfile);
947 return MS_SUCCESS;
948 }
949 }
950
951 msTryBuildPath3(szPath, map->mappath, map->shapepath, layer->data);
952 char* decrypted_path = msDecryptStringTokens( map, szPath );
953 if( !decrypted_path )
954 return MS_FAILURE;
955
956 GDALAllRegister();
957
958 char** connectionoptions = msGetStringListFromHashTable(&(layer->connectionoptions));
959 GDALDatasetH hDS = GDALOpenEx(decrypted_path,
960 GDAL_OF_RASTER,
961 NULL,
962 (const char* const*)connectionoptions,
963 NULL);
964 CSLDestroy(connectionoptions);
965 msFree( decrypted_path );
966 if( hDS == NULL ) {
967 return MS_FAILURE;
968 }
969
970 const int nXSize = GDALGetRasterXSize( hDS );
971 const int nYSize = GDALGetRasterYSize( hDS );
972 double adfGeoTransform[6] = {0};
973 const CPLErr eErr = GDALGetGeoTransform( hDS, adfGeoTransform );
974 if( eErr != CE_None ) {
975 GDALClose( hDS );
976 return MS_FAILURE;
977 }
978
979 /* If this appears to be an ungeoreferenced raster than flip it for
980 mapservers purposes. */
981 if( adfGeoTransform[5] == 1.0 && adfGeoTransform[3] == 0.0 ) {
982 adfGeoTransform[5] = -1.0;
983 adfGeoTransform[3] = nYSize;
984 }
985
986 extent->minx = adfGeoTransform[0];
987 extent->maxy = adfGeoTransform[3];
988
989 extent->maxx = adfGeoTransform[0] + nXSize * adfGeoTransform[1];
990 extent->miny = adfGeoTransform[3] + nYSize * adfGeoTransform[5];
991
992 return MS_SUCCESS;
993 }
994
995
996 /************************************************************************/
997 /* msUVRASTERLayerSetTimeFilter() */
998 /* */
999 /* This function is actually just used in the context of */
1000 /* setting a filter on the tileindex for time based queries. */
1001 /* For instance via WMS requests. So it isn't really related */
1002 /* to the "raster query" support at all. */
1003 /* */
1004 /* If a local shapefile tileindex is in use, we will set a */
1005 /* backtics filter (shapefile compatible). If another layer is */
1006 /* being used as the tileindex then we will forward the */
1007 /* SetTimeFilter call to it. If there is no tileindex in */
1008 /* place, we do nothing. */
1009 /************************************************************************/
1010
msUVRASTERLayerSetTimeFilter(layerObj * layer,const char * timestring,const char * timefield)1011 int msUVRASTERLayerSetTimeFilter(layerObj *layer, const char *timestring,
1012 const char *timefield)
1013 {
1014 int tilelayerindex;
1015
1016 /* -------------------------------------------------------------------- */
1017 /* If we don't have a tileindex the time filter has no effect. */
1018 /* -------------------------------------------------------------------- */
1019 if( layer->tileindex == NULL )
1020 return MS_SUCCESS;
1021
1022 /* -------------------------------------------------------------------- */
1023 /* Find the tileindex layer. */
1024 /* -------------------------------------------------------------------- */
1025 tilelayerindex = msGetLayerIndex(layer->map, layer->tileindex);
1026
1027 /* -------------------------------------------------------------------- */
1028 /* If we are using a local shapefile as our tileindex (that is */
1029 /* to say, the tileindex name is not of another layer), then we */
1030 /* just install a backtics style filter on the raster layer. */
1031 /* This is propogated to the "working layer" created for the */
1032 /* tileindex by code in mapraster.c. */
1033 /* -------------------------------------------------------------------- */
1034 if( tilelayerindex == -1 )
1035 return msLayerMakeBackticsTimeFilter( layer, timestring, timefield );
1036
1037 /* -------------------------------------------------------------------- */
1038 /* Otherwise we invoke the tileindex layers SetTimeFilter */
1039 /* method. */
1040 /* -------------------------------------------------------------------- */
1041 if ( msCheckParentPointer(layer->map,"map")==MS_FAILURE )
1042 return MS_FAILURE;
1043 return msLayerSetTimeFilter( layer->GET_LAYER(map,tilelayerindex),
1044 timestring, timefield );
1045 }
1046
1047
1048 /************************************************************************/
1049 /* msRASTERLayerInitializeVirtualTable() */
1050 /************************************************************************/
1051
1052 int
msUVRASTERLayerInitializeVirtualTable(layerObj * layer)1053 msUVRASTERLayerInitializeVirtualTable(layerObj *layer)
1054 {
1055 assert(layer != NULL);
1056 assert(layer->vtable != NULL);
1057
1058 layer->vtable->LayerInitItemInfo = msUVRASTERLayerInitItemInfo;
1059 layer->vtable->LayerFreeItemInfo = msUVRASTERLayerFreeItemInfo;
1060 layer->vtable->LayerOpen = msUVRASTERLayerOpen;
1061 layer->vtable->LayerIsOpen = msUVRASTERLayerIsOpen;
1062 layer->vtable->LayerWhichShapes = msUVRASTERLayerWhichShapes;
1063 layer->vtable->LayerNextShape = msUVRASTERLayerNextShape;
1064 layer->vtable->LayerGetShape = msUVRASTERLayerGetShape;
1065 /* layer->vtable->LayerGetShapeCount, use default */
1066 layer->vtable->LayerClose = msUVRASTERLayerClose;
1067 layer->vtable->LayerGetItems = msUVRASTERLayerGetItems;
1068 layer->vtable->LayerGetExtent = msUVRASTERLayerGetExtent;
1069 /* layer->vtable->LayerGetAutoStyle, use default */
1070 /* layer->vtable->LayerApplyFilterToLayer, use default */
1071 /* layer->vtable->LayerCloseConnection = msUVRASTERLayerClose; */
1072 /* we use backtics for proper tileindex shapefile functioning */
1073 layer->vtable->LayerSetTimeFilter = msUVRASTERLayerSetTimeFilter;
1074 /* layer->vtable->LayerCreateItems, use default */
1075 /* layer->vtable->LayerGetNumFeatures, use default */
1076
1077 return MS_SUCCESS;
1078 }
1079