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