PostGIS  2.5.7dev-r@@SVN_REVISION@@
rt_band.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) 2018 Bborie Park <dustymugs@gmail.com>
7  * Copyright (C) 2011-2013 Regents of the University of California
8  * <bkpark@ucdavis.edu>
9  * Copyright (C) 2010-2011 Jorge Arevalo <jorge.arevalo@deimos-space.com>
10  * Copyright (C) 2010-2011 David Zwarg <dzwarg@azavea.com>
11  * Copyright (C) 2009-2011 Pierre Racine <pierre.racine@sbf.ulaval.ca>
12  * Copyright (C) 2009-2011 Mateusz Loskot <mateusz@loskot.net>
13  * Copyright (C) 2008-2009 Sandro Santilli <strk@kbt.io>
14  *
15  * This program is free software; you can redistribute it and/or
16  * modify it under the terms of the GNU General Public License
17  * as published by the Free Software Foundation; either version 2
18  * of the License, or (at your option) any later version.
19  *
20  * This program is distributed in the hope that it will be useful,
21  * but WITHOUT ANY WARRANTY; without even the implied warranty of
22  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23  * GNU General Public License for more details.
24  *
25  * You should have received a copy of the GNU General Public License
26  * along with this program; if not, write to the Free Software Foundation,
27  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
28  *
29  */
30 
31 // For stat64()
32 #define _LARGEFILE64_SOURCE 1
33 
34 #include <stdio.h>
35 
36 #include "librtcore.h"
37 #include "librtcore_internal.h"
38 
39 #include "cpl_vsi.h"
40 #include "gdal_vrt.h"
41 
62 rt_band
64  uint16_t width, uint16_t height,
65  rt_pixtype pixtype,
66  uint32_t hasnodata, double nodataval,
67  uint8_t* data
68 ) {
69  rt_band band = NULL;
70 
71  assert(NULL != data);
72 
73  band = rtalloc(sizeof(struct rt_band_t));
74  if (band == NULL) {
75  rterror("rt_band_new_inline: Out of memory allocating rt_band");
76  return NULL;
77  }
78 
79  RASTER_DEBUGF(3, "Created rt_band @ %p with pixtype %s", band, rt_pixtype_name(pixtype));
80 
81  band->pixtype = pixtype;
82  band->offline = 0;
83  band->width = width;
84  band->height = height;
85  band->hasnodata = hasnodata ? 1 : 0;
86  band->isnodata = FALSE; /* we don't know what is in data, so must be FALSE */
87  band->nodataval = 0;
88  band->data.mem = data;
89  band->ownsdata = 0; /* we do NOT own this data!!! */
90  band->raster = NULL;
91 
92  RASTER_DEBUGF(3, "Created rt_band with dimensions %d x %d", band->width, band->height);
93 
94  /* properly set nodataval as it may need to be constrained to the data type */
95  if (hasnodata && rt_band_set_nodata(band, nodataval, NULL) != ES_NONE) {
96  rterror("rt_band_new_inline: Could not set NODATA value");
98  return NULL;
99  }
100 
101  return band;
102 }
103 
123 rt_band
125  uint16_t width, uint16_t height,
126  rt_pixtype pixtype,
127  uint32_t hasnodata, double nodataval,
128  uint8_t bandNum, const char* path
129 ) {
130  rt_band band = NULL;
131  int pathlen = 0;
132 
133  assert(NULL != path);
134 
135  band = rtalloc(sizeof(struct rt_band_t));
136  if (band == NULL) {
137  rterror("rt_band_new_offline: Out of memory allocating rt_band");
138  return NULL;
139  }
140 
141  RASTER_DEBUGF(3, "Created rt_band @ %p with pixtype %s",
142  band, rt_pixtype_name(pixtype)
143  );
144 
145  band->pixtype = pixtype;
146  band->offline = 1;
147  band->width = width;
148  band->height = height;
149  band->hasnodata = hasnodata ? 1 : 0;
150  band->nodataval = 0;
151  band->isnodata = FALSE; /* we don't know if the offline band is NODATA */
152  band->ownsdata = 0; /* offline, flag is useless as all offline data cache is owned internally */
153  band->raster = NULL;
154 
155  /* properly set nodataval as it may need to be constrained to the data type */
156  if (hasnodata && rt_band_set_nodata(band, nodataval, NULL) != ES_NONE) {
157  rterror("rt_band_new_offline: Could not set NODATA value");
159  return NULL;
160  }
161 
162  band->data.offline.bandNum = bandNum;
163 
164  /* memory for data.offline.path is managed internally */
165  pathlen = strlen(path);
166  band->data.offline.path = rtalloc(sizeof(char) * (pathlen + 1));
167  if (band->data.offline.path == NULL) {
168  rterror("rt_band_new_offline: Out of memory allocating offline path");
170  return NULL;
171  }
172  memcpy(band->data.offline.path, path, pathlen);
173  band->data.offline.path[pathlen] = '\0';
174 
175  band->data.offline.mem = NULL;
176 
177  return band;
178 }
179 
198 rt_band
200  uint16_t width,
201  uint16_t height,
202  int hasnodata,
203  double nodataval,
204  uint8_t bandNum,
205  const char* path,
206  int force
207 ) {
208  GDALDatasetH hdsSrc = NULL;
209  int nband = 0;
210  GDALRasterBandH hbandSrc = NULL;
211 
212  GDALDataType gdpixtype;
213  rt_pixtype pt = PT_END;
214 
215  /* open outdb raster file */
217  hdsSrc = rt_util_gdal_open(path, GA_ReadOnly, 1);
218  if (hdsSrc == NULL && !force) {
219  rterror("rt_band_new_offline_from_path: Cannot open offline raster: %s", path);
220  return NULL;
221  }
222 
223  nband = GDALGetRasterCount(hdsSrc);
224  if (!nband && !force) {
225  rterror("rt_band_new_offline_from_path: No bands found in offline raster: %s", path);
226  GDALClose(hdsSrc);
227  return NULL;
228  }
229  /* bandNum is 1-based */
230  else if (bandNum > nband && !force) {
231  rterror(
232  "rt_band_new_offline_from_path: Specified band %d not found in offline raster: %s",
233  bandNum,
234  path
235  );
236  GDALClose(hdsSrc);
237  return NULL;
238  }
239 
240  hbandSrc = GDALGetRasterBand(hdsSrc, bandNum);
241  if (hbandSrc == NULL && !force) {
242  rterror(
243  "rt_band_new_offline_from_path: Cannot get band %d from GDAL dataset",
244  bandNum
245  );
246  GDALClose(hdsSrc);
247  return NULL;
248  }
249 
250  gdpixtype = GDALGetRasterDataType(hbandSrc);
251  pt = rt_util_gdal_datatype_to_pixtype(gdpixtype);
252  if (pt == PT_END && !force) {
253  rterror(
254  "rt_band_new_offline_from_path: Unsupported pixel type %s of band %d from GDAL dataset",
255  GDALGetDataTypeName(gdpixtype),
256  bandNum
257  );
258  GDALClose(hdsSrc);
259  return NULL;
260  }
261 
262  /* use out-db band's nodata value if nodataval not already set */
263  if (!hasnodata)
264  nodataval = GDALGetRasterNoDataValue(hbandSrc, &hasnodata);
265 
266  GDALClose(hdsSrc);
267 
268  return rt_band_new_offline(
269  width, height,
270  pt,
271  hasnodata, nodataval,
272  bandNum - 1, path
273  );
274 }
275 
286 rt_band
288  rt_band rtn = NULL;
289 
290  assert(band != NULL);
291 
292  /* offline */
293  if (band->offline) {
294  rtn = rt_band_new_offline(
295  band->width, band->height,
296  band->pixtype,
297  band->hasnodata, band->nodataval,
298  band->data.offline.bandNum, (const char *) band->data.offline.path
299  );
300  }
301  /* online */
302  else {
303  uint8_t *data = NULL;
304  data = rtalloc(rt_pixtype_size(band->pixtype) * band->width * band->height);
305  if (data == NULL) {
306  rterror("rt_band_duplicate: Out of memory allocating online band data");
307  return NULL;
308  }
309  memcpy(data, band->data.mem, rt_pixtype_size(band->pixtype) * band->width * band->height);
310 
311  rtn = rt_band_new_inline(
312  band->width, band->height,
313  band->pixtype,
314  band->hasnodata, band->nodataval,
315  data
316  );
317  rt_band_set_ownsdata_flag(rtn, 1); /* we DO own this data!!! */
318  }
319 
320  if (rtn == NULL) {
321  rterror("rt_band_duplicate: Could not copy band");
322  return NULL;
323  }
324 
325  return rtn;
326 }
327 
328 int
330  assert(NULL != band);
331  return band->offline ? 1 : 0;
332 }
333 
339 void
341  if (band == NULL)
342  return;
343 
344  RASTER_DEBUGF(3, "Destroying rt_band @ %p", band);
345 
346  /* offline band */
347  if (band->offline) {
348  /* memory cache */
349  if (band->data.offline.mem != NULL)
350  rtdealloc(band->data.offline.mem);
351  /* offline file path */
352  if (band->data.offline.path != NULL)
353  rtdealloc(band->data.offline.path);
354  }
355  /* inline band and band owns the data */
356  else if (band->data.mem != NULL && band->ownsdata)
357  rtdealloc(band->data.mem);
358 
359  rtdealloc(band);
360 }
361 
362 const char*
364 
365  assert(NULL != band);
366 
367 
368  if (!band->offline) {
369  RASTER_DEBUG(3, "rt_band_get_ext_path: Band is not offline");
370  return NULL;
371  }
372  return band->data.offline.path;
373 }
374 
377  assert(NULL != band);
378  assert(NULL != bandnum);
379 
380  *bandnum = 0;
381 
382  if (!band->offline) {
383  RASTER_DEBUG(3, "rt_band_get_ext_band_num: Band is not offline");
384  return ES_ERROR;
385  }
386 
387  *bandnum = band->data.offline.bandNum;
388 
389  return ES_NONE;
390 }
391 
399 void *
401  assert(NULL != band);
402 
403  if (band->offline) {
404  if (band->data.offline.mem != NULL)
405  return band->data.offline.mem;
406 
408  return NULL;
409  else
410  return band->data.offline.mem;
411  }
412  else
413  return band->data.mem;
414 }
415 
416 /* variable for PostgreSQL GUC: postgis.enable_outdb_rasters */
418 
430  GDALDatasetH hdsSrc = NULL;
431  int nband = 0;
432  VRTDatasetH hdsDst = NULL;
433  VRTSourcedRasterBandH hbandDst = NULL;
434  double ogt[6] = {0};
435  double offset[2] = {0};
436 
437  rt_raster _rast = NULL;
438  rt_band _band = NULL;
439  int aligned = 0;
440  int err = ES_NONE;
441 
442  assert(band != NULL);
443  assert(band->raster != NULL);
444 
445  if (!band->offline) {
446  rterror("rt_band_load_offline_data: Band is not offline");
447  return ES_ERROR;
448  }
449  else if (!strlen(band->data.offline.path)) {
450  rterror("rt_band_load_offline_data: Offline band does not a have a specified file");
451  return ES_ERROR;
452  }
453 
454  /* offline_data is disabled */
455  if (!enable_outdb_rasters) {
456  rterror("rt_band_load_offline_data: Access to offline bands disabled");
457  return ES_ERROR;
458  }
459 
461  hdsSrc = rt_util_gdal_open(band->data.offline.path, GA_ReadOnly, 1);
462  if (hdsSrc == NULL) {
463  rterror("rt_band_load_offline_data: Cannot open offline raster: %s", band->data.offline.path);
464  return ES_ERROR;
465  }
466 
467  /* # of bands */
468  nband = GDALGetRasterCount(hdsSrc);
469  if (!nband) {
470  rterror("rt_band_load_offline_data: No bands found in offline raster: %s", band->data.offline.path);
471  GDALClose(hdsSrc);
472  return ES_ERROR;
473  }
474  /* bandNum is 0-based */
475  else if (band->data.offline.bandNum + 1 > nband) {
476  rterror("rt_band_load_offline_data: Specified band %d not found in offline raster: %s", band->data.offline.bandNum, band->data.offline.path);
477  GDALClose(hdsSrc);
478  return ES_ERROR;
479  }
480 
481  /* get offline raster's geotransform */
482  if (GDALGetGeoTransform(hdsSrc, ogt) != CE_None) {
483  RASTER_DEBUG(4, "Using default geotransform matrix (0, 1, 0, 0, 0, -1)");
484  ogt[0] = 0;
485  ogt[1] = 1;
486  ogt[2] = 0;
487  ogt[3] = 0;
488  ogt[4] = 0;
489  ogt[5] = -1;
490  }
491  RASTER_DEBUGF(3, "Offline geotransform (%f, %f, %f, %f, %f, %f)",
492  ogt[0], ogt[1], ogt[2], ogt[3], ogt[4], ogt[5]);
493 
494  /* are rasters aligned? */
495  _rast = rt_raster_new(1, 1);
497  rt_raster_set_srid(_rast, band->raster->srid);
498  err = rt_raster_same_alignment(band->raster, _rast, &aligned, NULL);
499  rt_raster_destroy(_rast);
500 
501  if (err != ES_NONE) {
502  rterror("rt_band_load_offline_data: Could not test alignment of in-db representation of out-db raster");
503  GDALClose(hdsSrc);
504  return ES_ERROR;
505  }
506  else if (!aligned) {
507  rtwarn("The in-db representation of the out-db raster is not aligned. Band data may be incorrect");
508  }
509 
510  /* get offsets */
512  band->raster,
513  ogt[0], ogt[3],
514  &(offset[0]), &(offset[1]),
515  NULL
516  );
517 
518  RASTER_DEBUGF(4, "offsets: (%f, %f)", offset[0], offset[1]);
519 
520  /* create VRT dataset */
521  hdsDst = VRTCreate(band->width, band->height);
522  GDALSetGeoTransform(hdsDst, ogt);
523  /*
524  GDALSetDescription(hdsDst, "/tmp/offline.vrt");
525  */
526 
527  /* add band as simple sources */
528  GDALAddBand(hdsDst, rt_util_pixtype_to_gdal_datatype(band->pixtype), NULL);
529  hbandDst = (VRTSourcedRasterBandH) GDALGetRasterBand(hdsDst, 1);
530 
531  if (band->hasnodata)
532  GDALSetRasterNoDataValue(hbandDst, band->nodataval);
533 
534  VRTAddSimpleSource(
535  hbandDst, GDALGetRasterBand(hdsSrc, band->data.offline.bandNum + 1),
536  fabs(offset[0]), fabs(offset[1]),
537  band->width, band->height,
538  0, 0,
539  band->width, band->height,
540  "near", VRT_NODATA_UNSET
541  );
542 
543  /* make sure VRT reflects all changes */
544  VRTFlushCache(hdsDst);
545 
546  /* convert VRT dataset to rt_raster */
547  _rast = rt_raster_from_gdal_dataset(hdsDst);
548 
549  GDALClose(hdsDst);
550  GDALClose(hdsSrc);
551  /*
552  {
553  FILE *fp;
554  fp = fopen("/tmp/gdal_open_files.log", "w");
555  GDALDumpOpenDatasets(fp);
556  fclose(fp);
557  }
558  */
559 
560  if (_rast == NULL) {
561  rterror("rt_band_load_offline_data: Cannot load data from offline raster: %s", band->data.offline.path);
562  return ES_ERROR;
563  }
564 
565  _band = rt_raster_get_band(_rast, 0);
566  if (_band == NULL) {
567  rterror("rt_band_load_offline_data: Cannot load data from offline raster: %s", band->data.offline.path);
568  rt_raster_destroy(_rast);
569  return ES_ERROR;
570  }
571 
572  /* band->data.offline.mem not NULL, free first */
573  if (band->data.offline.mem != NULL) {
574  rtdealloc(band->data.offline.mem);
575  band->data.offline.mem = NULL;
576  }
577 
578  band->data.offline.mem = _band->data.mem;
579 
580  rtdealloc(_band); /* cannot use rt_band_destroy */
581  rt_raster_destroy(_rast);
582 
583  return ES_NONE;
584 }
585 
587  VSIStatBufL sStat;
588 
589  assert(NULL != band);
590  if (!band->offline) {
591  rterror("rt_band_get_file_size: Band is not offline");
592  return 0;
593  }
594  /* offline_data is disabled */
595  if (!enable_outdb_rasters) {
596  rterror("rt_band_get_file_size: Access to offline bands disabled");
597  return 0;
598  }
599 
600  if( VSIStatL(band->data.offline.path, &sStat) != 0 ) {
601  rterror("rt_band_get_file_size: Cannot access file");
602  return 0;
603  }
604 
605  return sStat.st_size;
606 }
607 
609  VSIStatBufL sStat;
610 
611  assert(NULL != band);
612  if (!band->offline) {
613  rterror("rt_band_get_file_timestamp: Band is not offline");
614  return 0;
615  }
616  /* offline_data is disabled */
617  if (!enable_outdb_rasters) {
618  rterror("rt_band_get_file_timestamp: Access to offline bands disabled");
619  return 0;
620  }
621 
622  if( VSIStatL(band->data.offline.path, &sStat) != 0 ) {
623  rterror("rt_band_get_file_timestamp: Cannot access file");
624  return 0;
625  }
626 
627  return sStat.st_mtime;
628 }
629 
632 
633  assert(NULL != band);
634 
635 
636  return band->pixtype;
637 }
638 
639 uint16_t
641 
642  assert(NULL != band);
643 
644 
645  return band->width;
646 }
647 
648 uint16_t
650 
651  assert(NULL != band);
652 
653 
654  return band->height;
655 }
656 
657 /* Get ownsdata flag */
658 int
660  assert(NULL != band);
661 
662  return band->ownsdata ? 1 : 0;
663 }
664 
665 /* set ownsdata flag */
666 void
668  assert(NULL != band);
669 
670  band->ownsdata = flag ? 1 : 0;
671 }
672 
673 int
675  assert(NULL != band);
676 
677  return band->hasnodata ? 1 : 0;
678 }
679 
680 void
682 
683  assert(NULL != band);
684 
685  band->hasnodata = (flag) ? 1 : 0;
686 
687  /* isnodata depends on hasnodata */
688  if (!band->hasnodata && band->isnodata) {
689  RASTER_DEBUG(3, "Setting isnodata to FALSE as band no longer has NODATA");
690  band->isnodata = 0;
691  }
692 }
693 
696  assert(NULL != band);
697 
698  if (!band->hasnodata) {
699  /* silently permit setting isnodata flag to FALSE */
700  if (!flag)
701  band->isnodata = 0;
702  else {
703  rterror("rt_band_set_isnodata_flag: Cannot set isnodata flag as band has no NODATA");
704  return ES_ERROR;
705  }
706  }
707  else
708  band->isnodata = (flag) ? 1 : 0;
709 
710  return ES_NONE;
711 }
712 
713 int
715  assert(NULL != band);
716 
717  if (band->hasnodata)
718  return band->isnodata ? 1 : 0;
719  else
720  return 0;
721 }
722 
733 rt_band_set_nodata(rt_band band, double val, int *converted) {
734  rt_pixtype pixtype = PT_END;
735  int32_t checkvalint = 0;
736  uint32_t checkvaluint = 0;
737  float checkvalfloat = 0;
738  double checkvaldouble = 0;
739 
740  assert(NULL != band);
741 
742  if (converted != NULL)
743  *converted = 0;
744 
745  pixtype = band->pixtype;
746 
747  RASTER_DEBUGF(3, "rt_band_set_nodata: setting nodata value %g with band type %s", val, rt_pixtype_name(pixtype));
748 
749  /* return -1 on out of range */
750  switch (pixtype) {
751  case PT_1BB: {
752  band->nodataval = rt_util_clamp_to_1BB(val);
753  checkvalint = band->nodataval;
754  break;
755  }
756  case PT_2BUI: {
757  band->nodataval = rt_util_clamp_to_2BUI(val);
758  checkvalint = band->nodataval;
759  break;
760  }
761  case PT_4BUI: {
762  band->nodataval = rt_util_clamp_to_4BUI(val);
763  checkvalint = band->nodataval;
764  break;
765  }
766  case PT_8BSI: {
767  band->nodataval = rt_util_clamp_to_8BSI(val);
768  checkvalint = band->nodataval;
769  break;
770  }
771  case PT_8BUI: {
772  band->nodataval = rt_util_clamp_to_8BUI(val);
773  checkvalint = band->nodataval;
774  break;
775  }
776  case PT_16BSI: {
777  band->nodataval = rt_util_clamp_to_16BSI(val);
778  checkvalint = band->nodataval;
779  break;
780  }
781  case PT_16BUI: {
782  band->nodataval = rt_util_clamp_to_16BUI(val);
783  checkvalint = band->nodataval;
784  break;
785  }
786  case PT_32BSI: {
787  band->nodataval = rt_util_clamp_to_32BSI(val);
788  checkvalint = band->nodataval;
789  break;
790  }
791  case PT_32BUI: {
792  band->nodataval = rt_util_clamp_to_32BUI(val);
793  checkvaluint = band->nodataval;
794  break;
795  }
796  case PT_32BF: {
797  band->nodataval = rt_util_clamp_to_32F(val);
798  checkvalfloat = band->nodataval;
799  break;
800  }
801  case PT_64BF: {
802  band->nodataval = val;
803  checkvaldouble = band->nodataval;
804  break;
805  }
806  default: {
807  rterror("rt_band_set_nodata: Unknown pixeltype %d", pixtype);
808  band->hasnodata = 0;
809  return ES_ERROR;
810  }
811  }
812 
813  RASTER_DEBUGF(3, "rt_band_set_nodata: band->hasnodata = %d", band->hasnodata);
814  RASTER_DEBUGF(3, "rt_band_set_nodata: band->nodataval = %f", band->nodataval);
815  /* the nodata value was just set, so this band has NODATA */
816  band->hasnodata = 1;
817 
818  /* also set isnodata flag to false */
819  band->isnodata = 0;
820 
822  val,
823  checkvalint, checkvaluint,
824  checkvalfloat, checkvaldouble,
825  pixtype
826  ) && converted != NULL) {
827  *converted = 1;
828  }
829 
830  return ES_NONE;
831 }
832 
854  rt_band band,
855  int x, int y,
856  void *vals, uint32_t len
857 ) {
858  rt_pixtype pixtype = PT_END;
859  int size = 0;
860  uint8_t *data = NULL;
861  uint32_t offset = 0;
862 
863  assert(NULL != band);
864  assert(vals != NULL && len > 0);
865 
866  RASTER_DEBUGF(3, "length of values = %d", len);
867 
868  if (band->offline) {
869  rterror("rt_band_set_pixel_line not implemented yet for OFFDB bands");
870  return ES_ERROR;
871  }
872 
873  pixtype = band->pixtype;
874  size = rt_pixtype_size(pixtype);
875 
876  if (
877  x < 0 || x >= band->width ||
878  y < 0 || y >= band->height
879  ) {
880  rterror("rt_band_set_pixel_line: Coordinates out of range (%d, %d) vs (%d, %d)", x, y, band->width, band->height);
881  return ES_ERROR;
882  }
883 
885  offset = x + (y * band->width);
886  RASTER_DEBUGF(4, "offset = %d", offset);
887 
888  /* make sure len of values to copy don't exceed end of data */
889  if (len > (band->width * band->height) - offset) {
890  rterror("rt_band_set_pixel_line: Could not apply pixels as values length exceeds end of data");
891  return ES_ERROR;
892  }
893 
894  switch (pixtype) {
895  case PT_1BB:
896  case PT_2BUI:
897  case PT_4BUI:
898  case PT_8BUI:
899  case PT_8BSI: {
900  uint8_t *ptr = data;
901  ptr += offset;
902  memcpy(ptr, vals, size * len);
903  break;
904  }
905  case PT_16BUI: {
906  uint16_t *ptr = (uint16_t *) data;
907  ptr += offset;
908  memcpy(ptr, vals, size * len);
909  break;
910  }
911  case PT_16BSI: {
912  int16_t *ptr = (int16_t *) data;
913  ptr += offset;
914  memcpy(ptr, vals, size * len);
915  break;
916  }
917  case PT_32BUI: {
918  uint32_t *ptr = (uint32_t *) data;
919  ptr += offset;
920  memcpy(ptr, vals, size * len);
921  break;
922  }
923  case PT_32BSI: {
924  int32_t *ptr = (int32_t *) data;
925  ptr += offset;
926  memcpy(ptr, vals, size * len);
927  break;
928  }
929  case PT_32BF: {
930  float *ptr = (float *) data;
931  ptr += offset;
932  memcpy(ptr, vals, size * len);
933  break;
934  }
935  case PT_64BF: {
936  double *ptr = (double *) data;
937  ptr += offset;
938  memcpy(ptr, vals, size * len);
939  break;
940  }
941  default: {
942  rterror("rt_band_set_pixel_line: Unknown pixeltype %d", pixtype);
943  return ES_ERROR;
944  }
945  }
946 
947 #if POSTGIS_DEBUG_LEVEL > 0
948  {
949  double value;
950  rt_band_get_pixel(band, x, y, &value, NULL);
951  RASTER_DEBUGF(4, "pixel at (%d, %d) = %f", x, y, value);
952  }
953 #endif
954 
955  /* set band's isnodata flag to FALSE */
958 
959  return ES_NONE;
960 }
961 
975  rt_band band,
976  int x, int y,
977  double val,
978  int *converted
979 ) {
980  rt_pixtype pixtype = PT_END;
981  unsigned char* data = NULL;
982  uint32_t offset = 0;
983 
984  int32_t checkvalint = 0;
985  uint32_t checkvaluint = 0;
986  float checkvalfloat = 0;
987  double checkvaldouble = 0;
988 
989  assert(NULL != band);
990 
991  if (converted != NULL)
992  *converted = 0;
993 
994  if (band->offline) {
995  rterror("rt_band_set_pixel not implemented yet for OFFDB bands");
996  return ES_ERROR;
997  }
998 
999  pixtype = band->pixtype;
1000 
1001  if (
1002  x < 0 || x >= band->width ||
1003  y < 0 || y >= band->height
1004  ) {
1005  rterror("rt_band_set_pixel: Coordinates out of range");
1006  return ES_ERROR;
1007  }
1008 
1009  /* check that clamped value isn't clamped NODATA */
1010  if (band->hasnodata && pixtype != PT_64BF) {
1011  double newval;
1012  int corrected;
1013 
1014  rt_band_corrected_clamped_value(band, val, &newval, &corrected);
1015 
1016  if (corrected) {
1017 #if POSTGIS_RASTER_WARN_ON_TRUNCATION > 0
1018  rtwarn("Value for pixel %d x %d has been corrected as clamped value becomes NODATA", x, y);
1019 #endif
1020  val = newval;
1021 
1022  if (converted != NULL)
1023  *converted = 1;
1024  }
1025  }
1026 
1028  offset = x + (y * band->width);
1029 
1030  switch (pixtype) {
1031  case PT_1BB: {
1032  data[offset] = rt_util_clamp_to_1BB(val);
1033  checkvalint = data[offset];
1034  break;
1035  }
1036  case PT_2BUI: {
1037  data[offset] = rt_util_clamp_to_2BUI(val);
1038  checkvalint = data[offset];
1039  break;
1040  }
1041  case PT_4BUI: {
1042  data[offset] = rt_util_clamp_to_4BUI(val);
1043  checkvalint = data[offset];
1044  break;
1045  }
1046  case PT_8BSI: {
1047  data[offset] = rt_util_clamp_to_8BSI(val);
1048  checkvalint = (int8_t) data[offset];
1049  break;
1050  }
1051  case PT_8BUI: {
1052  data[offset] = rt_util_clamp_to_8BUI(val);
1053  checkvalint = data[offset];
1054  break;
1055  }
1056  case PT_16BSI: {
1057  int16_t *ptr = (int16_t*) data; /* we assume correct alignment */
1058  ptr[offset] = rt_util_clamp_to_16BSI(val);
1059  checkvalint = (int16_t) ptr[offset];
1060  break;
1061  }
1062  case PT_16BUI: {
1063  uint16_t *ptr = (uint16_t*) data; /* we assume correct alignment */
1064  ptr[offset] = rt_util_clamp_to_16BUI(val);
1065  checkvalint = ptr[offset];
1066  break;
1067  }
1068  case PT_32BSI: {
1069  int32_t *ptr = (int32_t*) data; /* we assume correct alignment */
1070  ptr[offset] = rt_util_clamp_to_32BSI(val);
1071  checkvalint = (int32_t) ptr[offset];
1072  break;
1073  }
1074  case PT_32BUI: {
1075  uint32_t *ptr = (uint32_t*) data; /* we assume correct alignment */
1076  ptr[offset] = rt_util_clamp_to_32BUI(val);
1077  checkvaluint = ptr[offset];
1078  break;
1079  }
1080  case PT_32BF: {
1081  float *ptr = (float*) data; /* we assume correct alignment */
1082  ptr[offset] = rt_util_clamp_to_32F(val);
1083  checkvalfloat = ptr[offset];
1084  break;
1085  }
1086  case PT_64BF: {
1087  double *ptr = (double*) data; /* we assume correct alignment */
1088  ptr[offset] = val;
1089  checkvaldouble = ptr[offset];
1090  break;
1091  }
1092  default: {
1093  rterror("rt_band_set_pixel: Unknown pixeltype %d", pixtype);
1094  return ES_ERROR;
1095  }
1096  }
1097 
1098  /* If the stored value is not NODATA, reset the isnodata flag */
1100  RASTER_DEBUG(3, "Band has a value that is not NODATA. Setting isnodata to FALSE");
1101  band->isnodata = FALSE;
1102  }
1103 
1104  /* Overflow checking */
1106  val,
1107  checkvalint, checkvaluint,
1108  checkvalfloat, checkvaldouble,
1109  pixtype
1110  ) && converted != NULL) {
1111  *converted = 1;
1112  }
1113 
1114  return ES_NONE;
1115 }
1116 
1138  rt_band band,
1139  int x, int y,
1140  uint16_t len,
1141  void **vals, uint16_t *nvals
1142 ) {
1143  uint8_t *_vals = NULL;
1144  int pixsize = 0;
1145  uint8_t *data = NULL;
1146  uint32_t offset = 0;
1147  uint16_t _nvals = 0;
1148  int maxlen = 0;
1149  uint8_t *ptr = NULL;
1150 
1151  assert(NULL != band);
1152  assert(vals != NULL && nvals != NULL);
1153 
1154  /* initialize to no values */
1155  *nvals = 0;
1156 
1157  if (
1158  x < 0 || x >= band->width ||
1159  y < 0 || y >= band->height
1160  ) {
1161  rtwarn("Attempting to get pixel values with out of range raster coordinates: (%d, %d)", x, y);
1162  return ES_ERROR;
1163  }
1164 
1165  if (len < 1)
1166  return ES_NONE;
1167 
1169  if (data == NULL) {
1170  rterror("rt_band_get_pixel_line: Cannot get band data");
1171  return ES_ERROR;
1172  }
1173 
1174  /* +1 for the nodata value */
1175  offset = x + (y * band->width);
1176  RASTER_DEBUGF(4, "offset = %d", offset);
1177 
1178  pixsize = rt_pixtype_size(band->pixtype);
1179  RASTER_DEBUGF(4, "pixsize = %d", pixsize);
1180 
1181  /* cap _nvals so that it doesn't overflow */
1182  _nvals = len;
1183  maxlen = band->width * band->height;
1184 
1185  if (((int) (offset + _nvals)) > maxlen) {
1186  _nvals = maxlen - offset;
1187  rtwarn("Limiting returning number values to %d", _nvals);
1188  }
1189  RASTER_DEBUGF(4, "_nvals = %d", _nvals);
1190 
1191  ptr = data + (offset * pixsize);
1192 
1193  _vals = rtalloc(_nvals * pixsize);
1194  if (_vals == NULL) {
1195  rterror("rt_band_get_pixel_line: Could not allocate memory for pixel values");
1196  return ES_ERROR;
1197  }
1198 
1199  /* copy pixels */
1200  memcpy(_vals, ptr, _nvals * pixsize);
1201 
1202  *vals = _vals;
1203  *nvals = _nvals;
1204 
1205  return ES_NONE;
1206 }
1207 
1222  rt_band band,
1223  int x, int y,
1224  double *value,
1225  int *nodata
1226 ) {
1227  rt_pixtype pixtype = PT_END;
1228  uint8_t* data = NULL;
1229  uint32_t offset = 0;
1230 
1231  assert(NULL != band);
1232  assert(NULL != value);
1233 
1234  /* set nodata to 0 */
1235  if (nodata != NULL)
1236  *nodata = 0;
1237 
1238  if (
1239  x < 0 || x >= band->width ||
1240  y < 0 || y >= band->height
1241  ) {
1242  rtwarn("Attempting to get pixel value with out of range raster coordinates: (%d, %d)", x, y);
1243  return ES_ERROR;
1244  }
1245 
1246  /* band is NODATA */
1247  if (band->isnodata) {
1248  RASTER_DEBUG(3, "Band's isnodata flag is TRUE. Returning NODATA value");
1249  *value = band->nodataval;
1250  if (nodata != NULL) *nodata = 1;
1251  return ES_NONE;
1252  }
1253 
1255  if (data == NULL) {
1256  rterror("rt_band_get_pixel: Cannot get band data");
1257  return ES_ERROR;
1258  }
1259 
1260  /* +1 for the nodata value */
1261  offset = x + (y * band->width);
1262 
1263  pixtype = band->pixtype;
1264 
1265  switch (pixtype) {
1266  case PT_1BB:
1267 #ifdef OPTIMIZE_SPACE
1268  {
1269  int byteOffset = offset / 8;
1270  int bitOffset = offset % 8;
1271  data += byteOffset;
1272 
1273  /* Bit to set is bitOffset into data */
1274  *value = getBits(data, val, 1, bitOffset);
1275  break;
1276  }
1277 #endif
1278  case PT_2BUI:
1279 #ifdef OPTIMIZE_SPACE
1280  {
1281  int byteOffset = offset / 4;
1282  int bitOffset = offset % 4;
1283  data += byteOffset;
1284 
1285  /* Bits to set start at bitOffset into data */
1286  *value = getBits(data, val, 2, bitOffset);
1287  break;
1288  }
1289 #endif
1290  case PT_4BUI:
1291 #ifdef OPTIMIZE_SPACE
1292  {
1293  int byteOffset = offset / 2;
1294  int bitOffset = offset % 2;
1295  data += byteOffset;
1296 
1297  /* Bits to set start at bitOffset into data */
1298  *value = getBits(data, val, 2, bitOffset);
1299  break;
1300  }
1301 #endif
1302  case PT_8BSI: {
1303  int8_t val = data[offset];
1304  *value = val;
1305  break;
1306  }
1307  case PT_8BUI: {
1308  uint8_t val = data[offset];
1309  *value = val;
1310  break;
1311  }
1312  case PT_16BSI: {
1313  int16_t *ptr = (int16_t*) data; /* we assume correct alignment */
1314  *value = ptr[offset];
1315  break;
1316  }
1317  case PT_16BUI: {
1318  uint16_t *ptr = (uint16_t*) data; /* we assume correct alignment */
1319  *value = ptr[offset];
1320  break;
1321  }
1322  case PT_32BSI: {
1323  int32_t *ptr = (int32_t*) data; /* we assume correct alignment */
1324  *value = ptr[offset];
1325  break;
1326  }
1327  case PT_32BUI: {
1328  uint32_t *ptr = (uint32_t*) data; /* we assume correct alignment */
1329  *value = ptr[offset];
1330  break;
1331  }
1332  case PT_32BF: {
1333  float *ptr = (float*) data; /* we assume correct alignment */
1334  *value = ptr[offset];
1335  break;
1336  }
1337  case PT_64BF: {
1338  double *ptr = (double*) data; /* we assume correct alignment */
1339  *value = ptr[offset];
1340  break;
1341  }
1342  default: {
1343  rterror("rt_band_get_pixel: Unknown pixeltype %d", pixtype);
1344  return ES_ERROR;
1345  }
1346  }
1347 
1348  /* set NODATA flag */
1349  if (band->hasnodata && nodata != NULL) {
1351  *nodata = 1;
1352  }
1353 
1354  return ES_NONE;
1355 }
1356 
1375  rt_band band,
1376  int x, int y,
1377  uint16_t distancex, uint16_t distancey,
1378  int exclude_nodata_value,
1379  rt_pixel *npixels
1380 ) {
1381  rt_pixel npixel = NULL;
1382  int extent[4] = {0};
1383  int max_extent[4] = {0};
1384  int d0 = 0;
1385  uint32_t distance[2] = {0};
1386  uint32_t _d[2] = {0};
1387  uint32_t i = 0;
1388  uint32_t j = 0;
1389  uint32_t k = 0;
1390  int _max = 0;
1391  int _x = 0;
1392  int _y = 0;
1393  int *_min = NULL;
1394  double pixval = 0;
1395  double minval = 0;
1396  uint32_t count = 0;
1397  int isnodata = 0;
1398 
1399  int inextent = 0;
1400 
1401  assert(NULL != band);
1402  assert(NULL != npixels);
1403 
1404  RASTER_DEBUG(3, "Starting");
1405 
1406  /* process distance */
1407  distance[0] = distancex;
1408  distance[1] = distancey;
1409 
1410  /* no distance, means get nearest pixels and return */
1411  if (!distance[0] && !distance[1])
1412  d0 = 1;
1413 
1414  RASTER_DEBUGF(4, "Selected pixel: %d x %d", x, y);
1415  RASTER_DEBUGF(4, "Distances: %d x %d", distance[0], distance[1]);
1416 
1417  /* shortcuts if outside band extent */
1418  if (
1419  exclude_nodata_value && (
1420  (x < 0 || x > band->width) ||
1421  (y < 0 || y > band->height)
1422  )
1423  ) {
1424  /* no distances specified, jump to pixel close to extent */
1425  if (d0) {
1426  if (x < 0)
1427  x = -1;
1428  else if (x > band->width)
1429  x = band->width;
1430 
1431  if (y < 0)
1432  y = -1;
1433  else if (y > band->height)
1434  y = band->height;
1435 
1436  RASTER_DEBUGF(4, "Moved selected pixel: %d x %d", x, y);
1437  }
1438  /*
1439  distances specified
1440  if distances won't capture extent of band, return 0
1441  */
1442  else if (
1443  ((x < 0 && (uint32_t) abs(x) > distance[0]) || (x - band->width >= (int)distance[0])) ||
1444  ((y < 0 && (uint32_t) abs(y) > distance[1]) || (y - band->height >= (int)distance[1]))
1445  ) {
1446  RASTER_DEBUG(4, "No nearest pixels possible for provided pixel and distances");
1447  return 0;
1448  }
1449  }
1450 
1451  /* no NODATA, exclude is FALSE */
1452  if (!band->hasnodata)
1453  exclude_nodata_value = FALSE;
1454  /* band is NODATA and excluding NODATA */
1455  else if (exclude_nodata_value && band->isnodata) {
1456  RASTER_DEBUG(4, "No nearest pixels possible as band is NODATA and excluding NODATA values");
1457  return 0;
1458  }
1459 
1460  /* determine the maximum distance to prevent an infinite loop */
1461  if (d0) {
1462  int a, b;
1463 
1464  /* X axis */
1465  a = abs(x);
1466  b = abs(x - band->width);
1467 
1468  if (a > b)
1469  distance[0] = a;
1470  else
1471  distance[0] = b;
1472 
1473  /* Y axis */
1474  a = abs(y);
1475  b = abs(y - band->height);
1476  if (a > b)
1477  distance[1] = a;
1478  else
1479  distance[1] = b;
1480 
1481  RASTER_DEBUGF(4, "Maximum distances: %d x %d", distance[0], distance[1]);
1482  }
1483 
1484  /* minimum possible value for pixel type */
1485  minval = rt_pixtype_get_min_value(band->pixtype);
1486  RASTER_DEBUGF(4, "pixtype: %s", rt_pixtype_name(band->pixtype));
1487  RASTER_DEBUGF(4, "minval: %f", minval);
1488 
1489  /* set variables */
1490  count = 0;
1491  *npixels = NULL;
1492 
1493  /* maximum extent */
1494  max_extent[0] = x - distance[0]; /* min X */
1495  max_extent[1] = y - distance[1]; /* min Y */
1496  max_extent[2] = x + distance[0]; /* max X */
1497  max_extent[3] = y + distance[1]; /* max Y */
1498  RASTER_DEBUGF(4, "Maximum Extent: (%d, %d, %d, %d)",
1499  max_extent[0], max_extent[1], max_extent[2], max_extent[3]);
1500 
1501  _d[0] = 0;
1502  _d[1] = 0;
1503  do {
1504  _d[0]++;
1505  _d[1]++;
1506 
1507  extent[0] = x - _d[0]; /* min x */
1508  extent[1] = y - _d[1]; /* min y */
1509  extent[2] = x + _d[0]; /* max x */
1510  extent[3] = y + _d[1]; /* max y */
1511 
1512  RASTER_DEBUGF(4, "Processing distances: %d x %d", _d[0], _d[1]);
1513  RASTER_DEBUGF(4, "Extent: (%d, %d, %d, %d)",
1514  extent[0], extent[1], extent[2], extent[3]);
1515 
1516  for (i = 0; i < 2; i++) {
1517 
1518  /* by row */
1519  if (i < 1)
1520  _max = extent[2] - extent[0] + 1;
1521  /* by column */
1522  else
1523  _max = extent[3] - extent[1] + 1;
1524  _max = abs(_max);
1525 
1526  for (j = 0; j < 2; j++) {
1527  /* by row */
1528  if (i < 1) {
1529  _x = extent[0];
1530  _min = &_x;
1531 
1532  /* top row */
1533  if (j < 1)
1534  _y = extent[1];
1535  /* bottom row */
1536  else
1537  _y = extent[3];
1538  }
1539  /* by column */
1540  else {
1541  _y = extent[1] + 1;
1542  _min = &_y;
1543 
1544  /* left column */
1545  if (j < 1) {
1546  _x = extent[0];
1547  _max -= 2;
1548  }
1549  /* right column */
1550  else
1551  _x = extent[2];
1552  }
1553 
1554  RASTER_DEBUGF(4, "_min, _max: %d, %d", *_min, _max);
1555  for (k = 0; k < (uint32_t) _max; k++) {
1556  /* check that _x and _y are not outside max extent */
1557  if (
1558  _x < max_extent[0] || _x > max_extent[2] ||
1559  _y < max_extent[1] || _y > max_extent[3]
1560  ) {
1561  (*_min)++;
1562  continue;
1563  }
1564 
1565  /* outside band extent, set to NODATA */
1566  if (
1567  (_x < 0 || _x >= band->width) ||
1568  (_y < 0 || _y >= band->height)
1569  ) {
1570  /* no NODATA, set to minimum possible value */
1571  if (!band->hasnodata)
1572  pixval = minval;
1573  /* has NODATA, use NODATA */
1574  else
1575  pixval = band->nodataval;
1576  RASTER_DEBUGF(4, "NODATA pixel outside band extent: (x, y, val) = (%d, %d, %f)", _x, _y, pixval);
1577  inextent = 0;
1578  isnodata = 1;
1579  }
1580  else {
1581  if (rt_band_get_pixel(
1582  band,
1583  _x, _y,
1584  &pixval,
1585  &isnodata
1586  ) != ES_NONE) {
1587  rterror("rt_band_get_nearest_pixel: Could not get pixel value");
1588  if (count) rtdealloc(*npixels);
1589  return -1;
1590  }
1591  RASTER_DEBUGF(4, "Pixel: (x, y, val) = (%d, %d, %f)", _x, _y, pixval);
1592  inextent = 1;
1593  }
1594 
1595  /* use pixval? */
1596  if (!exclude_nodata_value || (exclude_nodata_value && !isnodata)) {
1597  /* add pixel to result set */
1598  RASTER_DEBUGF(4, "Adding pixel to set of nearest pixels: (x, y, val) = (%d, %d, %f)", _x, _y, pixval);
1599  count++;
1600 
1601  if (*npixels == NULL)
1602  *npixels = (rt_pixel) rtalloc(sizeof(struct rt_pixel_t) * count);
1603  else
1604  *npixels = (rt_pixel) rtrealloc(*npixels, sizeof(struct rt_pixel_t) * count);
1605  if (*npixels == NULL) {
1606  rterror("rt_band_get_nearest_pixel: Could not allocate memory for nearest pixel(s)");
1607  return -1;
1608  }
1609 
1610  npixel = &((*npixels)[count - 1]);
1611  npixel->x = _x;
1612  npixel->y = _y;
1613  npixel->value = pixval;
1614 
1615  /* special case for when outside band extent */
1616  if (!inextent && !band->hasnodata)
1617  npixel->nodata = 1;
1618  else
1619  npixel->nodata = 0;
1620  }
1621 
1622  (*_min)++;
1623  }
1624  }
1625  }
1626 
1627  /* distance threshholds met */
1628  if (_d[0] >= distance[0] && _d[1] >= distance[1])
1629  break;
1630  else if (d0 && count)
1631  break;
1632  }
1633  while (1);
1634 
1635  RASTER_DEBUGF(3, "Nearest pixels in return: %d", count);
1636 
1637  return count;
1638 }
1639 
1651 int
1653  rt_band band, int exclude_nodata_value,
1654  double *searchset, int searchcount,
1655  rt_pixel *pixels
1656 ) {
1657  int x;
1658  int y;
1659  int i;
1660  double pixval;
1661  int err;
1662  int count = 0;
1663  int isnodata = 0;
1664  int isequal = 0;
1665 
1666  rt_pixel pixel = NULL;
1667 
1668  assert(NULL != band);
1669  assert(NULL != pixels);
1670  assert(NULL != searchset && searchcount > 0);
1671 
1672  if (!band->hasnodata)
1673  exclude_nodata_value = FALSE;
1674  /* band is NODATA and exclude_nodata_value = TRUE, nothing to search */
1675  else if (exclude_nodata_value && band->isnodata) {
1676  RASTER_DEBUG(4, "Pixels cannot be searched as band is NODATA and excluding NODATA values");
1677  return 0;
1678  }
1679 
1680  for (x = 0; x < band->width; x++) {
1681  for (y = 0; y < band->height; y++) {
1682  err = rt_band_get_pixel(band, x, y, &pixval, &isnodata);
1683  if (err != ES_NONE) {
1684  rterror("rt_band_get_pixel_of_value: Cannot get band pixel");
1685  return -1;
1686  }
1687  else if (exclude_nodata_value && isnodata)
1688  continue;
1689 
1690  for (i = 0; i < searchcount; i++) {
1691  if (rt_pixtype_compare_clamped_values(band->pixtype, searchset[i], pixval, &isequal) != ES_NONE) {
1692  continue;
1693  }
1694 
1695  if (FLT_NEQ(pixval, searchset[i]) || !isequal)
1696  continue;
1697 
1698  /* match found */
1699  count++;
1700  if (*pixels == NULL)
1701  *pixels = (rt_pixel) rtalloc(sizeof(struct rt_pixel_t) * count);
1702  else
1703  *pixels = (rt_pixel) rtrealloc(*pixels, sizeof(struct rt_pixel_t) * count);
1704  if (*pixels == NULL) {
1705  rterror("rt_band_get_pixel_of_value: Could not allocate memory for pixel(s)");
1706  return -1;
1707  }
1708 
1709  pixel = &((*pixels)[count - 1]);
1710  pixel->x = x;
1711  pixel->y = y;
1712  pixel->nodata = 0;
1713  pixel->value = pixval;
1714  }
1715  }
1716  }
1717 
1718  return count;
1719 }
1720 
1730 rt_band_get_nodata(rt_band band, double *nodata) {
1731  assert(NULL != band);
1732  assert(NULL != nodata);
1733 
1734  *nodata = band->nodataval;
1735 
1736  if (!band->hasnodata) {
1737  rterror("rt_band_get_nodata: Band has no NODATA value");
1738  return ES_ERROR;
1739  }
1740 
1741  return ES_NONE;
1742 }
1743 
1744 double
1746  assert(NULL != band);
1747 
1748  return rt_pixtype_get_min_value(band->pixtype);
1749 }
1750 
1751 int
1753  int i, j, err;
1754  double pxValue;
1755  int isnodata = 0;
1756 
1757  assert(NULL != band);
1758 
1759  /* Check if band has nodata value */
1760  if (!band->hasnodata) {
1761  RASTER_DEBUG(3, "Band has no NODATA value");
1762  band->isnodata = FALSE;
1763  return FALSE;
1764  }
1765 
1766  pxValue = band->nodataval;
1767 
1768  /* Check all pixels */
1769  for (i = 0; i < band->width; i++) {
1770  for (j = 0; j < band->height; j++) {
1771  err = rt_band_get_pixel(band, i, j, &pxValue, &isnodata);
1772  if (err != ES_NONE) {
1773  rterror("rt_band_check_is_nodata: Cannot get band pixel");
1774  return FALSE;
1775  }
1776  else if (!isnodata) {
1777  band->isnodata = FALSE;
1778  return FALSE;
1779  }
1780  }
1781  }
1782 
1783  band->isnodata = TRUE;
1784  return TRUE;
1785 }
1786 
1797 int
1799  int isequal = 0;
1800 
1801  assert(NULL != band);
1802 
1803  /* no NODATA, so never equal */
1804  if (!band->hasnodata)
1805  return 0;
1806 
1807  /* value is exactly NODATA */
1808  if (FLT_EQ(val, band->nodataval))
1809  return 2;
1810 
1811  /* ignore error from rt_pixtype_compare_clamped_values */
1813  band->pixtype,
1814  val, band->nodataval,
1815  &isequal
1816  );
1817 
1818  return isequal ? 1 : 0;
1819 }
1820 
1835  rt_band band,
1836  double val,
1837  double *newval, int *corrected
1838 ) {
1839  double minval = 0.;
1840 
1841  assert(NULL != band);
1842  assert(NULL != newval);
1843 
1844  if (corrected != NULL)
1845  *corrected = 0;
1846 
1847  /* no need to correct if clamped values IS NOT clamped NODATA */
1848  if (rt_band_clamped_value_is_nodata(band, val) != 1) {
1849  *newval = val;
1850  return ES_NONE;
1851  }
1852 
1853  minval = rt_pixtype_get_min_value(band->pixtype);
1854  *newval = val;
1855 
1856  switch (band->pixtype) {
1857  case PT_1BB:
1858  *newval = !band->nodataval;
1859  break;
1860  case PT_2BUI:
1861  if (rt_util_clamp_to_2BUI(val) == rt_util_clamp_to_2BUI(minval))
1862  (*newval)++;
1863  else
1864  (*newval)--;
1865  break;
1866  case PT_4BUI:
1867  if (rt_util_clamp_to_4BUI(val) == rt_util_clamp_to_4BUI(minval))
1868  (*newval)++;
1869  else
1870  (*newval)--;
1871  break;
1872  case PT_8BSI:
1873  if (rt_util_clamp_to_8BSI(val) == rt_util_clamp_to_8BSI(minval))
1874  (*newval)++;
1875  else
1876  (*newval)--;
1877  break;
1878  case PT_8BUI:
1879  if (rt_util_clamp_to_8BUI(val) == rt_util_clamp_to_8BUI(minval))
1880  (*newval)++;
1881  else
1882  (*newval)--;
1883  break;
1884  case PT_16BSI:
1885  if (rt_util_clamp_to_16BSI(val) == rt_util_clamp_to_16BSI(minval))
1886  (*newval)++;
1887  else
1888  (*newval)--;
1889  break;
1890  case PT_16BUI:
1891  if (rt_util_clamp_to_16BUI(val) == rt_util_clamp_to_16BUI(minval))
1892  (*newval)++;
1893  else
1894  (*newval)--;
1895  break;
1896  case PT_32BSI:
1897  if (rt_util_clamp_to_32BSI(val) == rt_util_clamp_to_32BSI(minval))
1898  (*newval)++;
1899  else
1900  (*newval)--;
1901  break;
1902  case PT_32BUI:
1903  if (rt_util_clamp_to_32BUI(val) == rt_util_clamp_to_32BUI(minval))
1904  (*newval)++;
1905  else
1906  (*newval)--;
1907  break;
1908  case PT_32BF:
1910  *newval += FLT_EPSILON;
1911  else
1912  *newval -= FLT_EPSILON;
1913  break;
1914  case PT_64BF:
1915  break;
1916  default:
1917  rterror("rt_band_corrected_clamped_value: Unknown pixeltype %d", band->pixtype);
1918  return ES_ERROR;
1919  }
1920 
1921  if (corrected != NULL)
1922  *corrected = 1;
1923 
1924  return ES_NONE;
1925 }
1926 
#define TRUE
Definition: dbfopen.c:169
#define FALSE
Definition: dbfopen.c:168
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
rt_pixtype rt_util_gdal_datatype_to_pixtype(GDALDataType gdt)
Convert GDALDataType to rt_pixtype.
Definition: rt_util.c:153
#define FLT_NEQ(x, y)
Definition: librtcore.h:2233
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:295
int rt_util_gdal_register_all(int force_register_all)
Definition: rt_util.c:336
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:299
void rt_raster_set_geotransform_matrix(rt_raster raster, double *gt)
Set raster's geotransform using 6-element array.
Definition: rt_raster.c:727
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
int32_t rt_util_clamp_to_32BSI(double value)
Definition: rt_util.c:69
rt_errorstate rt_raster_geopoint_to_cell(rt_raster raster, double xw, double yw, double *xr, double *yr, double *igt)
Convert an xw, yw map point to a xr, yr raster point.
Definition: rt_raster.c:806
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
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
GDALDataType rt_util_pixtype_to_gdal_datatype(rt_pixtype pt)
Convert rt_pixtype to GDALDataType.
Definition: rt_util.c:118
double rt_pixtype_get_min_value(rt_pixtype pixtype)
Return minimum value possible for pixel type.
Definition: rt_pixel.c:148
struct rt_pixel_t * rt_pixel
Definition: librtcore.h:147
#define FLT_EQ(x, y)
Definition: librtcore.h:2234
rt_raster rt_raster_from_gdal_dataset(GDALDatasetH ds)
Return a raster from a GDAL dataset.
Definition: rt_raster.c:2182
uint8_t rt_util_clamp_to_2BUI(double value)
Definition: rt_util.c:39
void rtwarn(const char *fmt,...)
Definition: rt_context.c:224
rt_errorstate rt_pixtype_compare_clamped_values(rt_pixtype pixtype, double val, double refval, int *isequal)
Test to see if two values are equal when clamped.
Definition: rt_pixel.c:200
uint8_t rt_util_clamp_to_8BUI(double value)
Definition: rt_util.c:54
GDALDatasetH rt_util_gdal_open(const char *fn, GDALAccess fn_access, int shared)
Definition: rt_util.c:381
rt_errorstate
Enum definitions.
Definition: librtcore.h:179
@ ES_NONE
Definition: librtcore.h:180
@ ES_ERROR
Definition: librtcore.h:181
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
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster's SRID.
Definition: rt_raster.c:363
uint8_t rt_util_clamp_to_4BUI(double value)
Definition: rt_util.c:44
const char * rt_pixtype_name(rt_pixtype pixtype)
Definition: rt_pixel.c:110
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
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)
int rt_pixtype_size(rt_pixtype pixtype)
Return size in bytes of a value in the given pixtype.
Definition: rt_pixel.c:39
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.
Datum distance(PG_FUNCTION_ARGS)
int value
Definition: genraster.py:61
int count
Definition: genraster.py:56
band
Definition: ovdump.py:57
data
Definition: ovdump.py:103
pixval
Definition: pixval.py:93
pixel
Definition: pixval.py:90
nband
Definition: pixval.py:52
Definition: pixval.py:1
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
int rt_band_get_pixel_of_value(rt_band band, int exclude_nodata_value, double *searchset, int searchcount, rt_pixel *pixels)
Search band for pixel(s) with search values.
Definition: rt_band.c:1652
rt_band rt_band_duplicate(rt_band band)
Create a new band duplicated from source band.
Definition: rt_band.c:287
uint16_t rt_band_get_width(rt_band band)
Return width of this band.
Definition: rt_band.c:640
void rt_band_set_hasnodata_flag(rt_band band, int flag)
Set hasnodata flag value.
Definition: rt_band.c:681
rt_errorstate rt_band_set_isnodata_flag(rt_band band, int flag)
Set isnodata flag value.
Definition: rt_band.c:695
const char * rt_band_get_ext_path(rt_band band)
Return band's external path (only valid when rt_band_is_offline returns non-zero).
Definition: rt_band.c:363
rt_errorstate rt_band_load_offline_data(rt_band band)
Load offline band's data.
Definition: rt_band.c:429
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
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition: rt_band.c:714
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
uint64_t rt_band_get_file_size(rt_band band)
Return file size in bytes.
Definition: rt_band.c:586
int rt_band_check_is_nodata(rt_band band)
Returns TRUE if the band is only nodata values.
Definition: rt_band.c:1752
rt_errorstate rt_band_set_nodata(rt_band band, double val, int *converted)
Set nodata value.
Definition: rt_band.c:733
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 rt_band_set_pixel_line(rt_band band, int x, int y, void *vals, uint32_t len)
Set values of multiple pixels.
Definition: rt_band.c:853
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition: rt_band.c:340
rt_errorstate rt_band_corrected_clamped_value(rt_band band, double val, double *newval, int *corrected)
Correct value when clamped value is equal to clamped NODATA value.
Definition: rt_band.c:1834
rt_band rt_band_new_offline(uint16_t width, uint16_t height, rt_pixtype pixtype, uint32_t hasnodata, double nodataval, uint8_t bandNum, const char *path)
Create an out-db rt_band.
Definition: rt_band.c:124
rt_errorstate rt_band_get_ext_band_num(rt_band band, uint8_t *bandnum)
Return bands' external band number (only valid when rt_band_is_offline returns non-zero).
Definition: rt_band.c:376
uint64_t rt_band_get_file_timestamp(rt_band band)
Return file timestamp.
Definition: rt_band.c:608
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1730
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition: rt_band.c:631
rt_band rt_band_new_offline_from_path(uint16_t width, uint16_t height, int hasnodata, double nodataval, uint8_t bandNum, const char *path, int force)
Create an out-db rt_band from path.
Definition: rt_band.c:199
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 * rt_band_get_data(rt_band band)
Get pointer to raster band data.
Definition: rt_band.c:400
int rt_band_clamped_value_is_nodata(rt_band band, double val)
Compare clamped value to band's clamped NODATA value.
Definition: rt_band.c:1798
char enable_outdb_rasters
Definition: rt_band.c:417
int rt_band_get_ownsdata_flag(rt_band band)
Return 0 (FALSE) or non-zero (TRUE) indicating if rt_band is responsible for managing the memory for ...
Definition: rt_band.c:659
int rt_band_is_offline(rt_band band)
Return non-zero if the given band data is on the filesystem.
Definition: rt_band.c:329
uint16_t rt_band_get_height(rt_band band)
Return height of this band.
Definition: rt_band.c:649
rt_errorstate rt_band_get_pixel_line(rt_band band, int x, int y, uint16_t len, void **vals, uint16_t *nvals)
Get values of multiple pixels.
Definition: rt_band.c:1137
union rt_band_t::@8 data
void * mem
Definition: librtcore.h:2327
double value
Definition: librtcore.h:2338
uint8_t nodata
Definition: librtcore.h:2337
unsigned int uint32_t
Definition: uthash.h:78
unsigned char uint8_t
Definition: uthash.h:79