PostGIS  3.7.0dev-r@@SVN_REVISION@@
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 
37 static 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 
68 rt_band
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");
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  */
266 static int
267 rt_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 
287 rt_band
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__);
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 
401  int *isempty;
402  double **offset;
403  int *width;
404  int *height;
405 
406  struct {
408  int *hasnodata;
409  int *isnodata;
410  double *nodataval;
411  double *minval;
412  } band;
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;
427  } empty;
428 
430 };
431 
432 static _rti_iterator_arg
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 
470 static 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 
541 static 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 
680 static 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 
710 static 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 
756 static 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 
1007  _rti_iterator_arg_destroy(_param);
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 
1022  _rti_iterator_arg_destroy(_param);
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 
1050  _rti_iterator_arg_destroy(_param);
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 
1059  _rti_iterator_arg_destroy(_param);
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 
1095  _rti_iterator_arg_destroy(_param);
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 
1106  _rti_iterator_arg_destroy(_param);
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 
1125  _rti_iterator_arg_destroy(_param);
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 
1163  _rti_iterator_arg_destroy(_param);
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 
1179  _rti_iterator_arg_destroy(_param);
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 
1190  _rti_iterator_arg_destroy(_param);
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 
1203  _rti_iterator_arg_destroy(_param);
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 
1220  _rti_iterator_arg_destroy(_param);
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 
1280  status = rt_band_get_nearest_pixel(
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 
1290  _rti_iterator_arg_destroy(_param);
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");
1305  if (rt_band_get_pixel(
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 
1313  _rti_iterator_arg_destroy(_param);
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 
1345  _rti_iterator_arg_destroy(_param);
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 
1376  _rti_iterator_arg_destroy(_param);
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 
1397  _rti_iterator_arg_destroy(_param);
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 
1420  _rti_iterator_arg_destroy(_param);
1421  rt_band_destroy(rtnband);
1422  rt_raster_destroy(rtnrast);
1423 
1424  return ES_ERROR;
1425  }
1426  }
1427  }
1428 
1429  /* lots of cleanup */
1430  _rti_iterator_arg_destroy(_param);
1431 
1432  *rtnraster = rtnrast;
1433  return ES_NONE;
1434 }
1435 
1436 /******************************************************************************
1437 * rt_raster_colormap()
1438 ******************************************************************************/
1439 
1444 
1447  double nodataval;
1448 
1449  int nexpr;
1451 
1452  int npos;
1453  uint16_t *pos;
1454 
1455 };
1456 
1457 static _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 
1492 static 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)
1503  }
1504 
1505  rt_raster_destroy(arg->raster);
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 */
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 
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 */
1571  if (arg == NULL) {
1572  rterror("rt_raster_colormap: Could not initialize internal variables");
1573  return NULL;
1574  }
1575 
1576  /* handle NODATA */
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:112
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:63
void rt_band_set_ownsdata_flag(rt_band band, int flag)
Definition: rt_band.c:818
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:791
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:302
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:306
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:825
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_band.c:1527
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:86
rt_pixtype
Definition: librtcore.h:187
@ PT_32BUI
Definition: librtcore.h:196
@ PT_2BUI
Definition: librtcore.h:189
@ PT_32BSI
Definition: librtcore.h:195
@ PT_END
Definition: librtcore.h:199
@ PT_4BUI
Definition: librtcore.h:190
@ PT_1BB
Definition: librtcore.h:188
@ PT_16BUI
Definition: librtcore.h:194
@ PT_8BSI
Definition: librtcore.h:191
@ PT_16BSI
Definition: librtcore.h:193
@ PT_8BUI
Definition: librtcore.h:192
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition: rt_band.c:865
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:2424
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:2053
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:1125
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:491
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:202
@ ET_CUSTOM
Definition: librtcore.h:208
@ ET_LAST
Definition: librtcore.h:207
@ ET_INTERSECTION
Definition: librtcore.h:203
@ ET_UNION
Definition: librtcore.h:204
@ ET_SECOND
Definition: librtcore.h:206
@ ET_FIRST
Definition: librtcore.h:205
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:2038
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:288
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:1680
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:800
int rt_pixtype_size(rt_pixtype pixtype)
Return size in bytes of a value in the given pixtype.
Definition: rt_pixel.c:39
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.
if(!(yy_init))
int value
Definition: genraster.py:62
band
Definition: ovdump.py:58
nband
Definition: pixval.py:54
raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
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.
Definition: rt_mapalgebra.c:69
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)
Definition: rt_mapalgebra.c:38
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:2688
double value
Definition: librtcore.h:2687
int isnodata
Definition: librtcore.h:2686
Definition: librtcore.h:2685
rt_colormap_entry entry
Definition: librtcore.h:2700
uint16_t nentry
Definition: librtcore.h:2699
enum rt_colormap_t::@13 method
double *** values
Definition: librtcore.h:2663
uint32_t columns
Definition: librtcore.h:2659
uint16_t rasters
Definition: librtcore.h:2655
rt_raster raster
Definition: librtcore.h:2647
uint8_t nbnodata
Definition: librtcore.h:2649
double value
Definition: librtcore.h:2528
uint8_t nodata
Definition: librtcore.h:2527
rt_band * bands
Definition: librtcore.h:2493
uint16_t numBands
Definition: librtcore.h:2480
struct rt_reclassexpr_t::rt_reclassrange src
struct rt_reclassexpr_t::rt_reclassrange dst
uint32_t count
Definition: librtcore.h:2639
rt_pixtype dsttype
Definition: librtcore.h:2641
struct rt_classpair_t * pairs
Definition: librtcore.h:2642