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
42 /**
43 * Create an in-db rt_band with no data
44 *
45 * @param width : number of pixel columns
46 * @param height : number of pixel rows
47 * @param pixtype : pixel type for the band
48 * @param hasnodata : indicates if the band has nodata value
49 * @param nodataval : the nodata value, will be appropriately
50 * truncated to fit the pixtype size.
51 * @param data : pointer to actual band data, required to
52 * be aligned accordingly to
53 * rt_pixtype_aligment(pixtype) and big enough
54 * to hold raster width*height values.
55 * Data will NOT be copied, ownership is left
56 * to caller which is responsible to keep it
57 * allocated for the whole lifetime of the returned
58 * rt_band.
59 *
60 * @return an rt_band, or 0 on failure
61 */
62 rt_band
rt_band_new_inline(uint16_t width,uint16_t height,rt_pixtype pixtype,uint32_t hasnodata,double nodataval,uint8_t * data)63 rt_band_new_inline(
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");
97 rt_band_destroy(band);
98 return NULL;
99 }
100
101 return band;
102 }
103
104 /**
105 * Create an out-db rt_band
106 *
107 * @param width : number of pixel columns
108 * @param height : number of pixel rows
109 * @param pixtype : pixel type for the band
110 * @param hasnodata : indicates if the band has nodata value
111 * @param nodataval : the nodata value, will be appropriately
112 * truncated to fit the pixtype size.
113 * @param bandNum : 0-based band number in the external file
114 * to associate this band with.
115 * @param path : NULL-terminated path string pointing to the file
116 * containing band data. The string will NOT be
117 * copied, ownership is left to caller which is
118 * responsible to keep it allocated for the whole
119 * lifetime of the returned rt_band.
120 *
121 * @return an rt_band, or 0 on failure
122 */
123 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)124 rt_band_new_offline(
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");
158 rt_band_destroy(band);
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");
169 rt_band_destroy(band);
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
180 /**
181 * Create an out-db rt_band from path
182 *
183 * @param width : number of pixel columns
184 * @param height : number of pixel rows
185 * @param hasnodata : indicates if the band has nodata value
186 * @param nodataval : the nodata value, will be appropriately
187 * @param bandNum : 1-based band number in the external file
188 * to associate this band with.
189 * @param path : NULL-terminated path string pointing to the file
190 * containing band data. The string will NOT be
191 * copied, ownership is left to caller which is
192 * responsible to keep it allocated for the whole
193 * lifetime of the returned rt_band.
194 * @param force : if True, ignore all validation checks
195 *
196 * @return an rt_band, or 0 on failure
197 */
198 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)199 rt_band_new_offline_from_path(
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 */
216 rt_util_gdal_register_all(0);
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
276 /**
277 * Create a new band duplicated from source band. Memory is allocated
278 * for band path (if band is offline) or band data (if band is online).
279 * The caller is responsible for freeing the memory when the returned
280 * rt_band is destroyed.
281 *
282 * @param : the band to copy
283 *
284 * @return an rt_band or NULL on failure
285 */
286 rt_band
rt_band_duplicate(rt_band band)287 rt_band_duplicate(rt_band 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
rt_band_is_offline(rt_band band)329 rt_band_is_offline(rt_band band) {
330 assert(NULL != band);
331 return band->offline ? 1 : 0;
332 }
333
334 /**
335 * Destroy a raster band
336 *
337 * @param band : the band to destroy
338 */
339 void
rt_band_destroy(rt_band band)340 rt_band_destroy(rt_band band) {
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*
rt_band_get_ext_path(rt_band band)363 rt_band_get_ext_path(rt_band band) {
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
375 rt_errorstate
rt_band_get_ext_band_num(rt_band band,uint8_t * bandnum)376 rt_band_get_ext_band_num(rt_band band, uint8_t *bandnum) {
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
392 /**
393 * Get pointer to raster band data
394 *
395 * @param band : the band who's data to get
396 *
397 * @return pointer to band data or NULL if error
398 */
399 void *
rt_band_get_data(rt_band band)400 rt_band_get_data(rt_band band) {
401 assert(NULL != band);
402
403 if (band->offline) {
404 if (band->data.offline.mem != NULL)
405 return band->data.offline.mem;
406
407 if (rt_band_load_offline_data(band) != ES_NONE)
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 */
417 char enable_outdb_rasters = 1;
418
419 /**
420 * Load offline band's data. Loaded data is internally owned
421 * and should not be released by the caller. Data will be
422 * released when band is destroyed with rt_band_destroy().
423 *
424 * @param band : the band who's data to get
425 *
426 * @return ES_NONE if success, ES_ERROR if failure
427 */
428 rt_errorstate
rt_band_load_offline_data(rt_band band)429 rt_band_load_offline_data(rt_band band) {
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
460 rt_util_gdal_register_all(0);
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);
496 rt_raster_set_geotransform_matrix(_rast, ogt);
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 */
511 rt_raster_geopoint_to_cell(
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
rt_band_get_file_size(rt_band band)586 uint64_t rt_band_get_file_size(rt_band band) {
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
rt_band_get_file_timestamp(rt_band band)608 uint64_t rt_band_get_file_timestamp(rt_band band) {
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
630 rt_pixtype
rt_band_get_pixtype(rt_band band)631 rt_band_get_pixtype(rt_band band) {
632
633 assert(NULL != band);
634
635
636 return band->pixtype;
637 }
638
639 uint16_t
rt_band_get_width(rt_band band)640 rt_band_get_width(rt_band band) {
641
642 assert(NULL != band);
643
644
645 return band->width;
646 }
647
648 uint16_t
rt_band_get_height(rt_band band)649 rt_band_get_height(rt_band band) {
650
651 assert(NULL != band);
652
653
654 return band->height;
655 }
656
657 /* Get ownsdata flag */
658 int
rt_band_get_ownsdata_flag(rt_band band)659 rt_band_get_ownsdata_flag(rt_band band) {
660 assert(NULL != band);
661
662 return band->ownsdata ? 1 : 0;
663 }
664
665 /* set ownsdata flag */
666 void
rt_band_set_ownsdata_flag(rt_band band,int flag)667 rt_band_set_ownsdata_flag(rt_band band, int flag) {
668 assert(NULL != band);
669
670 band->ownsdata = flag ? 1 : 0;
671 }
672
673 int
rt_band_get_hasnodata_flag(rt_band band)674 rt_band_get_hasnodata_flag(rt_band band) {
675 assert(NULL != band);
676
677 return band->hasnodata ? 1 : 0;
678 }
679
680 void
rt_band_set_hasnodata_flag(rt_band band,int flag)681 rt_band_set_hasnodata_flag(rt_band band, int flag) {
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
694 rt_errorstate
rt_band_set_isnodata_flag(rt_band band,int flag)695 rt_band_set_isnodata_flag(rt_band band, int flag) {
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
rt_band_get_isnodata_flag(rt_band band)714 rt_band_get_isnodata_flag(rt_band band) {
715 assert(NULL != band);
716
717 if (band->hasnodata)
718 return band->isnodata ? 1 : 0;
719 else
720 return 0;
721 }
722
723 /**
724 * Set nodata value
725 *
726 * @param band : the band to set nodata value to
727 * @param val : the nodata value
728 * @param converted : if non-zero, value was truncated/clamped/coverted
729 *
730 * @return ES_NONE or ES_ERROR
731 */
732 rt_errorstate
rt_band_set_nodata(rt_band band,double val,int * converted)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
821 if (rt_util_dbl_trunc_warning(
822 val,
823 checkvalint, checkvaluint,
824 checkvalfloat, checkvaldouble,
825 pixtype
826 ) && converted != NULL) {
827 *converted = 1;
828 }
829
830 return ES_NONE;
831 }
832
833 /**
834 * Set values of multiple pixels. Unlike rt_band_set_pixel,
835 * values in vals are expected to be of the band's pixel type
836 * as this function uses memcpy.
837 *
838 * It is important to be careful when using this function as
839 * the number of values being set may exceed a pixel "row".
840 * Remember that the band values are stored in a stream (1-D array)
841 * regardless of what the raster's width and height might be.
842 * So, setting a number of values may cross multiple pixel "rows".
843 *
844 * @param band : the band to set value to
845 * @param x : X coordinate (0-based)
846 * @param y : Y coordinate (0-based)
847 * @param vals : the pixel values to apply
848 * @param len : # of elements in vals
849 *
850 * @return ES_NONE on success, ES_ERROR on error
851 */
852 rt_errorstate
rt_band_set_pixel_line(rt_band band,int x,int y,void * vals,uint32_t len)853 rt_band_set_pixel_line(
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
884 data = rt_band_get_data(band);
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 */
956 if (rt_band_get_hasnodata_flag(band))
957 rt_band_set_isnodata_flag(band, 0);
958
959 return ES_NONE;
960 }
961
962 /**
963 * Set single pixel's value
964 *
965 * @param band : the band to set value to
966 * @param x : x ordinate (0-based)
967 * @param y : y ordinate (0-based)
968 * @param val : the pixel value
969 * @param converted : (optional) non-zero if value truncated/clamped/converted
970 *
971 * @return ES_NONE on success, ES_ERROR on error
972 */
973 rt_errorstate
rt_band_set_pixel(rt_band band,int x,int y,double val,int * converted)974 rt_band_set_pixel(
975 rt_band band,
976 int x, int y,
977 double val,
978 int *converted
979 ) {
980 rt_pixtype pixtype = PT_END;
981 uint8_t *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
1027 data = rt_band_get_data(band);
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] = (uint8_t)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 */
1099 if (!rt_band_clamped_value_is_nodata(band, val)) {
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 */
1105 if (rt_util_dbl_trunc_warning(
1106 val,
1107 checkvalint, checkvaluint,
1108 checkvalfloat, checkvaldouble,
1109 pixtype
1110 ) && converted != NULL) {
1111 *converted = 1;
1112 }
1113
1114 return ES_NONE;
1115 }
1116
1117 /**
1118 * Get values of multiple pixels. Unlike rt_band_get_pixel,
1119 * values in vals are of the band's pixel type so cannot be
1120 * assumed to be double. Function uses memcpy.
1121 *
1122 * It is important to be careful when using this function as
1123 * the number of values being fetched may exceed a pixel "row".
1124 * Remember that the band values are stored in a stream (1-D array)
1125 * regardless of what the raster's width and height might be.
1126 * So, getting a number of values may cross multiple pixel "rows".
1127 *
1128 * @param band : the band to get pixel value from
1129 * @param x : pixel column (0-based)
1130 * @param y : pixel row (0-based)
1131 * @param len : the number of pixels to get
1132 * @param **vals : the pixel values
1133 * @param *nvals : the number of pixel values being returned
1134 *
1135 * @return ES_NONE on success, ES_ERROR on error
1136 */
rt_band_get_pixel_line(rt_band band,int x,int y,uint16_t len,void ** vals,uint16_t * nvals)1137 rt_errorstate rt_band_get_pixel_line(
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
1168 data = rt_band_get_data(band);
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
1208 /**
1209 * Get pixel value. If band's isnodata flag is TRUE, value returned
1210 * will be the band's NODATA value
1211 *
1212 * @param band : the band to set nodata value to
1213 * @param x : x ordinate (0-based)
1214 * @param y : x ordinate (0-based)
1215 * @param *value : pixel value
1216 * @param *nodata : 0 if pixel is not NODATA
1217 *
1218 * @return 0 on success, -1 on error (value out of valid range).
1219 */
1220 rt_errorstate
rt_band_get_pixel(rt_band band,int x,int y,double * value,int * nodata)1221 rt_band_get_pixel(
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
1254 data = rt_band_get_data(band);
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 = (int8_t)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) {
1350 if (rt_band_clamped_value_is_nodata(band, *value))
1351 *nodata = 1;
1352 }
1353
1354 return ES_NONE;
1355 }
1356
1357 /**
1358 * Get nearest pixel(s) with value (not NODATA) to specified pixel
1359 *
1360 * @param band : the band to get nearest pixel(s) from
1361 * @param x : the column of the pixel (0-based)
1362 * @param y : the line of the pixel (0-based)
1363 * @param distancex : the number of pixels around the specified pixel
1364 * along the X axis
1365 * @param distancey : the number of pixels around the specified pixel
1366 * along the Y axis
1367 * @param exclude_nodata_value : if non-zero, ignore nodata values
1368 * to check for pixels with value
1369 * @param npixels : return set of rt_pixel object or NULL
1370 *
1371 * @return -1 on error, otherwise the number of rt_pixel objects
1372 * in npixels
1373 */
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)1374 uint32_t rt_band_get_nearest_pixel(
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 - (int)distance[0]; /* min X */
1495 max_extent[1] = y - (int)distance[1]; /* min Y */
1496 max_extent[2] = x + (int)distance[0]; /* max X */
1497 max_extent[3] = y + (int)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 - (int)_d[0]; /* min x */
1508 extent[1] = y - (int)_d[1]; /* min y */
1509 extent[2] = x + (int)_d[0]; /* max x */
1510 extent[3] = y + (int)_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
1640 /**
1641 * Search band for pixel(s) with search values
1642 *
1643 * @param band : the band to query for minimum and maximum pixel values
1644 * @param exclude_nodata_value : if non-zero, ignore nodata values
1645 * @param searchset : array of values to count
1646 * @param searchcount : the number of search values
1647 * @param pixels : pixels with the search value
1648 *
1649 * @return -1 on error, otherwise number of pixels
1650 */
1651 int
rt_band_get_pixel_of_value(rt_band band,int exclude_nodata_value,double * searchset,int searchcount,rt_pixel * pixels)1652 rt_band_get_pixel_of_value(
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
1721 /**
1722 * Get NODATA value
1723 *
1724 * @param band : the band whose NODATA value will be returned
1725 * @param nodata : the band's NODATA value
1726 *
1727 * @return ES_NONE or ES_ERROR
1728 */
1729 rt_errorstate
rt_band_get_nodata(rt_band band,double * nodata)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
rt_band_get_min_value(rt_band band)1745 rt_band_get_min_value(rt_band band) {
1746 assert(NULL != band);
1747
1748 return rt_pixtype_get_min_value(band->pixtype);
1749 }
1750
1751 int
rt_band_check_is_nodata(rt_band band)1752 rt_band_check_is_nodata(rt_band band) {
1753 int i, j, err;
1754 double pxValue;
1755 int isnodata = 0;
1756
1757 assert(NULL != band);
1758 band->isnodata = FALSE;
1759
1760 /* Check if band has nodata value */
1761 if (!band->hasnodata) {
1762 RASTER_DEBUG(3, "Band has no NODATA value");
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
1787 /**
1788 * Compare clamped value to band's clamped NODATA value.
1789 *
1790 * @param band : the band whose NODATA value will be used for comparison
1791 * @param val : the value to compare to the NODATA value
1792 *
1793 * @return 2 if unclamped value is unclamped NODATA
1794 * 1 if clamped value is clamped NODATA
1795 * 0 if clamped value is NOT clamped NODATA
1796 */
1797 int
rt_band_clamped_value_is_nodata(rt_band band,double val)1798 rt_band_clamped_value_is_nodata(rt_band band, double val) {
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 */
1812 rt_pixtype_compare_clamped_values(
1813 band->pixtype,
1814 val, band->nodataval,
1815 &isequal
1816 );
1817
1818 return isequal ? 1 : 0;
1819 }
1820
1821 /**
1822 * Correct value when clamped value is equal to clamped NODATA value.
1823 * Correction does NOT occur if unclamped value is exactly unclamped
1824 * NODATA value.
1825 *
1826 * @param band : the band whose NODATA value will be used for comparison
1827 * @param val : the value to compare to the NODATA value and correct
1828 * @param *newval : pointer to corrected value
1829 * @param *corrected : (optional) non-zero if val was corrected
1830 *
1831 * @return ES_NONE if success, ES_ERROR if error
1832 */
1833 rt_errorstate
rt_band_corrected_clamped_value(rt_band band,double val,double * newval,int * corrected)1834 rt_band_corrected_clamped_value(
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:
1909 if (FLT_EQ(rt_util_clamp_to_32F(val), rt_util_clamp_to_32F(minval)))
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
1927