PostGIS 3.7.0dev-r@@SVN_REVISION@@
Loading...
Searching...
No Matches
rt_mapalgebra.c
Go to the documentation of this file.
1/*
2 *
3 * WKTRaster - Raster Types for PostGIS
4 * http://trac.osgeo.org/postgis/wiki/WKTRaster
5 *
6 * Copyright (C) 2011-2013 Regents of the University of California
7 * <bkpark@ucdavis.edu>
8 * Copyright (C) 2010-2011 Jorge Arevalo <jorge.arevalo@deimos-space.com>
9 * Copyright (C) 2010-2011 David Zwarg <dzwarg@azavea.com>
10 * Copyright (C) 2009-2011 Pierre Racine <pierre.racine@sbf.ulaval.ca>
11 * Copyright (C) 2009-2011 Mateusz Loskot <mateusz@loskot.net>
12 * Copyright (C) 2008-2009 Sandro Santilli <strk@kbt.io>
13 *
14 * This program is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU General Public License
16 * as published by the Free Software Foundation; either version 2
17 * of the License, or (at your option) any later version.
18 *
19 * This program is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public License
25 * along with this program; if not, write to the Free Software Foundation,
26 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
27 *
28 */
29
30#include "librtcore.h"
31#include "librtcore_internal.h"
32
33/******************************************************************************
34* rt_band_reclass()
35******************************************************************************/
36
37static inline double
39{
40 switch (pixtype) {
41 case PT_1BB:
42 case PT_2BUI:
43 case PT_4BUI:
44 case PT_8BSI:
45 case PT_8BUI:
46 case PT_16BSI:
47 case PT_16BUI:
48 case PT_32BSI:
49 case PT_32BUI:
50 return round(nv);
51 default:
52 return nv;
53 }
54}
55
70 rt_band srcband, rt_pixtype pixtype,
71 uint32_t hasnodata, double nodataval,
72 rt_reclassexpr *exprset, int exprcount
73) {
74 rt_band band = NULL;
75 uint32_t width = 0;
76 uint32_t height = 0;
77 int numval = 0;
78 int memsize = 0;
79 void *mem = NULL;
80 uint32_t src_hasnodata = 0;
81 double src_nodataval = 0.0;
82 int isnodata = 0;
83
84 int rtn;
85 uint32_t x;
86 uint32_t y;
87 int i;
88 double or = 0;
89 double ov = 0;
90 double nr = 0;
91 double nv = 0;
92 int do_nv = 0;
93 rt_reclassexpr expr = NULL;
94
95 assert(NULL != srcband);
96 assert(NULL != exprset && exprcount > 0);
97 RASTER_DEBUGF(4, "exprcount = %d", exprcount);
98 RASTER_DEBUGF(4, "exprset @ %p", exprset);
99
100 /* source nodata */
101 src_hasnodata = rt_band_get_hasnodata_flag(srcband);
102 if (src_hasnodata)
103 rt_band_get_nodata(srcband, &src_nodataval);
104
105 /* size of memory block to allocate */
106 width = rt_band_get_width(srcband);
107 height = rt_band_get_height(srcband);
108 numval = width * height;
109 memsize = rt_pixtype_size(pixtype) * numval;
110 mem = (int *) rtalloc(memsize);
111 if (!mem) {
112 rterror("rt_band_reclass: Could not allocate memory for band");
113 return 0;
114 }
115
116 band = rt_band_new_inline(width, height, pixtype, hasnodata, nodataval, mem);
117 if (!band) {
118 rterror("rt_band_reclass: Could not create new band");
119 rtdealloc(mem);
120 return 0;
121 }
122 rt_band_set_ownsdata_flag(band, 1); /* we DO own this data!!! */
123 rt_band_init_value(band, hasnodata ? nodataval : 0.0);
124
125 RASTER_DEBUGF(3, "rt_band_reclass: new band @ %p", band);
126
127 for (x = 0; x < width; x++) {
128 for (y = 0; y < height; y++) {
129 rtn = rt_band_get_pixel(srcband, x, y, &ov, &isnodata);
130
131 /* error getting value, skip */
132 if (rtn != ES_NONE) {
133 RASTER_DEBUGF(3, "Cannot get value at %d, %d", x, y);
134 continue;
135 }
136 RASTER_DEBUGF(4, "(x, y, ov, isnodata) = (%d, %d, %f, %d)", x, y, ov, isnodata);
137
138 do {
139 do_nv = 0;
140
141 /* no data*/
142 if (hasnodata && isnodata) {
143 do_nv = 1;
144 break;
145 }
146
147 for (i = 0; i < exprcount; i++) {
148 expr = exprset[i];
149
150 /* ov matches min and max*/
151 if (
152 FLT_EQ(expr->src.min, ov) &&
153 FLT_EQ(expr->src.max, ov)
154 ) {
155 do_nv = 1;
156 break;
157 }
158
159 /* process min */
160 if ((
161 expr->src.exc_min && (
162 expr->src.min > ov ||
163 FLT_EQ(expr->src.min, ov)
164 )) || (
165 expr->src.inc_min && (
166 expr->src.min < ov ||
167 FLT_EQ(expr->src.min, ov)
168 )) || (
169 expr->src.min < ov
170 )) {
171 /* process max */
172 if ((
173 expr->src.exc_max && (
174 ov > expr->src.max ||
175 FLT_EQ(expr->src.max, ov)
176 )) || (
177 expr->src.inc_max && (
178 ov < expr->src.max ||
179 FLT_EQ(expr->src.max, ov)
180 )) || (
181 ov < expr->src.max
182 )) {
183 do_nv = 1;
184 break;
185 }
186 }
187 }
188 }
189 while (0);
190
191 /* no expression matched, do not continue */
192 if (!do_nv) continue;
193 RASTER_DEBUGF(3, "Using exprset[%d] unless NODATA", i);
194
195 /* converting a value from one range to another range
196 OldRange = (OldMax - OldMin)
197 NewRange = (NewMax - NewMin)
198 NewValue = (((OldValue - OldMin) * NewRange) / OldRange) + NewMin
199 */
200
201 /* NODATA */
202 if (hasnodata && isnodata) {
203 nv = nodataval;
204 }
205 /*
206 "src" min and max is the same, prevent division by zero
207 set nv to "dst" min, which should be the same as "dst" max
208 */
209 else if (FLT_EQ(expr->src.max, expr->src.min)) {
210 nv = expr->dst.min;
211 }
212 else {
213 or = expr->src.max - expr->src.min;
214 nr = expr->dst.max - expr->dst.min;
215 nv = (((ov - expr->src.min) * nr) / or) + expr->dst.min;
216
217 /* if dst range is from high to low */
218 if (expr->dst.min > expr->dst.max) {
219 if (nv > expr->dst.min)
220 nv = expr->dst.min;
221 else if (nv < expr->dst.max)
222 nv = expr->dst.max;
223 }
224 /* if dst range is from low to high */
225 else {
226 if (nv < expr->dst.min)
227 nv = expr->dst.min;
228 else if (nv > expr->dst.max)
229 nv = expr->dst.max;
230 }
231 }
232
233 /* round the value for integers */
234 nv = rt_band_reclass_round_integer(pixtype, nv);
235
236 RASTER_DEBUGF(4, "(%d, %d) ov: %f or: %f - %f nr: %f - %f nv: %f"
237 , x
238 , y
239 , ov
240 , (NULL != expr) ? expr->src.min : 0
241 , (NULL != expr) ? expr->src.max : 0
242 , (NULL != expr) ? expr->dst.min : 0
243 , (NULL != expr) ? expr->dst.max : 0
244 , nv
245 );
246 if (rt_band_set_pixel(band, x, y, nv, NULL) != ES_NONE) {
247 rterror("rt_band_reclass: Could not assign value to new band");
248 rt_band_destroy(band);
249 rtdealloc(mem);
250 return 0;
251 }
252
253 expr = NULL;
254 }
255 }
256
257 return band;
258}
259
260
261/*
262 * Pre-sorting the classpairs allows us to use bsearch
263 * on the array of classpairs in the pixel reclassification
264 * step later on.
265 */
266static int
267rt_classpair_cmp(const void *aptr, const void *bptr)
268{
269 struct rt_classpair_t *a = (struct rt_classpair_t*)aptr;
270 struct rt_classpair_t *b = (struct rt_classpair_t*)bptr;
271 if (a->src < b->src) return -1;
272 else if (a->src > b->src) return 1;
273 else return 0;
274}
275
276
289 rt_band srcband, rt_reclassmap map,
290 uint32_t hasnodata, double nodataval
291) {
292 rt_band band = NULL;
293 uint32_t width = 0;
294 uint32_t height = 0;
295 int numval = 0;
296 int memsize = 0;
297 void *mem = NULL;
298 uint32_t src_hasnodata = 0;
299 rt_pixtype pixtype;
300
301 assert(NULL != srcband);
302 assert(NULL != map);
303
304 /* Validate that the map has consistent array lengths */
305 if (map->count == 0)
306 rterror("%s: reclassification map must contain at least one class pair", __func__);
307
308 /* Need a nodataval from somewhere */
309 src_hasnodata = rt_band_get_hasnodata_flag(srcband);
310 if (!(src_hasnodata || hasnodata))
311 rterror("%s: source band missing nodata value and nodata value not supplied", __func__);
312
313 /*
314 * Use the nodata value of the input in the case where user
315 * does not supply a nodataval
316 */
317 if (!hasnodata && src_hasnodata) {
318 rt_band_get_nodata(srcband, &nodataval);
319 hasnodata = src_hasnodata;
320 }
321
322 /* size of memory block to allocate */
323 width = rt_band_get_width(srcband);
324 height = rt_band_get_height(srcband);
325 numval = width * height;
326 pixtype = map->dsttype;
327 memsize = rt_pixtype_size(pixtype) * numval;
328 mem = (int *) rtalloc(memsize);
329 if (!mem)
330 rterror("%s: Could not allocate memory for band", __func__);
331
332 band = rt_band_new_inline(width, height, pixtype, hasnodata, nodataval, mem);
333 if (! band) {
334 rterror("%s: Could not add band to raster. Aborting", __func__);
335 rtdealloc(mem);
336 return NULL;
337 }
338 rt_band_set_ownsdata_flag(band, 1); /* we own this data */
339
340 /*
341 * apply nodata as baseline value and then only
342 * map in values that match our map
343 */
344 rt_band_init_value(band, nodataval);
345
346 /*
347 * Put the rt_classpairs in order of source value
348 * so we can bsearch them
349 */
350 qsort(map->pairs, map->count, sizeof(struct rt_classpair_t), rt_classpair_cmp);
351
352 for (uint32_t x = 0; x < width; x++) {
353 for (uint32_t y = 0; y < height; y++) {
354 int isnodata;
355 double nv = nodataval;
356 struct rt_classpair_t query = {0.0, 0.0};
357 struct rt_classpair_t *rslt;
358
359 /* get pixel, skip on error */
360 if (rt_band_get_pixel(srcband, x, y, &query.src, &isnodata) != ES_NONE)
361 continue;
362
363 /* output was already initialized to nodataval */
364 if (isnodata)
365 continue;
366 /*
367 * Look for pixel value in map.
368 * If not in map, skip, leaving the nodata value in place.
369 * Otherwise continue and replace with new value
370 */
371 rslt = bsearch(&query, map->pairs, map->count, sizeof(struct rt_classpair_t), rt_classpair_cmp);
372 if (!rslt)
373 continue;
374 else
375 nv = rslt->dst;
376
377 /* round the new value for integer pixel types */
378 nv = rt_band_reclass_round_integer(pixtype, nv);
379
380 if (rt_band_set_pixel(band, x, y, nv, NULL) != ES_NONE) {
381 rterror("%s: Could not assign value to new band", __func__);
382 rt_band_destroy(band);
383 rtdealloc(mem);
384 return NULL;
385 }
386 }
387 }
388
389 return band;
390}
391
392/******************************************************************************
393* rt_raster_iterator()
394******************************************************************************/
395
398 uint32_t count;
399
402 double **offset;
403 int *width;
404 int *height;
405
406 struct {
410 double *nodataval;
411 double *minval;
413
414 struct {
415 uint16_t x;
416 uint16_t y;
418
419 struct {
420 uint32_t rows;
421 uint32_t columns;
423
424 struct {
425 double **values;
426 int **nodata;
428
430};
431
434 _rti_iterator_arg _param;
435
436 _param = rtalloc(sizeof(struct _rti_iterator_arg_t));
437 if (_param == NULL) {
438 rterror("_rti_iterator_arg_init: Could not allocate memory for _rti_iterator_arg");
439 return NULL;
440 }
441
442 _param->count = 0;
443
444 _param->raster = NULL;
445 _param->isempty = NULL;
446 _param->offset = NULL;
447 _param->width = NULL;
448 _param->height = NULL;
449
450 _param->band.rtband = NULL;
451 _param->band.hasnodata = NULL;
452 _param->band.isnodata = NULL;
453 _param->band.nodataval = NULL;
454 _param->band.minval = NULL;
455
456 _param->distance.x = 0;
457 _param->distance.y = 0;
458
459 _param->dimension.rows = 0;
460 _param->dimension.columns = 0;
461
462 _param->empty.values = NULL;
463 _param->empty.nodata = NULL;
464
465 _param->arg = NULL;
466
467 return _param;
468}
469
470static void
472 uint32_t i = 0;
473
474 if (_param->raster != NULL)
475 rtdealloc(_param->raster);
476 if (_param->isempty != NULL)
477 rtdealloc(_param->isempty);
478 if (_param->width != NULL)
479 rtdealloc(_param->width);
480 if (_param->height != NULL)
481 rtdealloc(_param->height);
482
483 if (_param->band.rtband != NULL)
484 rtdealloc(_param->band.rtband);
485 if (_param->band.hasnodata != NULL)
486 rtdealloc(_param->band.hasnodata);
487 if (_param->band.isnodata != NULL)
488 rtdealloc(_param->band.isnodata);
489 if (_param->band.nodataval != NULL)
490 rtdealloc(_param->band.nodataval);
491 if (_param->band.minval != NULL)
492 rtdealloc(_param->band.minval);
493
494 if (_param->offset != NULL) {
495 for (i = 0; i < _param->count; i++) {
496 if (_param->offset[i] == NULL)
497 continue;
498 rtdealloc(_param->offset[i]);
499 }
500 rtdealloc(_param->offset);
501 }
502
503 if (_param->empty.values != NULL) {
504 for (i = 0; i < _param->dimension.rows; i++) {
505 if (_param->empty.values[i] == NULL)
506 continue;
507 rtdealloc(_param->empty.values[i]);
508 }
509 rtdealloc(_param->empty.values);
510 }
511 if (_param->empty.nodata != NULL) {
512 for (i = 0; i < _param->dimension.rows; i++) {
513 if (_param->empty.nodata[i] == NULL)
514 continue;
515 rtdealloc(_param->empty.nodata[i]);
516 }
517 rtdealloc(_param->empty.nodata);
518 }
519
520 if (_param->arg != NULL) {
521 if (_param->arg->values != NULL)
522 rtdealloc(_param->arg->values);
523 if (_param->arg->nodata != NULL)
524 rtdealloc(_param->arg->nodata);
525 if (_param->arg->src_pixel != NULL) {
526 for (i = 0; i < _param->count; i++) {
527 if (_param->arg->src_pixel[i] == NULL)
528 continue;
529 rtdealloc(_param->arg->src_pixel[i]);
530 }
531
532 rtdealloc(_param->arg->src_pixel);
533 }
534
535 rtdealloc(_param->arg);
536 }
537
538 rtdealloc(_param);
539}
540
541static int
543 _rti_iterator_arg _param,
544 rt_iterator itrset, uint16_t itrcount,
545 uint16_t distancex, uint16_t distancey,
546 int *allnull, int *allempty
547) {
548 int i = 0;
549 int hasband = 0;
550
551 _param->count = itrcount;
552 _param->distance.x = distancex;
553 _param->distance.y = distancey;
554 _param->dimension.columns = distancex * 2 + 1;
555 _param->dimension.rows = distancey * 2 + 1;
556
557 /* allocate memory for children */
558 _param->raster = rtalloc(sizeof(rt_raster) * itrcount);
559 _param->isempty = rtalloc(sizeof(int) * itrcount);
560 _param->width = rtalloc(sizeof(int) * itrcount);
561 _param->height = rtalloc(sizeof(int) * itrcount);
562
563 _param->offset = rtalloc(sizeof(double *) * itrcount);
564
565 _param->band.rtband = rtalloc(sizeof(rt_band) * itrcount);
566 _param->band.hasnodata = rtalloc(sizeof(int) * itrcount);
567 _param->band.isnodata = rtalloc(sizeof(int) * itrcount);
568 _param->band.nodataval = rtalloc(sizeof(double) * itrcount);
569 _param->band.minval = rtalloc(sizeof(double) * itrcount);
570
571 if (
572 _param->raster == NULL ||
573 _param->isempty == NULL ||
574 _param->width == NULL ||
575 _param->height == NULL ||
576 _param->offset == NULL ||
577 _param->band.rtband == NULL ||
578 _param->band.hasnodata == NULL ||
579 _param->band.isnodata == NULL ||
580 _param->band.nodataval == NULL ||
581 _param->band.minval == NULL
582 ) {
583 rterror("_rti_iterator_arg_populate: Could not allocate memory for children of _rti_iterator_arg");
584 return 0;
585 }
586
587 *allnull = 0;
588 *allempty = 0;
589
590 /*
591 check input rasters
592 not empty, band # is valid
593 copy raster pointers and set flags
594 */
595 for (i = 0; i < itrcount; i++) {
596 /* initialize elements */
597 _param->raster[i] = NULL;
598 _param->isempty[i] = 0;
599 _param->width[i] = 0;
600 _param->height[i] = 0;
601
602 _param->offset[i] = NULL;
603
604 _param->band.rtband[i] = NULL;
605 _param->band.hasnodata[i] = 0;
606 _param->band.isnodata[i] = 0;
607 _param->band.nodataval[i] = 0;
608 _param->band.minval[i] = 0;
609
610 /* set isempty */
611 if (itrset[i].raster == NULL) {
612 _param->isempty[i] = 1;
613
614 (*allnull)++;
615 (*allempty)++;
616
617 continue;
618 }
619 else if (rt_raster_is_empty(itrset[i].raster)) {
620 _param->isempty[i] = 1;
621
622 (*allempty)++;
623
624 continue;
625 }
626
627 /* check band number */
628 hasband = rt_raster_has_band(itrset[i].raster, itrset[i].nband);
629 if (!hasband) {
630 if (!itrset[i].nbnodata) {
631 rterror("_rti_iterator_arg_populate: Band %d not found for raster %d", itrset[i].nband, i);
632 return 0;
633 }
634 else {
635 RASTER_DEBUGF(4, "Band %d not found for raster %d. Using NODATA", itrset[i].nband, i);
636 }
637 }
638
639 _param->raster[i] = itrset[i].raster;
640 if (hasband) {
641 _param->band.rtband[i] = rt_raster_get_band(itrset[i].raster, itrset[i].nband);
642 if (_param->band.rtband[i] == NULL) {
643 rterror("_rti_iterator_arg_populate: Could not get band %d for raster %d", itrset[i].nband, i);
644 return 0;
645 }
646
647 /* hasnodata */
648 _param->band.hasnodata[i] = rt_band_get_hasnodata_flag(_param->band.rtband[i]);
649
650 /* hasnodata = TRUE */
651 if (_param->band.hasnodata[i]) {
652 /* nodataval */
653 rt_band_get_nodata(_param->band.rtband[i], &(_param->band.nodataval[i]));
654
655 /* isnodata */
656 _param->band.isnodata[i] = rt_band_get_isnodata_flag(_param->band.rtband[i]);
657 }
658 /* hasnodata = FALSE */
659 else {
660 /* minval */
661 _param->band.minval[i] = rt_band_get_min_value(_param->band.rtband[i]);
662 }
663 }
664
665 /* width, height */
666 _param->width[i] = rt_raster_get_width(_param->raster[i]);
667 _param->height[i] = rt_raster_get_height(_param->raster[i]);
668
669 /* init offset */
670 _param->offset[i] = rtalloc(sizeof(double) * 2);
671 if (_param->offset[i] == NULL) {
672 rterror("_rti_iterator_arg_populate: Could not allocate memory for offsets");
673 return 0;
674 }
675 }
676
677 return 1;
678}
679
680static int
682 uint32_t x = 0;
683 uint32_t y = 0;
684
685 _param->empty.values = rtalloc(sizeof(double *) * _param->dimension.rows);
686 _param->empty.nodata = rtalloc(sizeof(int *) * _param->dimension.rows);
687 if (_param->empty.values == NULL || _param->empty.nodata == NULL) {
688 rterror("_rti_iterator_arg_empty_init: Could not allocate memory for empty values and NODATA");
689 return 0;
690 }
691
692 for (y = 0; y < _param->dimension.rows; y++) {
693 _param->empty.values[y] = rtalloc(sizeof(double) * _param->dimension.columns);
694 _param->empty.nodata[y] = rtalloc(sizeof(int) * _param->dimension.columns);
695
696 if (_param->empty.values[y] == NULL || _param->empty.nodata[y] == NULL) {
697 rterror("_rti_iterator_arg_empty_init: Could not allocate memory for elements of empty values and NODATA");
698 return 0;
699 }
700
701 for (x = 0; x < _param->dimension.columns; x++) {
702 _param->empty.values[y][x] = 0;
703 _param->empty.nodata[y][x] = 1;
704 }
705 }
706
707 return 1;
708}
709
710static int
712 uint32_t i = 0;
713
714 _param->arg = rtalloc(sizeof(struct rt_iterator_arg_t));
715 if (_param->arg == NULL) {
716 rterror("_rti_iterator_arg_callback_init: Could not allocate memory for rt_iterator_arg");
717 return 0;
718 }
719
720 _param->arg->values = NULL;
721 _param->arg->nodata = NULL;
722 _param->arg->src_pixel = NULL;
723
724 /* initialize argument components */
725 _param->arg->values = rtalloc(sizeof(double **) * _param->count);
726 _param->arg->nodata = rtalloc(sizeof(int **) * _param->count);
727 _param->arg->src_pixel = rtalloc(sizeof(int *) * _param->count);
728 if (_param->arg->values == NULL || _param->arg->nodata == NULL || _param->arg->src_pixel == NULL) {
729 rterror("_rti_iterator_arg_callback_init: Could not allocate memory for element of rt_iterator_arg");
730 return 0;
731 }
732 memset(_param->arg->values, 0, sizeof(double **) * _param->count);
733 memset(_param->arg->nodata, 0, sizeof(int **) * _param->count);
734
735 /* initialize pos */
736 for (i = 0; i < _param->count; i++) {
737
738 _param->arg->src_pixel[i] = rtalloc(sizeof(int) * 2);
739 if (_param->arg->src_pixel[i] == NULL) {
740 rterror("_rti_iterator_arg_callback_init: Could not allocate memory for position elements of rt_iterator_arg");
741 return 0;
742 }
743 memset(_param->arg->src_pixel[i], 0, sizeof(int) * 2);
744 }
745
746 _param->arg->rasters = _param->count;
747 _param->arg->rows = _param->dimension.rows;
748 _param->arg->columns = _param->dimension.columns;
749
750 _param->arg->dst_pixel[0] = 0;
751 _param->arg->dst_pixel[1] = 0;
752
753 return 1;
754}
755
756static void
758 uint32_t i = 0;
759 uint32_t y = 0;
760
761 for (i = 0; i < _param->count; i++) {
762 RASTER_DEBUGF(5, "empty at @ %p", _param->empty.values);
763 RASTER_DEBUGF(5, "values at @ %p", _param->arg->values[i]);
764
765 if (_param->arg->values[i] == _param->empty.values) {
766 _param->arg->values[i] = NULL;
767 _param->arg->nodata[i] = NULL;
768
769 continue;
770 }
771
772 for (y = 0; y < _param->dimension.rows; y++) {
773 rtdealloc(_param->arg->values[i][y]);
774 rtdealloc(_param->arg->nodata[i][y]);
775 }
776
777 rtdealloc(_param->arg->values[i]);
778 rtdealloc(_param->arg->nodata[i]);
779
780 _param->arg->values[i] = NULL;
781 _param->arg->nodata[i] = NULL;
782 }
783}
784
824 rt_iterator itrset, uint16_t itrcount,
825 rt_extenttype extenttype, rt_raster customextent,
826 rt_pixtype pixtype,
827 uint8_t hasnodata, double nodataval,
828 uint16_t distancex, uint16_t distancey,
829 rt_mask mask,
830 void *userarg,
831 int (*callback)(
832 rt_iterator_arg arg,
833 void *userarg,
834 double *value,
835 int *nodata
836 ),
837 rt_raster *rtnraster
838) {
839 /* output raster */
840 rt_raster rtnrast = NULL;
841 /* output raster's band */
842 rt_band rtnband = NULL;
843
844 /* working raster */
845 rt_raster rast = NULL;
846
847 _rti_iterator_arg _param = NULL;
848 int allnull = 0;
849 int allempty = 0;
850 int aligned = 0;
851 double offset[4] = {0.};
852 rt_pixel npixels;
853
854 int i = 0;
855 int status = 0;
856 int inextent = 0;
857 int x = 0;
858 int y = 0;
859 int _x = 0;
860 int _y = 0;
861
862 int _width = 0;
863 int _height = 0;
864
865 double minval;
866 double value;
867 int isnodata;
868 int nodata;
869
870 RASTER_DEBUG(3, "Starting...");
871
872 assert(itrset != NULL && itrcount > 0);
873 assert(rtnraster != NULL);
874
875 /* init rtnraster to NULL */
876 *rtnraster = NULL;
877
878 /* check that callback function is not NULL */
879 if (callback == NULL) {
880 rterror("rt_raster_iterator: Callback function not provided");
881 return ES_ERROR;
882 }
883
884 /* check that custom extent is provided if extenttype = ET_CUSTOM */
885 if (extenttype == ET_CUSTOM && rt_raster_is_empty(customextent)) {
886 rterror("rt_raster_iterator: Custom extent cannot be empty if extent type is ET_CUSTOM");
887 return ES_ERROR;
888 }
889
890 /* check that pixtype != PT_END */
891 if (pixtype == PT_END) {
892 rterror("rt_raster_iterator: Pixel type cannot be PT_END");
893 return ES_ERROR;
894 }
895
896 /* initialize _param */
897 if ((_param = _rti_iterator_arg_init()) == NULL) {
898 rterror("rt_raster_iterator: Could not initialize internal variables");
899 return ES_ERROR;
900 }
901
902 /* fill _param */
903 if (!_rti_iterator_arg_populate(_param, itrset, itrcount, distancex, distancey, &allnull, &allempty)) {
904 rterror("rt_raster_iterator: Could not populate for internal variables");
906 return ES_ERROR;
907 }
908
909 /* shortcut if all null, return NULL */
910 if (allnull == itrcount) {
911 RASTER_DEBUG(3, "all rasters are NULL, returning NULL");
912
914
915 return ES_NONE;
916 }
917 /* shortcut if all empty, return empty raster */
918 else if (allempty == itrcount) {
919 RASTER_DEBUG(3, "all rasters are empty, returning empty raster");
920
922
923 rtnrast = rt_raster_new(0, 0);
924 if (rtnrast == NULL) {
925 rterror("rt_raster_iterator: Could not create empty raster");
926 return ES_ERROR;
927 }
928 rt_raster_set_scale(rtnrast, 0, 0);
929
930 *rtnraster = rtnrast;
931 return ES_NONE;
932 }
933
934 /* check that all rasters are aligned */
935 RASTER_DEBUG(3, "checking alignment of all rasters");
936 rast = NULL;
937
938 /* find raster to use as reference */
939 /* use custom if provided */
940 if (extenttype == ET_CUSTOM) {
941 RASTER_DEBUG(4, "using custom extent as reference raster");
942 rast = customextent;
943 }
944 /* use first valid one in _param->raster */
945 else {
946 for (i = 0; i < itrcount; i++) {
947 if (!_param->isempty[i]) {
948 RASTER_DEBUGF(4, "using raster at index %d as reference raster", i);
949 rast = _param->raster[i];
950 break;
951 }
952 }
953 }
954
955 /* no rasters found, SHOULD NEVER BE HERE! */
956 if (rast == NULL) {
957 rterror("rt_raster_iterator: Could not find reference raster to use for alignment tests");
958
960
961 return ES_ERROR;
962 }
963
964 do {
965 aligned = 1;
966
967 /* check custom first if set. also skip if rasters are the same */
968 if (extenttype == ET_CUSTOM && rast != customextent) {
969 if (rt_raster_same_alignment(rast, customextent, &aligned, NULL) != ES_NONE) {
970 rterror("rt_raster_iterator: Could not test for alignment between reference raster and custom extent");
971
973
974 return ES_ERROR;
975 }
976
977 RASTER_DEBUGF(5, "custom extent alignment: %d", aligned);
978 if (!aligned)
979 break;
980 }
981
982 for (i = 0; i < itrcount; i++) {
983 /* skip NULL rasters and if rasters are the same */
984 if (_param->isempty[i] || rast == _param->raster[i])
985 continue;
986
987 if (rt_raster_same_alignment(rast, _param->raster[i], &aligned, NULL) != ES_NONE) {
988 rterror("rt_raster_iterator: Could not test for alignment between reference raster and raster %d", i);
989
991
992 return ES_ERROR;
993 }
994 RASTER_DEBUGF(5, "raster at index %d alignment: %d", i, aligned);
995
996 /* abort checking since a raster isn't aligned */
997 if (!aligned)
998 break;
999 }
1000 }
1001 while (0);
1002
1003 /* not aligned, error */
1004 if (!aligned) {
1005 rterror("rt_raster_iterator: The set of rasters provided (custom extent included, if appropriate) do not have the same alignment");
1006
1008
1009 return ES_ERROR;
1010 }
1011
1012 /* use extenttype to build output raster (no bands though) */
1013 i = -1;
1014 switch (extenttype) {
1015 case ET_INTERSECTION:
1016 case ET_UNION:
1017 /* make copy of first "real" raster */
1018 rtnrast = rtalloc(sizeof(struct rt_raster_t));
1019 if (rtnrast == NULL) {
1020 rterror("rt_raster_iterator: Could not allocate memory for output raster");
1021
1023
1024 return ES_ERROR;
1025 }
1026
1027 for (i = 0; i < itrcount; i++) {
1028 if (!_param->isempty[i]) {
1029 memcpy(rtnrast, _param->raster[i], sizeof(struct rt_raster_t));
1030 break;
1031 }
1032 }
1033 rtnrast->numBands = 0;
1034 rtnrast->bands = NULL;
1035
1036 /* get extent of output raster */
1037 rast = NULL;
1038 for (i = i + 1; i < itrcount; i++) {
1039 if (_param->isempty[i])
1040 continue;
1041
1042 status = rt_raster_from_two_rasters(rtnrast, _param->raster[i], extenttype, &rast, NULL);
1043 rtdealloc(rtnrast);
1044
1045 if (rast == NULL || status != ES_NONE) {
1046 rterror("rt_raster_iterator: Could not compute %s extent of rasters",
1047 extenttype == ET_UNION ? "union" : "intersection"
1048 );
1049
1051
1052 return ES_ERROR;
1053 }
1054 else if (rt_raster_is_empty(rast)) {
1055 rtinfo("rt_raster_iterator: Computed raster for %s extent is empty",
1056 extenttype == ET_UNION ? "union" : "intersection"
1057 );
1058
1060
1061 *rtnraster = rast;
1062 return ES_NONE;
1063 }
1064
1065 rtnrast = rast;
1066 rast = NULL;
1067 }
1068
1069 break;
1070 /*
1071 first, second and last have similar checks
1072 and continue into custom
1073 */
1074 case ET_FIRST:
1075 i = 0;
1076 /* FALLTHROUGH */
1077 case ET_SECOND:
1078 if (i < 0) {
1079 if (itrcount < 2)
1080 i = 0;
1081 else
1082 i = 1;
1083 }
1084 /* FALLTHROUGH */
1085 case ET_LAST:
1086 if (i < 0) i = itrcount - 1;
1087
1088 /* input raster is null, return NULL */
1089 if (_param->raster[i] == NULL) {
1090 RASTER_DEBUGF(3, "returning NULL as %s raster is NULL and extent type is ET_%s",
1091 (i == 0 ? "first" : (i == 1 ? "second" : "last")),
1092 (i == 0 ? "FIRST" : (i == 1 ? "SECOND" : "LAST"))
1093 );
1094
1096
1097 return ES_NONE;
1098 }
1099 /* input raster is empty, return empty raster */
1100 else if (_param->isempty[i]) {
1101 RASTER_DEBUGF(3, "returning empty raster as %s raster is empty and extent type is ET_%s",
1102 (i == 0 ? "first" : (i == 1 ? "second" : "last")),
1103 (i == 0 ? "FIRST" : (i == 1 ? "SECOND" : "LAST"))
1104 );
1105
1107
1108 rtnrast = rt_raster_new(0, 0);
1109 if (rtnrast == NULL) {
1110 rterror("rt_raster_iterator: Could not create empty raster");
1111 return ES_ERROR;
1112 }
1113 rt_raster_set_scale(rtnrast, 0, 0);
1114
1115 *rtnraster = rtnrast;
1116 return ES_NONE;
1117 }
1118 /* FALLTHROUGH */
1119 /* copy the custom extent raster */
1120 case ET_CUSTOM:
1121 rtnrast = rtalloc(sizeof(struct rt_raster_t));
1122 if (rtnrast == NULL) {
1123 rterror("rt_raster_iterator: Could not allocate memory for output raster");
1124
1126
1127 return ES_ERROR;
1128 }
1129
1130 switch (extenttype) {
1131 case ET_CUSTOM:
1132 memcpy(rtnrast, customextent, sizeof(struct rt_raster_t));
1133 break;
1134 /* first, second, last */
1135 default:
1136 memcpy(rtnrast, _param->raster[i], sizeof(struct rt_raster_t));
1137 break;
1138 }
1139 rtnrast->numBands = 0;
1140 rtnrast->bands = NULL;
1141 break;
1142 }
1143
1144 _width = rt_raster_get_width(rtnrast);
1145 _height = rt_raster_get_height(rtnrast);
1146
1147 RASTER_DEBUGF(4, "rtnrast (width, height, ulx, uly, scalex, scaley, skewx, skewy, srid) = (%d, %d, %f, %f, %f, %f, %f, %f, %d)",
1148 _width,
1149 _height,
1150 rt_raster_get_x_offset(rtnrast),
1151 rt_raster_get_y_offset(rtnrast),
1152 rt_raster_get_x_scale(rtnrast),
1153 rt_raster_get_y_scale(rtnrast),
1154 rt_raster_get_x_skew(rtnrast),
1155 rt_raster_get_y_skew(rtnrast),
1156 rt_raster_get_srid(rtnrast)
1157 );
1158
1159 /* init values and NODATA for use with empty rasters */
1160 if (!_rti_iterator_arg_empty_init(_param)) {
1161 rterror("rt_raster_iterator: Could not initialize empty values and NODATA");
1162
1164 rt_raster_destroy(rtnrast);
1165
1166 return ES_ERROR;
1167 }
1168
1169 /* create output band */
1171 rtnrast,
1172 pixtype,
1173 nodataval,
1174 hasnodata, nodataval,
1175 0
1176 ) < 0) {
1177 rterror("rt_raster_iterator: Could not add new band to output raster");
1178
1180 rt_raster_destroy(rtnrast);
1181
1182 return ES_ERROR;
1183 }
1184
1185 /* get output band */
1186 rtnband = rt_raster_get_band(rtnrast, 0);
1187 if (rtnband == NULL) {
1188 rterror("rt_raster_iterator: Could not get new band from output raster");
1189
1191 rt_raster_destroy(rtnrast);
1192
1193 return ES_ERROR;
1194 }
1195
1196 /* output band's minimum value */
1197 minval = rt_band_get_min_value(rtnband);
1198
1199 /* initialize argument for callback function */
1200 if (!_rti_iterator_arg_callback_init(_param)) {
1201 rterror("rt_raster_iterator: Could not initialize callback function argument");
1202
1204 rt_band_destroy(rtnband);
1205 rt_raster_destroy(rtnrast);
1206
1207 return ES_ERROR;
1208 }
1209
1210 /* fill _param->offset */
1211 for (i = 0; i < itrcount; i++) {
1212 if (_param->isempty[i])
1213 continue;
1214
1215 status = rt_raster_from_two_rasters(rtnrast, _param->raster[i], ET_FIRST, &rast, offset);
1216 rtdealloc(rast);
1217 if (status != ES_NONE) {
1218 rterror("rt_raster_iterator: Could not compute raster offsets");
1219
1221 rt_band_destroy(rtnband);
1222 rt_raster_destroy(rtnrast);
1223
1224 return ES_ERROR;
1225 }
1226
1227 _param->offset[i][0] = offset[2];
1228 _param->offset[i][1] = offset[3];
1229 RASTER_DEBUGF(4, "rast %d offset: %f %f", i, offset[2], offset[3]);
1230 }
1231
1232 /* loop over each pixel (POI) of output raster */
1233 /* _x,_y are for output raster */
1234 /* x,y are for input raster */
1235 for (_y = 0; _y < _height; _y++) {
1236 for (_x = 0; _x < _width; _x++) {
1237 RASTER_DEBUGF(4, "iterating output pixel (x, y) = (%d, %d)", _x, _y);
1238 _param->arg->dst_pixel[0] = _x;
1239 _param->arg->dst_pixel[1] = _y;
1240
1241 /* loop through each input raster */
1242 for (i = 0; i < itrcount; i++) {
1243 RASTER_DEBUGF(4, "raster %d", i);
1244
1245 /*
1246 empty raster
1247 OR band does not exist and flag set to use NODATA
1248 OR band is NODATA
1249 */
1250 if (
1251 _param->isempty[i] ||
1252 (_param->band.rtband[i] == NULL && itrset[i].nbnodata) ||
1253 _param->band.isnodata[i]
1254 ) {
1255 RASTER_DEBUG(4, "empty raster, band does not exist or band is NODATA. using empty values and NODATA");
1256
1257 x = _x;
1258 y = _y;
1259
1260 _param->arg->values[i] = _param->empty.values;
1261 _param->arg->nodata[i] = _param->empty.nodata;
1262
1263 continue;
1264 }
1265
1266 /* input raster's X,Y */
1267 x = _x - (int) _param->offset[i][0];
1268 y = _y - (int) _param->offset[i][1];
1269 RASTER_DEBUGF(4, "source pixel (x, y) = (%d, %d)", x, y);
1270
1271 _param->arg->src_pixel[i][0] = x;
1272 _param->arg->src_pixel[i][1] = y;
1273
1274 /* neighborhood */
1275 npixels = NULL;
1276 status = 0;
1277 if (distancex > 0 && distancey > 0) {
1278 RASTER_DEBUG(4, "getting neighborhood");
1279
1281 _param->band.rtband[i],
1282 x, y,
1283 distancex, distancey,
1284 1,
1285 &npixels
1286 );
1287 if (status < 0) {
1288 rterror("rt_raster_iterator: Could not get pixel neighborhood");
1289
1291 rt_band_destroy(rtnband);
1292 rt_raster_destroy(rtnrast);
1293
1294 return ES_ERROR;
1295 }
1296 }
1297
1298 /* get value of POI */
1299 /* get pixel's value */
1300 if (
1301 (x >= 0 && x < _param->width[i]) &&
1302 (y >= 0 && y < _param->height[i])
1303 ) {
1304 RASTER_DEBUG(4, "getting value of POI");
1306 _param->band.rtband[i],
1307 x, y,
1308 &value,
1309 &isnodata
1310 ) != ES_NONE) {
1311 rterror("rt_raster_iterator: Could not get the pixel value of band");
1312
1314 rt_band_destroy(rtnband);
1315 rt_raster_destroy(rtnrast);
1316
1317 return ES_ERROR;
1318 }
1319 inextent = 1;
1320 }
1321 /* outside band extent, set to NODATA */
1322 else {
1323 RASTER_DEBUG(4, "Outside band extent, setting value to NODATA");
1324 /* has NODATA, use NODATA */
1325 if (_param->band.hasnodata[i])
1326 value = _param->band.nodataval[i];
1327 /* no NODATA, use min possible value */
1328 else
1329 value = _param->band.minval[i];
1330
1331 inextent = 0;
1332 isnodata = 1;
1333 }
1334
1335 /* add pixel to neighborhood */
1336 status++;
1337 if (status > 1)
1338 npixels = (rt_pixel) rtrealloc(npixels, sizeof(struct rt_pixel_t) * status);
1339 else
1340 npixels = (rt_pixel) rtalloc(sizeof(struct rt_pixel_t));
1341
1342 if (npixels == NULL) {
1343 rterror("rt_raster_iterator: Could not reallocate memory for neighborhood");
1344
1346 rt_band_destroy(rtnband);
1347 rt_raster_destroy(rtnrast);
1348
1349 return ES_ERROR;
1350 }
1351
1352 npixels[status - 1].x = x;
1353 npixels[status - 1].y = y;
1354 npixels[status - 1].nodata = 1;
1355 npixels[status - 1].value = value;
1356
1357 /* set nodata flag */
1358 if ((!_param->band.hasnodata[i] && inextent) || !isnodata) {
1359 npixels[status - 1].nodata = 0;
1360 }
1361 RASTER_DEBUGF(4, "value, nodata: %f, %d", value, npixels[status - 1].nodata);
1362
1363 /* convert set of rt_pixel to 2D array */
1364 status = rt_pixel_set_to_array(
1365 npixels, status, mask,
1366 x, y,
1367 distancex, distancey,
1368 &(_param->arg->values[i]),
1369 &(_param->arg->nodata[i]),
1370 NULL, NULL
1371 );
1372 rtdealloc(npixels);
1373 if (status != ES_NONE) {
1374 rterror("rt_raster_iterator: Could not create 2D array of neighborhood");
1375
1377 rt_band_destroy(rtnband);
1378 rt_raster_destroy(rtnrast);
1379
1380 return ES_ERROR;
1381 }
1382 }
1383
1384 /* callback */
1385 RASTER_DEBUG(4, "calling callback function");
1386 value = 0;
1387 nodata = 0;
1388 status = callback(_param->arg, userarg, &value, &nodata);
1389
1390 /* free memory from callback */
1392
1393 /* handle callback status */
1394 if (status == 0) {
1395 rterror("rt_raster_iterator: Callback function returned an error");
1396
1398 rt_band_destroy(rtnband);
1399 rt_raster_destroy(rtnrast);
1400
1401 return ES_ERROR;
1402 }
1403
1404 /* burn value to pixel */
1405 status = 0;
1406 if (!nodata) {
1407 status = rt_band_set_pixel(rtnband, _x, _y, value, NULL);
1408 RASTER_DEBUGF(4, "burning pixel (%d, %d) with value: %f", _x, _y, value);
1409 }
1410 else if (!hasnodata) {
1411 status = rt_band_set_pixel(rtnband, _x, _y, minval, NULL);
1412 RASTER_DEBUGF(4, "burning pixel (%d, %d) with minval: %f", _x, _y, minval);
1413 }
1414 else {
1415 RASTER_DEBUGF(4, "NOT burning pixel (%d, %d)", _x, _y);
1416 }
1417 if (status != ES_NONE) {
1418 rterror("rt_raster_iterator: Could not set pixel value");
1419
1421 rt_band_destroy(rtnband);
1422 rt_raster_destroy(rtnrast);
1423
1424 return ES_ERROR;
1425 }
1426 }
1427 }
1428
1429 /* lots of cleanup */
1431
1432 *rtnraster = rtnrast;
1433 return ES_NONE;
1434}
1435
1436/******************************************************************************
1437* rt_raster_colormap()
1438******************************************************************************/
1439
1456
1457static _rti_colormap_arg
1459 _rti_colormap_arg arg = NULL;
1460
1461 arg = rtalloc(sizeof(struct _rti_colormap_arg_t));
1462 if (arg == NULL) {
1463 rterror("_rti_colormap_arg_init: Could not allocate memory for _rti_color_arg");
1464 return NULL;
1465 }
1466
1467 arg->band = NULL;
1468 arg->nodataentry = NULL;
1469 arg->hasnodata = 0;
1470 arg->nodataval = 0;
1471
1472 if (raster == NULL)
1473 arg->raster = NULL;
1474 /* raster provided */
1475 else {
1476 arg->raster = rt_raster_clone(raster, 0);
1477 if (arg->raster == NULL) {
1478 rterror("_rti_colormap_arg_init: Could not create output raster");
1479 return NULL;
1480 }
1481 }
1482
1483 arg->nexpr = 0;
1484 arg->expr = NULL;
1485
1486 arg->npos = 0;
1487 arg->pos = NULL;
1488
1489 return arg;
1490}
1491
1492static void
1494 int i = 0;
1495
1496 if (arg->raster != NULL) {
1497 rt_band band = NULL;
1498
1499 for (i = rt_raster_get_num_bands(arg->raster) - 1; i >= 0; i--) {
1500 band = rt_raster_get_band(arg->raster, i);
1501 if (band != NULL)
1502 rt_band_destroy(band);
1503 }
1504
1506 }
1507
1508 if (arg->nexpr) {
1509 for (i = 0; i < arg->nexpr; i++) {
1510 if (arg->expr[i] != NULL)
1511 rtdealloc(arg->expr[i]);
1512 }
1513 rtdealloc(arg->expr);
1514 }
1515
1516 if (arg->npos)
1517 rtdealloc(arg->pos);
1518
1519 rtdealloc(arg);
1520 arg = NULL;
1521}
1522
1535 rt_raster raster, int nband,
1536 rt_colormap colormap
1537) {
1538 _rti_colormap_arg arg = NULL;
1539 rt_raster rtnraster = NULL;
1540 rt_band band = NULL;
1541 int i = 0;
1542 int j = 0;
1543 int k = 0;
1544
1545 assert(colormap != NULL);
1546
1547 /* empty raster */
1548 if (rt_raster_is_empty(raster))
1549 return NULL;
1550
1551 /* no colormap entries */
1552 if (colormap->nentry < 1) {
1553 rterror("rt_raster_colormap: colormap must have at least one entry");
1554 return NULL;
1555 }
1556
1557 /* nband is valid */
1558 if (!rt_raster_has_band(raster, nband)) {
1559 rterror("rt_raster_colormap: raster has no band at index %d", nband);
1560 return NULL;
1561 }
1562
1563 band = rt_raster_get_band(raster, nband);
1564 if (band == NULL) {
1565 rterror("rt_raster_colormap: Could not get band at index %d", nband);
1566 return NULL;
1567 }
1568
1569 /* init internal variables */
1570 arg = _rti_colormap_arg_init(raster);
1571 if (arg == NULL) {
1572 rterror("rt_raster_colormap: Could not initialize internal variables");
1573 return NULL;
1574 }
1575
1576 /* handle NODATA */
1577 if (rt_band_get_hasnodata_flag(band)) {
1578 arg->hasnodata = 1;
1579 rt_band_get_nodata(band, &(arg->nodataval));
1580 }
1581
1582 /* # of colors */
1583 if (colormap->ncolor < 1) {
1584 rterror("rt_raster_colormap: At least one color must be provided");
1586 return NULL;
1587 }
1588 else if (colormap->ncolor > 4) {
1589 rtinfo("More than four colors indicated. Using only the first four colors");
1590 colormap->ncolor = 4;
1591 }
1592
1593 /* find non-NODATA entries */
1594 arg->npos = 0;
1595 arg->pos = rtalloc(sizeof(uint16_t) * colormap->nentry);
1596 if (arg->pos == NULL) {
1597 rterror("rt_raster_colormap: Could not allocate memory for valid entries");
1599 return NULL;
1600 }
1601 for (i = 0, j = 0; i < colormap->nentry; i++) {
1602 /* special handling for NODATA entries */
1603 if (colormap->entry[i].isnodata) {
1604 /* first NODATA entry found, use it */
1605 if (arg->nodataentry == NULL)
1606 arg->nodataentry = &(colormap->entry[i]);
1607 else
1608 rtwarn("More than one colormap entry found for NODATA value. Only using first NOTDATA entry");
1609
1610 continue;
1611 }
1612
1613 (arg->npos)++;
1614 arg->pos[j++] = i;
1615 }
1616
1617 /* INTERPOLATE and only one non-NODATA entry */
1618 if (colormap->method == CM_INTERPOLATE && arg->npos < 2) {
1619 rtwarn("Method INTERPOLATE requires at least two non-NODATA colormap entries. Using NEAREST instead");
1620 colormap->method = CM_NEAREST;
1621 }
1622
1623 /* NODATA entry but band has no NODATA value */
1624 if (!arg->hasnodata && arg->nodataentry != NULL) {
1625 rtinfo("Band at index %d has no NODATA value. Ignoring NODATA entry", nband);
1626 arg->nodataentry = NULL;
1627 }
1628
1629 /* allocate expr */
1630 arg->nexpr = arg->npos;
1631
1632 /* INTERPOLATE needs one less than the number of entries */
1633 if (colormap->method == CM_INTERPOLATE)
1634 arg->nexpr -= 1;
1635 /* EXACT requires a no matching expression */
1636 else if (colormap->method == CM_EXACT)
1637 arg->nexpr += 1;
1638
1639 /* NODATA entry exists, add expression */
1640 if (arg->nodataentry != NULL)
1641 arg->nexpr += 1;
1642 arg->expr = rtalloc(sizeof(rt_reclassexpr) * arg->nexpr);
1643 if (arg->expr == NULL) {
1644 rterror("rt_raster_colormap: Could not allocate memory for reclass expressions");
1646 return NULL;
1647 }
1648 RASTER_DEBUGF(4, "nexpr = %d", arg->nexpr);
1649 RASTER_DEBUGF(4, "expr @ %p", arg->expr);
1650
1651 for (i = 0; i < arg->nexpr; i++) {
1652 arg->expr[i] = rtalloc(sizeof(struct rt_reclassexpr_t));
1653 if (arg->expr[i] == NULL) {
1654 rterror("rt_raster_colormap: Could not allocate memory for reclass expression");
1656 return NULL;
1657 }
1658 }
1659
1660 /* reclassify bands */
1661 /* by # of colors */
1662 for (i = 0; i < colormap->ncolor; i++) {
1663 k = 0;
1664
1665 /* handle NODATA entry first */
1666 if (arg->nodataentry != NULL) {
1667 arg->expr[k]->src.min = arg->nodataentry->value;
1668 arg->expr[k]->src.max = arg->nodataentry->value;
1669 arg->expr[k]->src.inc_min = 1;
1670 arg->expr[k]->src.inc_max = 1;
1671 arg->expr[k]->src.exc_min = 0;
1672 arg->expr[k]->src.exc_max = 0;
1673
1674 arg->expr[k]->dst.min = arg->nodataentry->color[i];
1675 arg->expr[k]->dst.max = arg->nodataentry->color[i];
1676
1677 arg->expr[k]->dst.inc_min = 1;
1678 arg->expr[k]->dst.inc_max = 1;
1679 arg->expr[k]->dst.exc_min = 0;
1680 arg->expr[k]->dst.exc_max = 0;
1681
1682 RASTER_DEBUGF(4, "NODATA expr[%d]->src (min, max, in, ix, en, ex) = (%f, %f, %d, %d, %d, %d)",
1683 k,
1684 arg->expr[k]->src.min,
1685 arg->expr[k]->src.max,
1686 arg->expr[k]->src.inc_min,
1687 arg->expr[k]->src.inc_max,
1688 arg->expr[k]->src.exc_min,
1689 arg->expr[k]->src.exc_max
1690 );
1691 RASTER_DEBUGF(4, "NODATA expr[%d]->dst (min, max, in, ix, en, ex) = (%f, %f, %d, %d, %d, %d)",
1692 k,
1693 arg->expr[k]->dst.min,
1694 arg->expr[k]->dst.max,
1695 arg->expr[k]->dst.inc_min,
1696 arg->expr[k]->dst.inc_max,
1697 arg->expr[k]->dst.exc_min,
1698 arg->expr[k]->dst.exc_max
1699 );
1700
1701 k++;
1702 }
1703
1704 /* by non-NODATA entry */
1705 for (j = 0; j < arg->npos; j++) {
1706 if (colormap->method == CM_INTERPOLATE) {
1707 if (j == arg->npos - 1)
1708 continue;
1709
1710 arg->expr[k]->src.min = colormap->entry[arg->pos[j + 1]].value;
1711 arg->expr[k]->src.inc_min = 1;
1712 arg->expr[k]->src.exc_min = 0;
1713
1714 arg->expr[k]->src.max = colormap->entry[arg->pos[j]].value;
1715 arg->expr[k]->src.inc_max = 1;
1716 arg->expr[k]->src.exc_max = 0;
1717
1718 arg->expr[k]->dst.min = colormap->entry[arg->pos[j + 1]].color[i];
1719 arg->expr[k]->dst.max = colormap->entry[arg->pos[j]].color[i];
1720
1721 arg->expr[k]->dst.inc_min = 1;
1722 arg->expr[k]->dst.exc_min = 0;
1723
1724 arg->expr[k]->dst.inc_max = 1;
1725 arg->expr[k]->dst.exc_max = 0;
1726 }
1727 else if (colormap->method == CM_NEAREST) {
1728
1729 /* NOT last entry */
1730 if (j != arg->npos - 1) {
1731 arg->expr[k]->src.min = ((colormap->entry[arg->pos[j]].value - colormap->entry[arg->pos[j + 1]].value) / 2.) + colormap->entry[arg->pos[j + 1]].value;
1732 arg->expr[k]->src.inc_min = 0;
1733 arg->expr[k]->src.exc_min = 0;
1734 }
1735 /* last entry */
1736 else {
1737 arg->expr[k]->src.min = colormap->entry[arg->pos[j]].value;
1738 arg->expr[k]->src.inc_min = 1;
1739 arg->expr[k]->src.exc_min = 1;
1740 }
1741
1742 /* NOT first entry */
1743 if (j > 0) {
1744 arg->expr[k]->src.max = arg->expr[k - 1]->src.min;
1745 arg->expr[k]->src.inc_max = 1;
1746 arg->expr[k]->src.exc_max = 0;
1747 }
1748 /* first entry */
1749 else {
1750 arg->expr[k]->src.max = colormap->entry[arg->pos[j]].value;
1751 arg->expr[k]->src.inc_max = 1;
1752 arg->expr[k]->src.exc_max = 1;
1753 }
1754
1755 arg->expr[k]->dst.min = colormap->entry[arg->pos[j]].color[i];
1756 arg->expr[k]->dst.inc_min = 1;
1757 arg->expr[k]->dst.exc_min = 0;
1758
1759 arg->expr[k]->dst.max = colormap->entry[arg->pos[j]].color[i];
1760 arg->expr[k]->dst.inc_max = 1;
1761 arg->expr[k]->dst.exc_max = 0;
1762 }
1763 else if (colormap->method == CM_EXACT) {
1764 arg->expr[k]->src.min = colormap->entry[arg->pos[j]].value;
1765 arg->expr[k]->src.inc_min = 1;
1766 arg->expr[k]->src.exc_min = 0;
1767
1768 arg->expr[k]->src.max = colormap->entry[arg->pos[j]].value;
1769 arg->expr[k]->src.inc_max = 1;
1770 arg->expr[k]->src.exc_max = 0;
1771
1772 arg->expr[k]->dst.min = colormap->entry[arg->pos[j]].color[i];
1773 arg->expr[k]->dst.inc_min = 1;
1774 arg->expr[k]->dst.exc_min = 0;
1775
1776 arg->expr[k]->dst.max = colormap->entry[arg->pos[j]].color[i];
1777 arg->expr[k]->dst.inc_max = 1;
1778 arg->expr[k]->dst.exc_max = 0;
1779 }
1780
1781 RASTER_DEBUGF(4, "expr[%d]->src (min, max, in, ix, en, ex) = (%f, %f, %d, %d, %d, %d)",
1782 k,
1783 arg->expr[k]->src.min,
1784 arg->expr[k]->src.max,
1785 arg->expr[k]->src.inc_min,
1786 arg->expr[k]->src.inc_max,
1787 arg->expr[k]->src.exc_min,
1788 arg->expr[k]->src.exc_max
1789 );
1790
1791 RASTER_DEBUGF(4, "expr[%d]->dst (min, max, in, ix, en, ex) = (%f, %f, %d, %d, %d, %d)",
1792 k,
1793 arg->expr[k]->dst.min,
1794 arg->expr[k]->dst.max,
1795 arg->expr[k]->dst.inc_min,
1796 arg->expr[k]->dst.inc_max,
1797 arg->expr[k]->dst.exc_min,
1798 arg->expr[k]->dst.exc_max
1799 );
1800
1801 k++;
1802 }
1803
1804 /* EXACT has one last expression for catching all uncaught values */
1805 if (colormap->method == CM_EXACT) {
1806 arg->expr[k]->src.min = 0;
1807 arg->expr[k]->src.inc_min = 1;
1808 arg->expr[k]->src.exc_min = 1;
1809
1810 arg->expr[k]->src.max = 0;
1811 arg->expr[k]->src.inc_max = 1;
1812 arg->expr[k]->src.exc_max = 1;
1813
1814 arg->expr[k]->dst.min = 0;
1815 arg->expr[k]->dst.inc_min = 1;
1816 arg->expr[k]->dst.exc_min = 0;
1817
1818 arg->expr[k]->dst.max = 0;
1819 arg->expr[k]->dst.inc_max = 1;
1820 arg->expr[k]->dst.exc_max = 0;
1821
1822 RASTER_DEBUGF(4, "expr[%d]->src (min, max, in, ix, en, ex) = (%f, %f, %d, %d, %d, %d)",
1823 k,
1824 arg->expr[k]->src.min,
1825 arg->expr[k]->src.max,
1826 arg->expr[k]->src.inc_min,
1827 arg->expr[k]->src.inc_max,
1828 arg->expr[k]->src.exc_min,
1829 arg->expr[k]->src.exc_max
1830 );
1831
1832 RASTER_DEBUGF(4, "expr[%d]->dst (min, max, in, ix, en, ex) = (%f, %f, %d, %d, %d, %d)",
1833 k,
1834 arg->expr[k]->dst.min,
1835 arg->expr[k]->dst.max,
1836 arg->expr[k]->dst.inc_min,
1837 arg->expr[k]->dst.inc_max,
1838 arg->expr[k]->dst.exc_min,
1839 arg->expr[k]->dst.exc_max
1840 );
1841
1842 k++;
1843 }
1844
1845 /* call rt_band_reclass */
1846 arg->band = rt_band_reclass(band, PT_8BUI, 0, 0, arg->expr, arg->nexpr);
1847 if (arg->band == NULL) {
1848 rterror("rt_raster_colormap: Could not reclassify band");
1850 return NULL;
1851 }
1852
1853 /* add reclassified band to raster */
1854 if (rt_raster_add_band(arg->raster, arg->band, rt_raster_get_num_bands(arg->raster)) < 0) {
1855 rterror("rt_raster_colormap: Could not add reclassified band to output raster");
1857 return NULL;
1858 }
1859 }
1860
1861 rtnraster = arg->raster;
1862 arg->raster = NULL;
1864
1865 return rtnraster;
1866}
void rt_band_init_value(rt_band band, double initval)
Fill in the cells of a band with a starting value frequently used to init with nodata value.
Definition rt_band.c:113
rt_band rt_band_new_inline(uint16_t width, uint16_t height, rt_pixtype pixtype, uint32_t hasnodata, double nodataval, uint8_t *data)
Create an in-db rt_band with no data.
Definition rt_band.c:64
void rt_band_set_ownsdata_flag(rt_band band, int flag)
Definition rt_band.c:826
void rterror(const char *fmt,...) __attribute__((format(printf
Wrappers used for reporting errors and info.
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition rt_context.c:191
uint16_t rt_band_get_width(rt_band band)
Return width of this band.
Definition rt_band.c:799
#define RASTER_DEBUG(level, msg)
Definition librtcore.h:304
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition rt_raster.c:360
#define RASTER_DEBUGF(level, msg,...)
Definition librtcore.h:308
double rt_raster_get_x_skew(rt_raster raster)
Get skew about the X axis.
Definition rt_raster.c:185
double rt_raster_get_x_offset(rt_raster raster)
Get raster x offset, in projection units.
Definition rt_raster.c:217
int rt_raster_generate_new_band(rt_raster raster, rt_pixtype pixtype, double initialvalue, uint32_t hasnodata, double nodatavalue, int index)
Generate a new inline band and add it to a raster.
Definition rt_raster.c:489
void void rtinfo(const char *fmt,...) __attribute__((format(printf
void void void rtwarn(const char *fmt,...) __attribute__((format(printf
void rt_raster_set_scale(rt_raster raster, double scaleX, double scaleY)
Set scale in projection units.
Definition rt_raster.c:141
int rt_raster_add_band(rt_raster raster, rt_band band, int index)
Add band data to a raster.
Definition rt_raster.c:409
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition rt_band.c:833
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition rt_band.c:1551
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition rt_raster.c:86
rt_pixtype
Definition librtcore.h:188
@ PT_32BUI
Definition librtcore.h:197
@ PT_2BUI
Definition librtcore.h:190
@ PT_32BSI
Definition librtcore.h:196
@ PT_END
Definition librtcore.h:201
@ PT_4BUI
Definition librtcore.h:191
@ PT_1BB
Definition librtcore.h:189
@ PT_16BUI
Definition librtcore.h:195
@ PT_8BSI
Definition librtcore.h:192
@ PT_16BSI
Definition librtcore.h:194
@ PT_8BUI
Definition librtcore.h:193
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition rt_band.c:873
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition rt_raster.c:52
int rt_raster_has_band(rt_raster raster, int nband)
Return TRUE if the raster has a band of this number.
Definition rt_raster.c:1253
struct rt_pixel_t * rt_pixel
Definition librtcore.h:148
#define FLT_EQ(x, y)
Definition librtcore.h:2436
double rt_raster_get_x_scale(rt_raster raster)
Get scale X in projection units.
Definition rt_raster.c:154
double rt_band_get_min_value(rt_band band)
Returns the minimal possible value for the band according to the pixel type.
Definition rt_band.c:2082
rt_errorstate rt_band_set_pixel(rt_band band, int x, int y, double val, int *converted)
Set single pixel's value.
Definition rt_band.c:1140
rt_errorstate
Enum definitions.
Definition librtcore.h:181
@ ES_NONE
Definition librtcore.h:182
@ ES_ERROR
Definition librtcore.h:183
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition rt_band.c:499
void * rtrealloc(void *mem, size_t size)
Definition rt_context.c:199
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition rt_raster.c:376
rt_raster rt_raster_clone(rt_raster raster, uint8_t deep)
Clone an existing raster.
Definition rt_raster.c:1446
uint16_t rt_raster_get_height(rt_raster raster)
Definition rt_raster.c:133
rt_extenttype
Definition librtcore.h:204
@ ET_CUSTOM
Definition librtcore.h:210
@ ET_LAST
Definition librtcore.h:209
@ ET_INTERSECTION
Definition librtcore.h:205
@ ET_UNION
Definition librtcore.h:206
@ ET_SECOND
Definition librtcore.h:208
@ ET_FIRST
Definition librtcore.h:207
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition rt_band.c:2067
rt_errorstate rt_pixel_set_to_array(rt_pixel npixel, uint32_t count, rt_mask mask, int x, int y, uint16_t distancex, uint16_t distancey, double ***value, int ***nodata, int *dimx, int *dimy)
Definition rt_pixel.c:301
uint16_t rt_raster_get_width(rt_raster raster)
Definition rt_raster.c:125
uint32_t rt_band_get_nearest_pixel(rt_band band, int x, int y, uint16_t distancex, uint16_t distancey, int exclude_nodata_value, rt_pixel *npixels)
Get nearest pixel(s) with value (not NODATA) to specified pixel.
Definition rt_band.c:1709
void rtdealloc(void *mem)
Definition rt_context.c:206
rt_errorstate rt_raster_from_two_rasters(rt_raster rast1, rt_raster rast2, rt_extenttype extenttype, rt_raster *rtnraster, double *offset)
Definition rt_raster.c:3348
double rt_raster_get_y_scale(rt_raster raster)
Get scale Y in projection units.
Definition rt_raster.c:163
rt_errorstate rt_raster_same_alignment(rt_raster rast1, rt_raster rast2, int *aligned, char **reason)
double rt_raster_get_y_skew(rt_raster raster)
Get skew about the Y axis.
Definition rt_raster.c:194
uint16_t rt_band_get_height(rt_band band)
Return height of this band.
Definition rt_band.c:808
int rt_pixtype_size(rt_pixtype pixtype)
Return size in bytes of a value in the given pixtype.
Definition rt_pixel.c:40
double rt_raster_get_y_offset(rt_raster raster)
Get raster y offset, in projection units.
Definition rt_raster.c:226
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
Definition rt_raster.c:1240
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition rt_raster.c:385
This library is the generic raster handling section of PostGIS.
rt_band rt_band_reclass(rt_band srcband, rt_pixtype pixtype, uint32_t hasnodata, double nodataval, rt_reclassexpr *exprset, int exprcount)
Returns new band with values reclassified.
static void _rti_colormap_arg_destroy(_rti_colormap_arg arg)
static int _rti_iterator_arg_callback_init(_rti_iterator_arg _param)
rt_raster rt_raster_colormap(rt_raster raster, int nband, rt_colormap colormap)
Returns a new raster with up to four 8BUI bands (RGBA) from applying a colormap to the user-specified...
static int _rti_iterator_arg_empty_init(_rti_iterator_arg _param)
struct _rti_colormap_arg_t * _rti_colormap_arg
static _rti_iterator_arg _rti_iterator_arg_init()
struct _rti_iterator_arg_t * _rti_iterator_arg
static int _rti_iterator_arg_populate(_rti_iterator_arg _param, rt_iterator itrset, uint16_t itrcount, uint16_t distancex, uint16_t distancey, int *allnull, int *allempty)
rt_errorstate rt_raster_iterator(rt_iterator itrset, uint16_t itrcount, rt_extenttype extenttype, rt_raster customextent, rt_pixtype pixtype, uint8_t hasnodata, double nodataval, uint16_t distancex, uint16_t distancey, rt_mask mask, void *userarg, int(*callback)(rt_iterator_arg arg, void *userarg, double *value, int *nodata), rt_raster *rtnraster)
n-raster iterator.
rt_band rt_band_reclass_exact(rt_band srcband, rt_reclassmap map, uint32_t hasnodata, double nodataval)
Returns new band with values reclassified.
static void _rti_iterator_arg_callback_clean(_rti_iterator_arg _param)
static int rt_classpair_cmp(const void *aptr, const void *bptr)
static void _rti_iterator_arg_destroy(_rti_iterator_arg _param)
static double rt_band_reclass_round_integer(rt_pixtype pixtype, double nv)
static _rti_colormap_arg _rti_colormap_arg_init(rt_raster raster)
rt_colormap_entry nodataentry
rt_reclassexpr * expr
struct _rti_iterator_arg_t::@17 distance
struct _rti_iterator_arg_t::@16 band
struct _rti_iterator_arg_t::@18 dimension
struct _rti_iterator_arg_t::@19 empty
rt_iterator_arg arg
uint8_t color[4]
Definition librtcore.h:2700
double value
Definition librtcore.h:2699
int isnodata
Definition librtcore.h:2698
Definition librtcore.h:2697
rt_colormap_entry entry
Definition librtcore.h:2712
uint16_t nentry
Definition librtcore.h:2711
enum rt_colormap_t::@13 method
double *** values
Definition librtcore.h:2675
rt_raster raster
Definition librtcore.h:2659
uint8_t nbnodata
Definition librtcore.h:2661
double value
Definition librtcore.h:2540
uint8_t nodata
Definition librtcore.h:2539
rt_band * bands
Definition librtcore.h:2505
uint16_t numBands
Definition librtcore.h:2492
struct rt_reclassexpr_t::rt_reclassrange src
struct rt_reclassexpr_t::rt_reclassrange dst
rt_pixtype dsttype
Definition librtcore.h:2653
struct rt_classpair_t * pairs
Definition librtcore.h:2654