1 /*
2 *
3 * WKTRaster - Raster Types for PostGIS
4 * http://trac.osgeo.org/postgis/wiki/WKTRaster
5 *
6 * Copyright (C) 2013 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 #include "librtcore.h"
32 #include "librtcore_internal.h"
33
34 #include <math.h>
35
36 /**
37 * Construct a raster with given dimensions.
38 *
39 * Transform will be set to identity.
40 * Will contain no bands.
41 *
42 * @param width : number of pixel columns
43 * @param height : number of pixel rows
44 *
45 * @return an rt_raster or NULL if out of memory
46 */
47 rt_raster
rt_raster_new(uint32_t width,uint32_t height)48 rt_raster_new(uint32_t width, uint32_t height) {
49 rt_raster ret = NULL;
50
51 ret = (rt_raster) rtalloc(sizeof (struct rt_raster_t));
52 if (!ret) {
53 rterror("rt_raster_new: Out of virtual memory creating an rt_raster");
54 return NULL;
55 }
56
57 RASTER_DEBUGF(3, "Created rt_raster @ %p", ret);
58
59 if (width > 65535 || height > 65535) {
60 rterror("rt_raster_new: Dimensions requested exceed the maximum (65535 x 65535) permitted for a raster");
61 rt_raster_destroy(ret);
62 return NULL;
63 }
64
65 ret->width = width;
66 ret->height = height;
67 ret->scaleX = 1;
68 ret->scaleY = -1;
69 ret->ipX = 0.0;
70 ret->ipY = 0.0;
71 ret->skewX = 0.0;
72 ret->skewY = 0.0;
73 ret->srid = SRID_UNKNOWN;
74
75 ret->numBands = 0;
76 ret->bands = NULL;
77
78 return ret;
79 }
80
81 void
rt_raster_destroy(rt_raster raster)82 rt_raster_destroy(rt_raster raster) {
83 if (raster == NULL)
84 return;
85
86 RASTER_DEBUGF(3, "Destroying rt_raster @ %p", raster);
87
88 if (raster->bands)
89 rtdealloc(raster->bands);
90
91 rtdealloc(raster);
92 }
93
94 static void
_rt_raster_geotransform_warn_offline_band(rt_raster raster)95 _rt_raster_geotransform_warn_offline_band(rt_raster raster) {
96 int numband = 0;
97 int i = 0;
98 rt_band band = NULL;
99
100 if (raster == NULL)
101 return;
102
103 numband = rt_raster_get_num_bands(raster);
104 if (numband < 1)
105 return;
106
107 for (i = 0; i < numband; i++) {
108 band = rt_raster_get_band(raster, i);
109 if (NULL == band)
110 continue;
111
112 if (!rt_band_is_offline(band))
113 continue;
114
115 rtwarn("Changes made to raster geotransform matrix may affect out-db band data. Returned band data may be incorrect");
116 break;
117 }
118 }
119
120 uint16_t
rt_raster_get_width(rt_raster raster)121 rt_raster_get_width(rt_raster raster) {
122
123 assert(NULL != raster);
124
125 return raster->width;
126 }
127
128 uint16_t
rt_raster_get_height(rt_raster raster)129 rt_raster_get_height(rt_raster raster) {
130
131 assert(NULL != raster);
132
133 return raster->height;
134 }
135
136 void
rt_raster_set_scale(rt_raster raster,double scaleX,double scaleY)137 rt_raster_set_scale(
138 rt_raster raster,
139 double scaleX, double scaleY
140 ) {
141 assert(NULL != raster);
142
143 raster->scaleX = scaleX;
144 raster->scaleY = scaleY;
145
146 _rt_raster_geotransform_warn_offline_band(raster);
147 }
148
149 double
rt_raster_get_x_scale(rt_raster raster)150 rt_raster_get_x_scale(rt_raster raster) {
151
152
153 assert(NULL != raster);
154
155 return raster->scaleX;
156 }
157
158 double
rt_raster_get_y_scale(rt_raster raster)159 rt_raster_get_y_scale(rt_raster raster) {
160
161
162 assert(NULL != raster);
163
164 return raster->scaleY;
165 }
166
167 void
rt_raster_set_skews(rt_raster raster,double skewX,double skewY)168 rt_raster_set_skews(
169 rt_raster raster,
170 double skewX, double skewY
171 ) {
172 assert(NULL != raster);
173
174 raster->skewX = skewX;
175 raster->skewY = skewY;
176
177 _rt_raster_geotransform_warn_offline_band(raster);
178 }
179
180 double
rt_raster_get_x_skew(rt_raster raster)181 rt_raster_get_x_skew(rt_raster raster) {
182
183
184 assert(NULL != raster);
185
186 return raster->skewX;
187 }
188
189 double
rt_raster_get_y_skew(rt_raster raster)190 rt_raster_get_y_skew(rt_raster raster) {
191
192
193 assert(NULL != raster);
194
195 return raster->skewY;
196 }
197
198 void
rt_raster_set_offsets(rt_raster raster,double x,double y)199 rt_raster_set_offsets(
200 rt_raster raster,
201 double x, double y
202 ) {
203
204 assert(NULL != raster);
205
206 raster->ipX = x;
207 raster->ipY = y;
208
209 _rt_raster_geotransform_warn_offline_band(raster);
210 }
211
212 double
rt_raster_get_x_offset(rt_raster raster)213 rt_raster_get_x_offset(rt_raster raster) {
214
215
216 assert(NULL != raster);
217
218 return raster->ipX;
219 }
220
221 double
rt_raster_get_y_offset(rt_raster raster)222 rt_raster_get_y_offset(rt_raster raster) {
223
224
225 assert(NULL != raster);
226
227 return raster->ipY;
228 }
229
230 void
rt_raster_get_phys_params(rt_raster rast,double * i_mag,double * j_mag,double * theta_i,double * theta_ij)231 rt_raster_get_phys_params(rt_raster rast,
232 double *i_mag, double *j_mag, double *theta_i, double *theta_ij)
233 {
234 double o11, o12, o21, o22 ; /* geotransform coefficients */
235
236 if (rast == NULL) return ;
237 if ( (i_mag==NULL) || (j_mag==NULL) || (theta_i==NULL) || (theta_ij==NULL))
238 return ;
239
240 /* retrieve coefficients from raster */
241 o11 = rt_raster_get_x_scale(rast) ;
242 o12 = rt_raster_get_x_skew(rast) ;
243 o21 = rt_raster_get_y_skew(rast) ;
244 o22 = rt_raster_get_y_scale(rast) ;
245
246 rt_raster_calc_phys_params(o11, o12, o21, o22, i_mag, j_mag, theta_i, theta_ij);
247 }
248
249 void
rt_raster_calc_phys_params(double xscale,double xskew,double yskew,double yscale,double * i_mag,double * j_mag,double * theta_i,double * theta_ij)250 rt_raster_calc_phys_params(double xscale, double xskew, double yskew, double yscale,
251 double *i_mag, double *j_mag, double *theta_i, double *theta_ij)
252
253 {
254 double theta_test ;
255
256 if ( (i_mag==NULL) || (j_mag==NULL) || (theta_i==NULL) || (theta_ij==NULL))
257 return ;
258
259 /* pixel size in the i direction */
260 *i_mag = sqrt(xscale*xscale + yskew*yskew) ;
261
262 /* pixel size in the j direction */
263 *j_mag = sqrt(xskew*xskew + yscale*yscale) ;
264
265 /* Rotation
266 * ========
267 * Two steps:
268 * 1] calculate the magnitude of the angle between the x axis and
269 * the i basis vector.
270 * 2] Calculate the sign of theta_i based on the angle between the y axis
271 * and the i basis vector.
272 */
273 *theta_i = acos(xscale/(*i_mag)) ; /* magnitude */
274 theta_test = acos(yskew/(*i_mag)) ; /* sign */
275 if (theta_test < M_PI_2){
276 *theta_i = -(*theta_i) ;
277 }
278
279
280 /* Angular separation of basis vectors
281 * ===================================
282 * Two steps:
283 * 1] calculate the magnitude of the angle between the j basis vector and
284 * the i basis vector.
285 * 2] Calculate the sign of theta_ij based on the angle between the
286 * perpendicular of the i basis vector and the j basis vector.
287 */
288 *theta_ij = acos(((xscale*xskew) + (yskew*yscale))/((*i_mag)*(*j_mag))) ;
289 theta_test = acos( ((-yskew*xskew)+(xscale*yscale)) /
290 ((*i_mag)*(*j_mag)));
291 if (theta_test > M_PI_2) {
292 *theta_ij = -(*theta_ij) ;
293 }
294 }
295
296 void
rt_raster_set_phys_params(rt_raster rast,double i_mag,double j_mag,double theta_i,double theta_ij)297 rt_raster_set_phys_params(rt_raster rast,double i_mag, double j_mag, double theta_i, double theta_ij)
298 {
299 double o11, o12, o21, o22 ; /* calculated geotransform coefficients */
300 int success ;
301
302 if (rast == NULL) return ;
303
304 success = rt_raster_calc_gt_coeff(i_mag, j_mag, theta_i, theta_ij,
305 &o11, &o12, &o21, &o22) ;
306
307 if (success) {
308 rt_raster_set_scale(rast, o11, o22) ;
309 rt_raster_set_skews(rast, o12, o21) ;
310 }
311 }
312
313 int
rt_raster_calc_gt_coeff(double i_mag,double j_mag,double theta_i,double theta_ij,double * xscale,double * xskew,double * yskew,double * yscale)314 rt_raster_calc_gt_coeff(double i_mag, double j_mag, double theta_i, double theta_ij,
315 double *xscale, double *xskew, double *yskew, double *yscale)
316 {
317 double f ; /* reflection flag 1.0 or -1.0 */
318 double k_i ; /* shearing coefficient */
319 double s_i, s_j ; /* scaling coefficients */
320 double cos_theta_i, sin_theta_i ;
321
322 if ( (xscale==NULL) || (xskew==NULL) || (yskew==NULL) || (yscale==NULL)) {
323 return 0;
324 }
325
326 if ( (theta_ij == 0.0) || (theta_ij == M_PI)) {
327 return 0;
328 }
329
330 /* Reflection across the i axis */
331 f=1.0 ;
332 if (theta_ij < 0) {
333 f = -1.0;
334 }
335
336 /* scaling along i axis */
337 s_i = i_mag ;
338
339 /* shearing parallel to i axis */
340 k_i = tan(f*M_PI_2 - theta_ij) ;
341
342 /* scaling along j axis */
343 s_j = j_mag / (sqrt(k_i*k_i + 1)) ;
344
345 /* putting it altogether */
346 cos_theta_i = cos(theta_i) ;
347 sin_theta_i = sin(theta_i) ;
348 *xscale = s_i * cos_theta_i ;
349 *xskew = k_i * s_j * f * cos_theta_i + s_j * f * sin_theta_i ;
350 *yskew = -s_i * sin_theta_i ;
351 *yscale = -k_i * s_j * f * sin_theta_i + s_j * f * cos_theta_i ;
352 return 1;
353 }
354
355 int32_t
rt_raster_get_srid(rt_raster raster)356 rt_raster_get_srid(rt_raster raster) {
357 assert(NULL != raster);
358
359 return clamp_srid(raster->srid);
360 }
361
362 void
rt_raster_set_srid(rt_raster raster,int32_t srid)363 rt_raster_set_srid(rt_raster raster, int32_t srid) {
364 assert(NULL != raster);
365
366 raster->srid = clamp_srid(srid);
367
368 _rt_raster_geotransform_warn_offline_band(raster);
369 }
370
371 uint16_t
rt_raster_get_num_bands(rt_raster raster)372 rt_raster_get_num_bands(rt_raster raster) {
373
374
375 assert(NULL != raster);
376
377 return raster->numBands;
378 }
379
380 rt_band
rt_raster_get_band(rt_raster raster,int n)381 rt_raster_get_band(rt_raster raster, int n) {
382 assert(NULL != raster);
383
384 if (n >= raster->numBands || n < 0)
385 return NULL;
386
387 return raster->bands[n];
388 }
389
390 /******************************************************************************
391 * rt_raster_add_band()
392 ******************************************************************************/
393
394 /**
395 * Add band data to a raster.
396 *
397 * @param raster : the raster to add a band to
398 * @param band : the band to add, ownership left to caller.
399 * Band dimensions are required to match with raster ones.
400 * @param index : the position where to insert the new band (0 based)
401 *
402 * @return identifier (position) for the just-added raster, or -1 on error
403 */
404 int
rt_raster_add_band(rt_raster raster,rt_band band,int index)405 rt_raster_add_band(rt_raster raster, rt_band band, int index) {
406 rt_band *oldbands = NULL;
407 rt_band oldband = NULL;
408 rt_band tmpband = NULL;
409 uint16_t i = 0;
410
411 assert(NULL != raster);
412 assert(NULL != band);
413
414 RASTER_DEBUGF(3, "Adding band %p to raster %p", band, raster);
415
416 if (band->width != raster->width || band->height != raster->height) {
417 rterror("rt_raster_add_band: Can't add a %dx%d band to a %dx%d raster",
418 band->width, band->height, raster->width, raster->height);
419 return -1;
420 }
421
422 if (index > raster->numBands)
423 index = raster->numBands;
424
425 if (index < 0)
426 index = 0;
427
428 oldbands = raster->bands;
429
430 RASTER_DEBUGF(3, "Oldbands at %p", oldbands);
431
432 raster->bands = (rt_band*) rtrealloc(raster->bands,
433 sizeof (rt_band)*(raster->numBands + 1)
434 );
435
436 RASTER_DEBUG(3, "Checking bands");
437
438 if (NULL == raster->bands) {
439 rterror("rt_raster_add_band: Out of virtual memory "
440 "reallocating band pointers");
441 raster->bands = oldbands;
442 return -1;
443 }
444
445 RASTER_DEBUGF(4, "realloc returned %p", raster->bands);
446
447 for (i = 0; i <= raster->numBands; ++i) {
448 if (i == index) {
449 oldband = raster->bands[i];
450 raster->bands[i] = band;
451 } else if (i > index) {
452 tmpband = raster->bands[i];
453 raster->bands[i] = oldband;
454 oldband = tmpband;
455 }
456 }
457
458 band->raster = raster;
459
460 raster->numBands++;
461
462 RASTER_DEBUGF(4, "Raster now has %d bands", raster->numBands);
463
464 return index;
465 }
466
467 /******************************************************************************
468 * rt_raster_generate_new_band()
469 ******************************************************************************/
470
471 /**
472 * Generate a new inline band and add it to a raster.
473 * Memory is allocated in this function for band data.
474 *
475 * @param raster : the raster to add a band to
476 * @param pixtype : the pixel type for the new band
477 * @param initialvalue : initial value for pixels
478 * @param hasnodata : indicates if the band has a nodata value
479 * @param nodatavalue : nodata value for the new band
480 * @param index : position to add the new band in the raster
481 *
482 * @return identifier (position) for the just-added raster, or -1 on error
483 */
484 int
rt_raster_generate_new_band(rt_raster raster,rt_pixtype pixtype,double initialvalue,uint32_t hasnodata,double nodatavalue,int index)485 rt_raster_generate_new_band(
486 rt_raster raster, rt_pixtype pixtype,
487 double initialvalue, uint32_t hasnodata, double nodatavalue,
488 int index
489 ) {
490 rt_band band = NULL;
491 int width = 0;
492 int height = 0;
493 int numval = 0;
494 int datasize = 0;
495 int oldnumbands = 0;
496 int numbands = 0;
497 void * mem = NULL;
498 int32_t checkvalint = 0;
499 uint32_t checkvaluint = 0;
500 double checkvaldouble = 0;
501 float checkvalfloat = 0;
502 int i;
503
504
505 assert(NULL != raster);
506
507 /* Make sure index is in a valid range */
508 oldnumbands = rt_raster_get_num_bands(raster);
509 if (index < 0)
510 index = 0;
511 else if (index > oldnumbands + 1)
512 index = oldnumbands + 1;
513
514 /* Determine size of memory block to allocate and allocate it */
515 width = rt_raster_get_width(raster);
516 height = rt_raster_get_height(raster);
517 numval = width * height;
518 datasize = rt_pixtype_size(pixtype) * numval;
519
520 mem = (int *)rtalloc(datasize);
521 if (!mem) {
522 rterror("rt_raster_generate_new_band: Could not allocate memory for band");
523 return -1;
524 }
525
526 if (FLT_EQ(initialvalue, 0.0))
527 memset(mem, 0, datasize);
528 else {
529 switch (pixtype)
530 {
531 case PT_1BB:
532 {
533 uint8_t *ptr = mem;
534 uint8_t clamped_initval = rt_util_clamp_to_1BB(initialvalue);
535 for (i = 0; i < numval; i++)
536 ptr[i] = clamped_initval;
537 checkvalint = ptr[0];
538 break;
539 }
540 case PT_2BUI:
541 {
542 uint8_t *ptr = mem;
543 uint8_t clamped_initval = rt_util_clamp_to_2BUI(initialvalue);
544 for (i = 0; i < numval; i++)
545 ptr[i] = clamped_initval;
546 checkvalint = ptr[0];
547 break;
548 }
549 case PT_4BUI:
550 {
551 uint8_t *ptr = mem;
552 uint8_t clamped_initval = rt_util_clamp_to_4BUI(initialvalue);
553 for (i = 0; i < numval; i++)
554 ptr[i] = clamped_initval;
555 checkvalint = ptr[0];
556 break;
557 }
558 case PT_8BSI:
559 {
560 int8_t *ptr = mem;
561 int8_t clamped_initval = rt_util_clamp_to_8BSI(initialvalue);
562 for (i = 0; i < numval; i++)
563 ptr[i] = clamped_initval;
564 checkvalint = ptr[0];
565 break;
566 }
567 case PT_8BUI:
568 {
569 uint8_t *ptr = mem;
570 uint8_t clamped_initval = rt_util_clamp_to_8BUI(initialvalue);
571 for (i = 0; i < numval; i++)
572 ptr[i] = clamped_initval;
573 checkvalint = ptr[0];
574 break;
575 }
576 case PT_16BSI:
577 {
578 int16_t *ptr = mem;
579 int16_t clamped_initval = rt_util_clamp_to_16BSI(initialvalue);
580 for (i = 0; i < numval; i++)
581 ptr[i] = clamped_initval;
582 checkvalint = ptr[0];
583 break;
584 }
585 case PT_16BUI:
586 {
587 uint16_t *ptr = mem;
588 uint16_t clamped_initval = rt_util_clamp_to_16BUI(initialvalue);
589 for (i = 0; i < numval; i++)
590 ptr[i] = clamped_initval;
591 checkvalint = ptr[0];
592 break;
593 }
594 case PT_32BSI:
595 {
596 int32_t *ptr = mem;
597 int32_t clamped_initval = rt_util_clamp_to_32BSI(initialvalue);
598 for (i = 0; i < numval; i++)
599 ptr[i] = clamped_initval;
600 checkvalint = ptr[0];
601 break;
602 }
603 case PT_32BUI:
604 {
605 uint32_t *ptr = mem;
606 uint32_t clamped_initval = rt_util_clamp_to_32BUI(initialvalue);
607 for (i = 0; i < numval; i++)
608 ptr[i] = clamped_initval;
609 checkvaluint = ptr[0];
610 break;
611 }
612 case PT_32BF:
613 {
614 float *ptr = mem;
615 float clamped_initval = rt_util_clamp_to_32F(initialvalue);
616 for (i = 0; i < numval; i++)
617 ptr[i] = clamped_initval;
618 checkvalfloat = ptr[0];
619 break;
620 }
621 case PT_64BF:
622 {
623 double *ptr = mem;
624 for (i = 0; i < numval; i++)
625 ptr[i] = initialvalue;
626 checkvaldouble = ptr[0];
627 break;
628 }
629 default:
630 {
631 rterror("rt_raster_generate_new_band: Unknown pixeltype %d", pixtype);
632 rtdealloc(mem);
633 return -1;
634 }
635 }
636 }
637
638 /* Overflow checking */
639 rt_util_dbl_trunc_warning(
640 initialvalue,
641 checkvalint, checkvaluint,
642 checkvalfloat, checkvaldouble,
643 pixtype
644 );
645
646 band = rt_band_new_inline(width, height, pixtype, hasnodata, nodatavalue, mem);
647 if (! band) {
648 rterror("rt_raster_generate_new_band: Could not add band to raster. Aborting");
649 rtdealloc(mem);
650 return -1;
651 }
652 rt_band_set_ownsdata_flag(band, 1); /* we DO own this data!!! */
653 index = rt_raster_add_band(raster, band, index);
654 numbands = rt_raster_get_num_bands(raster);
655 if (numbands == oldnumbands || index == -1) {
656 rterror("rt_raster_generate_new_band: Could not add band to raster. Aborting");
657 rt_band_destroy(band);
658 }
659
660 /* set isnodata if hasnodata = TRUE and initial value = nodatavalue */
661 if (hasnodata && FLT_EQ(initialvalue, nodatavalue))
662 rt_band_set_isnodata_flag(band, 1);
663
664 return index;
665 }
666
667 /**
668 * Get 6-element array of raster inverse geotransform matrix
669 *
670 * @param raster : the raster to get matrix of
671 * @param gt : optional input parameter, 6-element geotransform matrix
672 * @param igt : output parameter, 6-element inverse geotransform matrix
673 *
674 * @return ES_NONE if success, ES_ERROR if error
675 */
rt_raster_get_inverse_geotransform_matrix(rt_raster raster,double * gt,double * igt)676 rt_errorstate rt_raster_get_inverse_geotransform_matrix(
677 rt_raster raster,
678 double *gt, double *igt
679 ) {
680 double _gt[6] = {0};
681
682 assert((raster != NULL || gt != NULL));
683 assert(igt != NULL);
684
685 if (gt == NULL)
686 rt_raster_get_geotransform_matrix(raster, _gt);
687 else
688 memcpy(_gt, gt, sizeof(double) * 6);
689
690 if (!GDALInvGeoTransform(_gt, igt)) {
691 rterror("rt_raster_get_inverse_geotransform_matrix: Could not compute inverse geotransform matrix");
692 return ES_ERROR;
693 }
694
695 return ES_NONE;
696 }
697
698 /**
699 * Get 6-element array of raster geotransform matrix
700 *
701 * @param raster : the raster to get matrix of
702 * @param gt : output parameter, 6-element geotransform matrix
703 *
704 */
705 void
rt_raster_get_geotransform_matrix(rt_raster raster,double * gt)706 rt_raster_get_geotransform_matrix(rt_raster raster,
707 double *gt) {
708 assert(NULL != raster);
709 assert(NULL != gt);
710
711 gt[0] = raster->ipX;
712 gt[1] = raster->scaleX;
713 gt[2] = raster->skewX;
714 gt[3] = raster->ipY;
715 gt[4] = raster->skewY;
716 gt[5] = raster->scaleY;
717 }
718
719 /**
720 * Set raster's geotransform using 6-element array
721 *
722 * @param raster : the raster to set matrix of
723 * @param gt : intput parameter, 6-element geotransform matrix
724 *
725 */
726 void
rt_raster_set_geotransform_matrix(rt_raster raster,double * gt)727 rt_raster_set_geotransform_matrix(rt_raster raster,
728 double *gt) {
729 assert(NULL != raster);
730 assert(NULL != gt);
731
732 raster->ipX = gt[0];
733 raster->scaleX = gt[1];
734 raster->skewX = gt[2];
735 raster->ipY = gt[3];
736 raster->skewY = gt[4];
737 raster->scaleY = gt[5];
738
739 _rt_raster_geotransform_warn_offline_band(raster);
740 }
741
742 /**
743 * Convert an xr, yr raster point to an xw, yw point on map
744 *
745 * @param raster : the raster to get info from
746 * @param xr : the pixel's column
747 * @param yr : the pixel's row
748 * @param xw : output parameter, X ordinate of the geographical point
749 * @param yw : output parameter, Y ordinate of the geographical point
750 * @param gt : input/output parameter, 3x2 geotransform matrix
751 *
752 * @return ES_NONE if success, ES_ERROR if error
753 */
754 rt_errorstate
rt_raster_cell_to_geopoint(rt_raster raster,double xr,double yr,double * xw,double * yw,double * gt)755 rt_raster_cell_to_geopoint(
756 rt_raster raster,
757 double xr, double yr,
758 double *xw, double *yw,
759 double *gt
760 ) {
761 double _gt[6] = {0};
762
763 assert(NULL != raster);
764 assert(NULL != xw && NULL != yw);
765
766 if (NULL != gt)
767 memcpy(_gt, gt, sizeof(double) * 6);
768
769 /* scale of matrix is not set */
770 if (FLT_EQ(_gt[1], 0.0) || FLT_EQ(_gt[5], 0.0))
771 {
772 rt_raster_get_geotransform_matrix(raster, _gt);
773 }
774
775 RASTER_DEBUGF(4, "gt = (%f, %f, %f, %f, %f, %f)",
776 _gt[0],
777 _gt[1],
778 _gt[2],
779 _gt[3],
780 _gt[4],
781 _gt[5]
782 );
783
784 GDALApplyGeoTransform(_gt, xr, yr, xw, yw);
785 RASTER_DEBUGF(4, "GDALApplyGeoTransform (c -> g) for (%f, %f) = (%f, %f)",
786 xr, yr, *xw, *yw);
787
788 return ES_NONE;
789 }
790
791 /**
792 * Convert an xw,yw map point to a xr,yr raster point
793 *
794 * @param raster : the raster to get info from
795 * @param xw : X ordinate of the geographical point
796 * @param yw : Y ordinate of the geographical point
797 * @param xr : output parameter, the pixel's column
798 * @param yr : output parameter, the pixel's row
799 * @param igt : input/output parameter, inverse geotransform matrix
800 *
801 * @return ES_NONE if success, ES_ERROR if error
802 */
803 rt_errorstate
rt_raster_geopoint_to_cell(rt_raster raster,double xw,double yw,double * xr,double * yr,double * igt)804 rt_raster_geopoint_to_cell(
805 rt_raster raster,
806 double xw, double yw,
807 double *xr, double *yr,
808 double *igt
809 ) {
810 double _igt[6] = {0};
811 double rnd = 0;
812
813 assert(NULL != raster);
814 assert(NULL != xr && NULL != yr);
815
816 if (igt != NULL)
817 memcpy(_igt, igt, sizeof(double) * 6);
818
819 /* matrix is not set */
820 if (
821 FLT_EQ(_igt[0], 0.) &&
822 FLT_EQ(_igt[1], 0.) &&
823 FLT_EQ(_igt[2], 0.) &&
824 FLT_EQ(_igt[3], 0.) &&
825 FLT_EQ(_igt[4], 0.) &&
826 FLT_EQ(_igt[5], 0.)
827 ) {
828 if (rt_raster_get_inverse_geotransform_matrix(raster, NULL, _igt) != ES_NONE) {
829 rterror("rt_raster_geopoint_to_cell: Could not get inverse geotransform matrix");
830 return ES_ERROR;
831 }
832 }
833
834 GDALApplyGeoTransform(_igt, xw, yw, xr, yr);
835 RASTER_DEBUGF(4, "GDALApplyGeoTransform (g -> c) for (%f, %f) = (%f, %f)",
836 xw, yw, *xr, *yr);
837
838 rnd = ROUND(*xr, 0);
839 if (FLT_EQ(rnd, *xr))
840 *xr = rnd;
841 else
842 *xr = floor(*xr);
843
844 rnd = ROUND(*yr, 0);
845 if (FLT_EQ(rnd, *yr))
846 *yr = rnd;
847 else
848 *yr = floor(*yr);
849
850 RASTER_DEBUGF(4, "Corrected GDALApplyGeoTransform (g -> c) for (%f, %f) = (%f, %f)",
851 xw, yw, *xr, *yr);
852
853 return ES_NONE;
854 }
855
856 /******************************************************************************
857 * rt_raster_get_envelope()
858 ******************************************************************************/
859
860 /**
861 * Get raster's envelope.
862 *
863 * The envelope is the minimum bounding rectangle of the raster
864 *
865 * @param raster : the raster to get envelope of
866 * @param env : pointer to rt_envelope
867 *
868 * @return ES_NONE if success, ES_ERROR if error
869 */
870 rt_errorstate
rt_raster_get_envelope(rt_raster raster,rt_envelope * env)871 rt_raster_get_envelope(
872 rt_raster raster,
873 rt_envelope *env
874 ) {
875 int i;
876 int rtn;
877 int set = 0;
878 double _r[2] = {0.};
879 double _w[2] = {0.};
880 double _gt[6] = {0.};
881
882 assert(raster != NULL);
883 assert(env != NULL);
884
885 rt_raster_get_geotransform_matrix(raster, _gt);
886
887 for (i = 0; i < 4; i++) {
888 switch (i) {
889 case 0:
890 _r[0] = 0;
891 _r[1] = 0;
892 break;
893 case 1:
894 _r[0] = 0;
895 _r[1] = raster->height;
896 break;
897 case 2:
898 _r[0] = raster->width;
899 _r[1] = raster->height;
900 break;
901 case 3:
902 _r[0] = raster->width;
903 _r[1] = 0;
904 break;
905 }
906
907 rtn = rt_raster_cell_to_geopoint(
908 raster,
909 _r[0], _r[1],
910 &(_w[0]), &(_w[1]),
911 _gt
912 );
913 if (rtn != ES_NONE) {
914 rterror("rt_raster_get_envelope: Could not compute spatial coordinates for raster pixel");
915 return ES_ERROR;
916 }
917
918 if (!set) {
919 set = 1;
920 env->MinX = _w[0];
921 env->MaxX = _w[0];
922 env->MinY = _w[1];
923 env->MaxY = _w[1];
924 }
925 else {
926 if (_w[0] < env->MinX)
927 env->MinX = _w[0];
928 else if (_w[0] > env->MaxX)
929 env->MaxX = _w[0];
930
931 if (_w[1] < env->MinY)
932 env->MinY = _w[1];
933 else if (_w[1] > env->MaxY)
934 env->MaxY = _w[1];
935 }
936 }
937
938 return ES_NONE;
939 }
940
941 /******************************************************************************
942 * rt_raster_compute_skewed_raster()
943 ******************************************************************************/
944
945 /*
946 * Compute skewed extent that covers unskewed extent.
947 *
948 * @param envelope : unskewed extent of type rt_envelope
949 * @param skew : pointer to 2-element array (x, y) of skew
950 * @param scale : pointer to 2-element array (x, y) of scale
951 * @param tolerance : value between 0 and 1 where the smaller the tolerance
952 * results in an extent approaching the "minimum" skewed extent.
953 * If value <= 0, tolerance = 0.1. If value > 1, tolerance = 1.
954 *
955 * @return skewed raster who's extent covers unskewed extent, NULL on error
956 */
957 rt_raster
rt_raster_compute_skewed_raster(rt_envelope extent,double * skew,double * scale,double tolerance)958 rt_raster_compute_skewed_raster(
959 rt_envelope extent,
960 double *skew,
961 double *scale,
962 double tolerance
963 ) {
964 uint32_t run = 0;
965 uint32_t max_run = 1;
966 double dbl_run = 0;
967
968 int rtn;
969 int covers = 0;
970 rt_raster raster;
971 double _gt[6] = {0};
972 double _igt[6] = {0};
973 int _d[2] = {1, -1};
974 int _dlast = 0;
975 int _dlastpos = 0;
976 double _w[2] = {0};
977 double _r[2] = {0};
978 double _xy[2] = {0};
979 int i;
980 int j;
981 int x;
982 int y;
983
984 LWGEOM *geom = NULL;
985 GEOSGeometry *sgeom = NULL;
986 GEOSGeometry *ngeom = NULL;
987
988 if (
989 (tolerance < 0.) ||
990 FLT_EQ(tolerance, 0.)
991 ) {
992 tolerance = 0.1;
993 }
994 else if (tolerance > 1.)
995 tolerance = 1;
996
997 dbl_run = tolerance;
998 while (dbl_run < 10) {
999 dbl_run *= 10.;
1000 max_run *= 10;
1001 }
1002
1003 /* scale must be provided */
1004 if (scale == NULL)
1005 return NULL;
1006 for (i = 0; i < 2; i++) {
1007 if (FLT_EQ(scale[i], 0.0))
1008 {
1009 rterror("rt_raster_compute_skewed_raster: Scale cannot be zero");
1010 return 0;
1011 }
1012
1013 if (i < 1)
1014 _gt[1] = fabs(scale[i] * tolerance);
1015 else
1016 _gt[5] = fabs(scale[i] * tolerance);
1017 }
1018 /* conform scale-y to be negative */
1019 _gt[5] *= -1;
1020
1021 /* skew not provided or skew is zero, return raster of correct dim and spatial attributes */
1022 if ((skew == NULL) || (FLT_EQ(skew[0], 0.0) && FLT_EQ(skew[1], 0.0)))
1023 {
1024 int _dim[2] = {
1025 (int) fmax((fabs(extent.MaxX - extent.MinX) + (fabs(scale[0]) / 2.)) / fabs(scale[0]), 1),
1026 (int) fmax((fabs(extent.MaxY - extent.MinY) + (fabs(scale[1]) / 2.)) / fabs(scale[1]), 1)
1027 };
1028
1029 raster = rt_raster_new(_dim[0], _dim[1]);
1030 if (raster == NULL) {
1031 rterror("rt_raster_compute_skewed_raster: Could not create output raster");
1032 return NULL;
1033 }
1034
1035 rt_raster_set_offsets(raster, extent.MinX, extent.MaxY);
1036 rt_raster_set_scale(raster, fabs(scale[0]), -1 * fabs(scale[1]));
1037 rt_raster_set_skews(raster, skew[0], skew[1]);
1038
1039 return raster;
1040 }
1041
1042 /* direction to shift upper-left corner */
1043 if (skew[0] > 0.)
1044 _d[0] = -1;
1045 if (skew[1] < 0.)
1046 _d[1] = 1;
1047
1048 /* geotransform */
1049 _gt[0] = extent.UpperLeftX;
1050 _gt[2] = skew[0] * tolerance;
1051 _gt[3] = extent.UpperLeftY;
1052 _gt[4] = skew[1] * tolerance;
1053
1054 RASTER_DEBUGF(4, "Initial geotransform: %f, %f, %f, %f, %f, %f",
1055 _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]
1056 );
1057 RASTER_DEBUGF(4, "Delta: %d, %d", _d[0], _d[1]);
1058
1059 /* simple raster */
1060 if ((raster = rt_raster_new(1, 1)) == NULL) {
1061 rterror("rt_raster_compute_skewed_raster: Out of memory allocating extent raster");
1062 return NULL;
1063 }
1064 rt_raster_set_geotransform_matrix(raster, _gt);
1065
1066 /* get inverse geotransform matrix */
1067 if (!GDALInvGeoTransform(_gt, _igt)) {
1068 rterror("rt_raster_compute_skewed_raster: Could not compute inverse geotransform matrix");
1069 rt_raster_destroy(raster);
1070 return NULL;
1071 }
1072 RASTER_DEBUGF(4, "Inverse geotransform: %f, %f, %f, %f, %f, %f",
1073 _igt[0], _igt[1], _igt[2], _igt[3], _igt[4], _igt[5]
1074 );
1075
1076 /* shift along axis */
1077 for (i = 0; i < 2; i++) {
1078 covers = 0;
1079 run = 0;
1080
1081 RASTER_DEBUGF(3, "Shifting along %s axis", i < 1 ? "X" : "Y");
1082
1083 do {
1084
1085 /* prevent possible infinite loop */
1086 if (run > max_run) {
1087 rterror("rt_raster_compute_skewed_raster: Could not compute skewed extent due to check preventing infinite loop");
1088 rt_raster_destroy(raster);
1089 return NULL;
1090 }
1091
1092 /*
1093 check the four corners that they are covered along the specific axis
1094 pixel column should be >= 0
1095 */
1096 for (j = 0; j < 4; j++) {
1097 switch (j) {
1098 /* upper-left */
1099 case 0:
1100 _xy[0] = extent.MinX;
1101 _xy[1] = extent.MaxY;
1102 break;
1103 /* lower-left */
1104 case 1:
1105 _xy[0] = extent.MinX;
1106 _xy[1] = extent.MinY;
1107 break;
1108 /* lower-right */
1109 case 2:
1110 _xy[0] = extent.MaxX;
1111 _xy[1] = extent.MinY;
1112 break;
1113 /* upper-right */
1114 case 3:
1115 _xy[0] = extent.MaxX;
1116 _xy[1] = extent.MaxY;
1117 break;
1118 }
1119
1120 rtn = rt_raster_geopoint_to_cell(
1121 raster,
1122 _xy[0], _xy[1],
1123 &(_r[0]), &(_r[1]),
1124 _igt
1125 );
1126 if (rtn != ES_NONE) {
1127 rterror("rt_raster_compute_skewed_raster: Could not compute raster pixel for spatial coordinates");
1128 rt_raster_destroy(raster);
1129 return NULL;
1130 }
1131
1132 RASTER_DEBUGF(4, "Point %d at cell %d x %d", j, (int) _r[0], (int) _r[1]);
1133
1134 /* raster doesn't cover point */
1135 if ((int) _r[i] < 0) {
1136 RASTER_DEBUGF(4, "Point outside of skewed extent: %d", j);
1137 covers = 0;
1138
1139 if (_dlastpos != j) {
1140 _dlast = (int) _r[i];
1141 _dlastpos = j;
1142 }
1143 else if ((int) _r[i] < _dlast) {
1144 RASTER_DEBUG(4, "Point going in wrong direction. Reversing direction");
1145 _d[i] *= -1;
1146 _dlastpos = -1;
1147 run = 0;
1148 }
1149
1150 break;
1151 }
1152
1153 covers++;
1154 }
1155
1156 if (!covers) {
1157 x = 0;
1158 y = 0;
1159 if (i < 1)
1160 x = _d[i] * fabs(_r[i]);
1161 else
1162 y = _d[i] * fabs(_r[i]);
1163
1164 rtn = rt_raster_cell_to_geopoint(
1165 raster,
1166 x, y,
1167 &(_w[0]), &(_w[1]),
1168 _gt
1169 );
1170 if (rtn != ES_NONE) {
1171 rterror("rt_raster_compute_skewed_raster: Could not compute spatial coordinates for raster pixel");
1172 rt_raster_destroy(raster);
1173 return NULL;
1174 }
1175
1176 /* adjust ul */
1177 if (i < 1)
1178 _gt[0] = _w[i];
1179 else
1180 _gt[3] = _w[i];
1181 rt_raster_set_geotransform_matrix(raster, _gt);
1182 RASTER_DEBUGF(4, "Shifted geotransform: %f, %f, %f, %f, %f, %f",
1183 _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]
1184 );
1185
1186 /* get inverse geotransform matrix */
1187 if (!GDALInvGeoTransform(_gt, _igt)) {
1188 rterror("rt_raster_compute_skewed_raster: Could not compute inverse geotransform matrix");
1189 rt_raster_destroy(raster);
1190 return NULL;
1191 }
1192 RASTER_DEBUGF(4, "Inverse geotransform: %f, %f, %f, %f, %f, %f",
1193 _igt[0], _igt[1], _igt[2], _igt[3], _igt[4], _igt[5]
1194 );
1195 }
1196
1197 run++;
1198 }
1199 while (!covers);
1200 }
1201
1202 /* covers test */
1203 rtn = rt_raster_geopoint_to_cell(
1204 raster,
1205 extent.MaxX, extent.MinY,
1206 &(_r[0]), &(_r[1]),
1207 _igt
1208 );
1209 if (rtn != ES_NONE) {
1210 rterror("rt_raster_compute_skewed_raster: Could not compute raster pixel for spatial coordinates");
1211 rt_raster_destroy(raster);
1212 return NULL;
1213 }
1214
1215 RASTER_DEBUGF(4, "geopoint %f x %f at cell %d x %d", extent.MaxX, extent.MinY, (int) _r[0], (int) _r[1]);
1216
1217 raster->width = _r[0];
1218 raster->height = _r[1];
1219
1220 /* initialize GEOS */
1221 initGEOS(rtinfo, lwgeom_geos_error);
1222
1223 /* create reference LWPOLY */
1224 {
1225 LWPOLY *npoly = rt_util_envelope_to_lwpoly(extent);
1226 if (npoly == NULL) {
1227 rterror("rt_raster_compute_skewed_raster: Could not build extent's geometry for covers test");
1228 rt_raster_destroy(raster);
1229 return NULL;
1230 }
1231
1232 ngeom = (GEOSGeometry *) LWGEOM2GEOS(lwpoly_as_lwgeom(npoly), 0);
1233 lwpoly_free(npoly);
1234 }
1235
1236 do {
1237 covers = 0;
1238
1239 /* construct sgeom from raster */
1240 if ((rt_raster_get_convex_hull(raster, &geom) != ES_NONE) || geom == NULL) {
1241 rterror("rt_raster_compute_skewed_raster: Could not build skewed extent's geometry for covers test");
1242 GEOSGeom_destroy(ngeom);
1243 rt_raster_destroy(raster);
1244 return NULL;
1245 }
1246
1247 sgeom = (GEOSGeometry *) LWGEOM2GEOS(geom, 0);
1248 lwgeom_free(geom);
1249
1250 covers = GEOSRelatePattern(sgeom, ngeom, "******FF*");
1251 GEOSGeom_destroy(sgeom);
1252
1253 if (covers == 2) {
1254 rterror("rt_raster_compute_skewed_raster: Could not run covers test");
1255 GEOSGeom_destroy(ngeom);
1256 rt_raster_destroy(raster);
1257 return NULL;
1258 }
1259
1260 if (!covers)
1261 {
1262 raster->width++;
1263 raster->height++;
1264 }
1265 }
1266 while (!covers);
1267
1268 RASTER_DEBUGF(4, "Skewed extent does cover normal extent with dimensions %d x %d", raster->width, raster->height);
1269
1270 raster->width = (int) ((((double) raster->width) * fabs(_gt[1]) + fabs(scale[0] / 2.)) / fabs(scale[0]));
1271 raster->height = (int) ((((double) raster->height) * fabs(_gt[5]) + fabs(scale[1] / 2.)) / fabs(scale[1]));
1272 _gt[1] = fabs(scale[0]);
1273 _gt[5] = -1 * fabs(scale[1]);
1274 _gt[2] = skew[0];
1275 _gt[4] = skew[1];
1276 rt_raster_set_geotransform_matrix(raster, _gt);
1277
1278 /* minimize width/height */
1279 for (i = 0; i < 2; i++) {
1280 covers = 1;
1281 do {
1282 if (i < 1)
1283 raster->width--;
1284 else
1285 raster->height--;
1286
1287 /* construct sgeom from raster */
1288 if ((rt_raster_get_convex_hull(raster, &geom) != ES_NONE) || geom == NULL) {
1289 rterror("rt_raster_compute_skewed_raster: Could not build skewed extent's geometry for minimizing dimensions");
1290 GEOSGeom_destroy(ngeom);
1291 rt_raster_destroy(raster);
1292 return NULL;
1293 }
1294
1295 sgeom = (GEOSGeometry *) LWGEOM2GEOS(geom, 0);
1296 lwgeom_free(geom);
1297
1298 covers = GEOSRelatePattern(sgeom, ngeom, "******FF*");
1299 GEOSGeom_destroy(sgeom);
1300
1301 if (covers == 2) {
1302 rterror("rt_raster_compute_skewed_raster: Could not run covers test for minimizing dimensions");
1303 GEOSGeom_destroy(ngeom);
1304 rt_raster_destroy(raster);
1305 return NULL;
1306 }
1307 } while (covers);
1308
1309 if (i < 1)
1310 raster->width++;
1311 else
1312 raster->height++;
1313
1314 }
1315
1316 GEOSGeom_destroy(ngeom);
1317
1318 return raster;
1319 }
1320
1321 /**
1322 * Return TRUE if the raster is empty. i.e. is NULL, width = 0 or height = 0
1323 *
1324 * @param raster : the raster to get info from
1325 *
1326 * @return TRUE if the raster is empty, FALSE otherwise
1327 */
1328 int
rt_raster_is_empty(rt_raster raster)1329 rt_raster_is_empty(rt_raster raster) {
1330 return (!raster || raster->height == 0 || raster->width == 0);
1331 }
1332
1333 /**
1334 * Return TRUE if the raster has a band of this number.
1335 *
1336 * @param raster : the raster to get info from
1337 * @param nband : the band number. 0-based
1338 *
1339 * @return TRUE if the raster has a band of this number, FALSE otherwise
1340 */
1341 int
rt_raster_has_band(rt_raster raster,int nband)1342 rt_raster_has_band(rt_raster raster, int nband) {
1343 return !(NULL == raster || nband >= raster->numBands || nband < 0);
1344 }
1345
1346 /******************************************************************************
1347 * rt_raster_copy_band()
1348 ******************************************************************************/
1349
1350 /**
1351 * Copy one band from one raster to another. Bands are duplicated from
1352 * fromrast to torast using rt_band_duplicate. The caller will need
1353 * to ensure that the copied band's data or path remains allocated
1354 * for the lifetime of the copied bands.
1355 *
1356 * @param torast : raster to copy band to
1357 * @param fromrast : raster to copy band from
1358 * @param fromindex : index of band in source raster, 0-based
1359 * @param toindex : index of new band in destination raster, 0-based
1360 *
1361 * @return The band index of the second raster where the new band is copied.
1362 * -1 if error
1363 */
1364 int
rt_raster_copy_band(rt_raster torast,rt_raster fromrast,int fromindex,int toindex)1365 rt_raster_copy_band(
1366 rt_raster torast, rt_raster fromrast,
1367 int fromindex, int toindex
1368 ) {
1369 rt_band srcband = NULL;
1370 rt_band dstband = NULL;
1371
1372 assert(NULL != torast);
1373 assert(NULL != fromrast);
1374
1375 /* Check raster dimensions */
1376 if (torast->height != fromrast->height || torast->width != fromrast->width) {
1377 rtwarn("rt_raster_copy_band: Attempting to add a band with different width or height");
1378 return -1;
1379 }
1380
1381 /* Check bands limits */
1382 if (fromrast->numBands < 1) {
1383 rtwarn("rt_raster_copy_band: Second raster has no band");
1384 return -1;
1385 }
1386 else if (fromindex < 0) {
1387 rtwarn("rt_raster_copy_band: Band index for second raster < 0. Defaulted to 0");
1388 fromindex = 0;
1389 }
1390 else if (fromindex >= fromrast->numBands) {
1391 rtwarn("rt_raster_copy_band: Band index for second raster > number of bands, truncated from %u to %u", fromindex, fromrast->numBands - 1);
1392 fromindex = fromrast->numBands - 1;
1393 }
1394
1395 if (toindex < 0) {
1396 rtwarn("rt_raster_copy_band: Band index for first raster < 0. Defaulted to 0");
1397 toindex = 0;
1398 }
1399 else if (toindex > torast->numBands) {
1400 rtwarn("rt_raster_copy_band: Band index for first raster > number of bands, truncated from %u to %u", toindex, torast->numBands);
1401 toindex = torast->numBands;
1402 }
1403
1404 /* Get band from source raster */
1405 srcband = rt_raster_get_band(fromrast, fromindex);
1406
1407 /* duplicate band */
1408 dstband = rt_band_duplicate(srcband);
1409
1410 /* Add band to the second raster */
1411 return rt_raster_add_band(torast, dstband, toindex);
1412 }
1413
1414 /******************************************************************************
1415 * rt_raster_from_band()
1416 ******************************************************************************/
1417
1418 /**
1419 * Construct a new rt_raster from an existing rt_raster and an array
1420 * of band numbers
1421 *
1422 * @param raster : the source raster
1423 * @param bandNums : array of band numbers to extract from source raster
1424 * and add to the new raster (0 based)
1425 * @param count : number of elements in bandNums
1426 *
1427 * @return a new rt_raster or NULL on error
1428 */
1429 rt_raster
rt_raster_from_band(rt_raster raster,uint32_t * bandNums,int count)1430 rt_raster_from_band(rt_raster raster, uint32_t *bandNums, int count) {
1431 rt_raster rast = NULL;
1432 int i = 0;
1433 int j = 0;
1434 int idx;
1435 int32_t flag;
1436 double gt[6] = {0.};
1437
1438 assert(NULL != raster);
1439 assert(NULL != bandNums);
1440
1441 RASTER_DEBUGF(3, "rt_raster_from_band: source raster has %d bands",
1442 rt_raster_get_num_bands(raster));
1443
1444 /* create new raster */
1445 rast = rt_raster_new(raster->width, raster->height);
1446 if (NULL == rast) {
1447 rterror("rt_raster_from_band: Out of memory allocating new raster");
1448 return NULL;
1449 }
1450
1451 /* copy raster attributes */
1452 rt_raster_get_geotransform_matrix(raster, gt);
1453 rt_raster_set_geotransform_matrix(rast, gt);
1454
1455 /* srid */
1456 rt_raster_set_srid(rast, raster->srid);
1457
1458 /* copy bands */
1459 for (i = 0; i < count; i++) {
1460 idx = bandNums[i];
1461 flag = rt_raster_copy_band(rast, raster, idx, i);
1462
1463 if (flag < 0) {
1464 rterror("rt_raster_from_band: Could not copy band");
1465 for (j = 0; j < i; j++) rt_band_destroy(rast->bands[j]);
1466 rt_raster_destroy(rast);
1467 return NULL;
1468 }
1469
1470 RASTER_DEBUGF(3, "rt_raster_from_band: band created at index %d",
1471 flag);
1472 }
1473
1474 RASTER_DEBUGF(3, "rt_raster_from_band: new raster has %d bands",
1475 rt_raster_get_num_bands(rast));
1476 return rast;
1477 }
1478
1479 /******************************************************************************
1480 * rt_raster_replace_band()
1481 ******************************************************************************/
1482
1483 /**
1484 * Replace band at provided index with new band
1485 *
1486 * @param raster: raster of band to be replaced
1487 * @param band : new band to add to raster
1488 * @param index : index of band to replace (0-based)
1489 *
1490 * @return NULL on error or replaced band
1491 */
1492 rt_band
rt_raster_replace_band(rt_raster raster,rt_band band,int index)1493 rt_raster_replace_band(rt_raster raster, rt_band band, int index) {
1494 rt_band oldband = NULL;
1495 assert(NULL != raster);
1496 assert(NULL != band);
1497
1498 if (band->width != raster->width || band->height != raster->height) {
1499 rterror("rt_raster_replace_band: Band does not match raster's dimensions: %dx%d band to %dx%d raster",
1500 band->width, band->height, raster->width, raster->height);
1501 return 0;
1502 }
1503
1504 if (index >= raster->numBands || index < 0) {
1505 rterror("rt_raster_replace_band: Band index is not valid");
1506 return 0;
1507 }
1508
1509 oldband = rt_raster_get_band(raster, index);
1510 RASTER_DEBUGF(3, "rt_raster_replace_band: old band at %p", oldband);
1511 RASTER_DEBUGF(3, "rt_raster_replace_band: new band at %p", band);
1512
1513 raster->bands[index] = band;
1514 RASTER_DEBUGF(3, "rt_raster_replace_band: new band at %p", raster->bands[index]);
1515
1516 band->raster = raster;
1517 oldband->raster = NULL;
1518
1519 return oldband;
1520 }
1521
1522 /******************************************************************************
1523 * rt_raster_clone()
1524 ******************************************************************************/
1525
1526 /**
1527 * Clone an existing raster
1528 *
1529 * @param raster : raster to clone
1530 * @param deep : flag indicating if bands should be cloned
1531 *
1532 * @return a new rt_raster or NULL on error
1533 */
1534 rt_raster
rt_raster_clone(rt_raster raster,uint8_t deep)1535 rt_raster_clone(rt_raster raster, uint8_t deep) {
1536 rt_raster rtn = NULL;
1537 double gt[6] = {0};
1538
1539 assert(NULL != raster);
1540
1541 if (deep) {
1542 int numband = rt_raster_get_num_bands(raster);
1543 uint32_t *nband = NULL;
1544 int i = 0;
1545
1546 nband = rtalloc(sizeof(uint32_t) * numband);
1547 if (nband == NULL) {
1548 rterror("rt_raster_clone: Could not allocate memory for deep clone");
1549 return NULL;
1550 }
1551 for (i = 0; i < numband; i++)
1552 nband[i] = i;
1553
1554 rtn = rt_raster_from_band(raster, nband, numband);
1555 rtdealloc(nband);
1556
1557 return rtn;
1558 }
1559
1560 rtn = rt_raster_new(
1561 rt_raster_get_width(raster),
1562 rt_raster_get_height(raster)
1563 );
1564 if (rtn == NULL) {
1565 rterror("rt_raster_clone: Could not create cloned raster");
1566 return NULL;
1567 }
1568
1569 rt_raster_get_geotransform_matrix(raster, gt);
1570 rt_raster_set_geotransform_matrix(rtn, gt);
1571 rt_raster_set_srid(rtn, rt_raster_get_srid(raster));
1572
1573 return rtn;
1574 }
1575
1576 /******************************************************************************
1577 * rt_raster_to_gdal()
1578 ******************************************************************************/
1579
1580 /**
1581 * Return formatted GDAL raster from raster
1582 *
1583 * @param raster : the raster to convert
1584 * @param srs : the raster's coordinate system in OGC WKT
1585 * @param format : format to convert to. GDAL driver short name
1586 * @param options : list of format creation options. array of strings
1587 * @param gdalsize : will be set to the size of returned bytea
1588 *
1589 * @return formatted GDAL raster. the calling function is responsible
1590 * for freeing the returned data using CPLFree()
1591 */
1592 uint8_t*
rt_raster_to_gdal(rt_raster raster,const char * srs,char * format,char ** options,uint64_t * gdalsize)1593 rt_raster_to_gdal(
1594 rt_raster raster, const char *srs,
1595 char *format, char **options, uint64_t *gdalsize
1596 ) {
1597 const char *cc;
1598 const char *vio;
1599
1600 GDALDriverH src_drv = NULL;
1601 int destroy_src_drv = 0;
1602 GDALDatasetH src_ds = NULL;
1603
1604 vsi_l_offset rtn_lenvsi;
1605 uint64_t rtn_len = 0;
1606
1607 GDALDriverH rtn_drv = NULL;
1608 GDALDatasetH rtn_ds = NULL;
1609 uint8_t *rtn = NULL;
1610
1611 assert(NULL != raster);
1612 assert(NULL != gdalsize);
1613
1614 /* any supported format is possible */
1615 rt_util_gdal_register_all(0);
1616 RASTER_DEBUG(3, "loaded all supported GDAL formats");
1617
1618 /* output format not specified */
1619 if (format == NULL || !strlen(format))
1620 format = "GTiff";
1621 RASTER_DEBUGF(3, "output format is %s", format);
1622
1623 /* load raster into a GDAL MEM raster */
1624 src_ds = rt_raster_to_gdal_mem(raster, srs, NULL, NULL, 0, &src_drv, &destroy_src_drv);
1625 if (NULL == src_ds) {
1626 rterror("rt_raster_to_gdal: Could not convert raster to GDAL MEM format");
1627 return 0;
1628 }
1629
1630 /* load driver */
1631 rtn_drv = GDALGetDriverByName(format);
1632 if (NULL == rtn_drv) {
1633 rterror("rt_raster_to_gdal: Could not load the output GDAL driver");
1634 GDALClose(src_ds);
1635 if (destroy_src_drv) GDALDestroyDriver(src_drv);
1636 return 0;
1637 }
1638 RASTER_DEBUG(3, "Output driver loaded");
1639
1640 /* CreateCopy support */
1641 cc = GDALGetMetadataItem(rtn_drv, GDAL_DCAP_CREATECOPY, NULL);
1642 /* VirtualIO support */
1643 vio = GDALGetMetadataItem(rtn_drv, GDAL_DCAP_VIRTUALIO, NULL);
1644
1645 if (cc == NULL || vio == NULL) {
1646 rterror("rt_raster_to_gdal: Output GDAL driver does not support CreateCopy and/or VirtualIO");
1647 GDALClose(src_ds);
1648 if (destroy_src_drv) GDALDestroyDriver(src_drv);
1649 return 0;
1650 }
1651
1652 /* convert GDAL MEM raster to output format */
1653 RASTER_DEBUG(3, "Copying GDAL MEM raster to memory file in output format");
1654 rtn_ds = GDALCreateCopy(
1655 rtn_drv,
1656 "/vsimem/out.dat", /* should be fine assuming this is in a process */
1657 src_ds,
1658 FALSE, /* should copy be strictly equivelent? */
1659 options, /* format options */
1660 NULL, /* progress function */
1661 NULL /* progress data */
1662 );
1663
1664 /* close source dataset */
1665 GDALClose(src_ds);
1666 if (destroy_src_drv) GDALDestroyDriver(src_drv);
1667 RASTER_DEBUG(3, "Closed GDAL MEM raster");
1668
1669 if (NULL == rtn_ds) {
1670 rterror("rt_raster_to_gdal: Could not create the output GDAL dataset");
1671 return 0;
1672 }
1673
1674 RASTER_DEBUGF(4, "dataset SRS: %s", GDALGetProjectionRef(rtn_ds));
1675
1676 /* close dataset, this also flushes any pending writes */
1677 GDALClose(rtn_ds);
1678 RASTER_DEBUG(3, "Closed GDAL output raster");
1679
1680 RASTER_DEBUG(3, "Done copying GDAL MEM raster to memory file in output format");
1681
1682 /* from memory file to buffer */
1683 RASTER_DEBUG(3, "Copying GDAL memory file to buffer");
1684 rtn = VSIGetMemFileBuffer("/vsimem/out.dat", &rtn_lenvsi, TRUE);
1685 RASTER_DEBUG(3, "Done copying GDAL memory file to buffer");
1686 if (NULL == rtn) {
1687 rterror("rt_raster_to_gdal: Could not create the output GDAL raster");
1688 return 0;
1689 }
1690
1691 rtn_len = (uint64_t) rtn_lenvsi;
1692 *gdalsize = rtn_len;
1693
1694 return rtn;
1695 }
1696
1697 /******************************************************************************
1698 * rt_raster_gdal_drivers()
1699 ******************************************************************************/
1700
1701 /**
1702 * Returns a set of available GDAL drivers
1703 *
1704 * @param drv_count : number of GDAL drivers available
1705 * @param can_write : if non-zero, filter drivers to only those
1706 * with support for CreateCopy and VirtualIO
1707 *
1708 * @return set of "gdaldriver" values of available GDAL drivers
1709 */
1710 rt_gdaldriver
rt_raster_gdal_drivers(uint32_t * drv_count,uint8_t can_write)1711 rt_raster_gdal_drivers(uint32_t *drv_count, uint8_t can_write) {
1712 const char *cc;
1713 const char *vio;
1714 const char *txt;
1715 int txt_len;
1716 GDALDriverH *drv = NULL;
1717 rt_gdaldriver rtn = NULL;
1718 int count;
1719 int i;
1720 uint32_t j;
1721
1722 assert(drv_count != NULL);
1723
1724 rt_util_gdal_register_all(0);
1725 count = GDALGetDriverCount();
1726 RASTER_DEBUGF(3, "%d drivers found", count);
1727
1728 rtn = (rt_gdaldriver) rtalloc(count * sizeof(struct rt_gdaldriver_t));
1729 if (NULL == rtn) {
1730 rterror("rt_raster_gdal_drivers: Could not allocate memory for gdaldriver structure");
1731 return 0;
1732 }
1733
1734 for (i = 0, j = 0; i < count; i++) {
1735 drv = GDALGetDriver(i);
1736
1737 #ifdef GDAL_DCAP_RASTER
1738 /* Starting with GDAL 2.0, vector drivers can also be returned */
1739 /* Only keep raster drivers */
1740 const char *is_raster;
1741 is_raster = GDALGetMetadataItem(drv, GDAL_DCAP_RASTER, NULL);
1742 if (is_raster == NULL || !EQUAL(is_raster, "YES"))
1743 continue;
1744 #endif
1745
1746 /* CreateCopy support */
1747 cc = GDALGetMetadataItem(drv, GDAL_DCAP_CREATECOPY, NULL);
1748
1749 /* VirtualIO support */
1750 vio = GDALGetMetadataItem(drv, GDAL_DCAP_VIRTUALIO, NULL);
1751
1752 if (can_write && (cc == NULL || vio == NULL))
1753 continue;
1754
1755 /* we can always read what GDAL can load */
1756 rtn[j].can_read = 1;
1757 /* we require CreateCopy and VirtualIO support to write to GDAL */
1758 rtn[j].can_write = (cc != NULL && vio != NULL);
1759
1760 if (rtn[j].can_write) {
1761 RASTER_DEBUGF(3, "driver %s (%d) supports CreateCopy() and VirtualIO()", txt, i);
1762 }
1763
1764 /* index of driver */
1765 rtn[j].idx = i;
1766
1767 /* short name */
1768 txt = GDALGetDriverShortName(drv);
1769 txt_len = strlen(txt);
1770
1771 txt_len = (txt_len + 1) * sizeof(char);
1772 rtn[j].short_name = (char *) rtalloc(txt_len);
1773 memcpy(rtn[j].short_name, txt, txt_len);
1774
1775 /* long name */
1776 txt = GDALGetDriverLongName(drv);
1777 txt_len = strlen(txt);
1778
1779 txt_len = (txt_len + 1) * sizeof(char);
1780 rtn[j].long_name = (char *) rtalloc(txt_len);
1781 memcpy(rtn[j].long_name, txt, txt_len);
1782
1783 /* creation options */
1784 txt = GDALGetDriverCreationOptionList(drv);
1785 txt_len = strlen(txt);
1786
1787 txt_len = (txt_len + 1) * sizeof(char);
1788 rtn[j].create_options = (char *) rtalloc(txt_len);
1789 memcpy(rtn[j].create_options, txt, txt_len);
1790
1791 j++;
1792 }
1793
1794 /* free unused memory */
1795 rtn = rtrealloc(rtn, j * sizeof(struct rt_gdaldriver_t));
1796 *drv_count = j;
1797
1798 return rtn;
1799 }
1800
1801 /******************************************************************************
1802 * rt_raster_to_gdal_mem()
1803 ******************************************************************************/
1804
1805 /**
1806 * Return GDAL dataset using GDAL MEM driver from raster.
1807 *
1808 * @param raster : raster to convert to GDAL MEM
1809 * @param srs : the raster's coordinate system in OGC WKT
1810 * @param bandNums : array of band numbers to extract from raster
1811 * and include in the GDAL dataset (0 based)
1812 * @param excludeNodataValues : array of zero, nonzero where if non-zero,
1813 * ignore nodata values for the band
1814 * @param count : number of elements in bandNums
1815 * @param rtn_drv : is set to the GDAL driver object
1816 * @param destroy_rtn_drv : if non-zero, caller must destroy the MEM driver
1817 *
1818 * @return GDAL dataset using GDAL MEM driver
1819 */
1820 GDALDatasetH
rt_raster_to_gdal_mem(rt_raster raster,const char * srs,uint32_t * bandNums,int * excludeNodataValues,int count,GDALDriverH * rtn_drv,int * destroy_rtn_drv)1821 rt_raster_to_gdal_mem(
1822 rt_raster raster,
1823 const char *srs,
1824 uint32_t *bandNums,
1825 int *excludeNodataValues,
1826 int count,
1827 GDALDriverH *rtn_drv, int *destroy_rtn_drv
1828 ) {
1829 GDALDriverH drv = NULL;
1830 GDALDatasetH ds = NULL;
1831 double gt[6] = {0.0};
1832 CPLErr cplerr;
1833 GDALDataType gdal_pt = GDT_Unknown;
1834 GDALRasterBandH band;
1835 void *pVoid;
1836 char *pszDataPointer;
1837 char szGDALOption[50];
1838 char *apszOptions[4];
1839 double nodata = 0.0;
1840 int allocBandNums = 0;
1841 int allocNodataValues = 0;
1842
1843 int i;
1844 uint32_t numBands;
1845 uint32_t width = 0;
1846 uint32_t height = 0;
1847 rt_band rtband = NULL;
1848 rt_pixtype pt = PT_END;
1849
1850 assert(NULL != raster);
1851 assert(NULL != rtn_drv);
1852 assert(NULL != destroy_rtn_drv);
1853
1854 *destroy_rtn_drv = 0;
1855
1856 /* store raster in GDAL MEM raster */
1857 if (!rt_util_gdal_driver_registered("MEM")) {
1858 RASTER_DEBUG(4, "Registering MEM driver");
1859 GDALRegister_MEM();
1860 *destroy_rtn_drv = 1;
1861 }
1862 drv = GDALGetDriverByName("MEM");
1863 if (NULL == drv) {
1864 rterror("rt_raster_to_gdal_mem: Could not load the MEM GDAL driver");
1865 return 0;
1866 }
1867 *rtn_drv = drv;
1868
1869 /* unload driver from GDAL driver manager */
1870 if (*destroy_rtn_drv) {
1871 RASTER_DEBUG(4, "Deregistering MEM driver");
1872 GDALDeregisterDriver(drv);
1873 }
1874
1875 width = rt_raster_get_width(raster);
1876 height = rt_raster_get_height(raster);
1877 ds = GDALCreate(
1878 drv, "",
1879 width, height,
1880 0, GDT_Byte, NULL
1881 );
1882 if (NULL == ds) {
1883 rterror("rt_raster_to_gdal_mem: Could not create a GDALDataset to convert into");
1884 return 0;
1885 }
1886
1887 /* add geotransform */
1888 rt_raster_get_geotransform_matrix(raster, gt);
1889 cplerr = GDALSetGeoTransform(ds, gt);
1890 if (cplerr != CE_None) {
1891 rterror("rt_raster_to_gdal_mem: Could not set geotransformation");
1892 GDALClose(ds);
1893 return 0;
1894 }
1895
1896 /* set spatial reference */
1897 if (NULL != srs && strlen(srs)) {
1898 char *_srs = rt_util_gdal_convert_sr(srs, 0);
1899 if (_srs == NULL) {
1900 rterror("rt_raster_to_gdal_mem: Could not convert srs to GDAL accepted format");
1901 GDALClose(ds);
1902 return 0;
1903 }
1904
1905 cplerr = GDALSetProjection(ds, _srs);
1906 CPLFree(_srs);
1907 if (cplerr != CE_None) {
1908 rterror("rt_raster_to_gdal_mem: Could not set projection");
1909 GDALClose(ds);
1910 return 0;
1911 }
1912 RASTER_DEBUGF(3, "Projection set to: %s", GDALGetProjectionRef(ds));
1913 }
1914
1915 /* process bandNums */
1916 numBands = rt_raster_get_num_bands(raster);
1917 if (NULL != bandNums && count > 0) {
1918 for (i = 0; i < count; i++) {
1919 if (bandNums[i] >= numBands) {
1920 rterror("rt_raster_to_gdal_mem: The band index %d is invalid", bandNums[i]);
1921 GDALClose(ds);
1922 return 0;
1923 }
1924 }
1925 }
1926 else {
1927 count = numBands;
1928 bandNums = (uint32_t *) rtalloc(sizeof(uint32_t) * count);
1929 if (NULL == bandNums) {
1930 rterror("rt_raster_to_gdal_mem: Could not allocate memory for band indices");
1931 GDALClose(ds);
1932 return 0;
1933 }
1934 allocBandNums = 1;
1935 for (i = 0; i < count; i++) bandNums[i] = i;
1936 }
1937
1938 /* process exclude_nodata_values */
1939 if (NULL == excludeNodataValues) {
1940 excludeNodataValues = (int *) rtalloc(sizeof(int) * count);
1941 if (NULL == excludeNodataValues) {
1942 rterror("rt_raster_to_gdal_mem: Could not allocate memory for NODATA flags");
1943 GDALClose(ds);
1944 return 0;
1945 }
1946 allocNodataValues = 1;
1947 for (i = 0; i < count; i++) excludeNodataValues[i] = 1;
1948 }
1949
1950 /* add band(s) */
1951 for (i = 0; i < count; i++) {
1952 rtband = rt_raster_get_band(raster, bandNums[i]);
1953 if (NULL == rtband) {
1954 rterror("rt_raster_to_gdal_mem: Could not get requested band index %d", bandNums[i]);
1955 if (allocBandNums) rtdealloc(bandNums);
1956 if (allocNodataValues) rtdealloc(excludeNodataValues);
1957 GDALClose(ds);
1958 return 0;
1959 }
1960
1961 pt = rt_band_get_pixtype(rtband);
1962 gdal_pt = rt_util_pixtype_to_gdal_datatype(pt);
1963 if (gdal_pt == GDT_Unknown)
1964 rtwarn("rt_raster_to_gdal_mem: Unknown pixel type for band");
1965
1966 /*
1967 For all pixel types other than PT_8BSI, set pointer to start of data
1968 */
1969 if (pt != PT_8BSI) {
1970 pVoid = rt_band_get_data(rtband);
1971 RASTER_DEBUGF(4, "Band data is at pos %p", pVoid);
1972
1973 pszDataPointer = (char *) rtalloc(20 * sizeof (char));
1974 sprintf(pszDataPointer, "%p", pVoid);
1975 RASTER_DEBUGF(4, "rt_raster_to_gdal_mem: szDatapointer is %p",
1976 pszDataPointer);
1977
1978 if (strncasecmp(pszDataPointer, "0x", 2) == 0)
1979 sprintf(szGDALOption, "DATAPOINTER=%s", pszDataPointer);
1980 else
1981 sprintf(szGDALOption, "DATAPOINTER=0x%s", pszDataPointer);
1982
1983 RASTER_DEBUG(3, "Storing info for GDAL MEM raster band");
1984
1985 apszOptions[0] = szGDALOption;
1986 apszOptions[1] = NULL; /* pixel offset, not needed */
1987 apszOptions[2] = NULL; /* line offset, not needed */
1988 apszOptions[3] = NULL;
1989
1990 /* free */
1991 rtdealloc(pszDataPointer);
1992
1993 /* add band */
1994 if (GDALAddBand(ds, gdal_pt, apszOptions) == CE_Failure) {
1995 rterror("rt_raster_to_gdal_mem: Could not add GDAL raster band");
1996 if (allocBandNums) rtdealloc(bandNums);
1997 GDALClose(ds);
1998 return 0;
1999 }
2000 }
2001 /*
2002 PT_8BSI is special as GDAL has no equivalent pixel type.
2003 Must convert 8BSI to 16BSI so create basic band
2004 */
2005 else {
2006 /* add band */
2007 if (GDALAddBand(ds, gdal_pt, NULL) == CE_Failure) {
2008 rterror("rt_raster_to_gdal_mem: Could not add GDAL raster band");
2009 if (allocBandNums) rtdealloc(bandNums);
2010 if (allocNodataValues) rtdealloc(excludeNodataValues);
2011 GDALClose(ds);
2012 return 0;
2013 }
2014 }
2015
2016 /* check band count */
2017 if (GDALGetRasterCount(ds) != i + 1) {
2018 rterror("rt_raster_to_gdal_mem: Error creating GDAL MEM raster band");
2019 if (allocBandNums) rtdealloc(bandNums);
2020 if (allocNodataValues) rtdealloc(excludeNodataValues);
2021 GDALClose(ds);
2022 return 0;
2023 }
2024
2025 /* get new band */
2026 band = NULL;
2027 band = GDALGetRasterBand(ds, i + 1);
2028 if (NULL == band) {
2029 rterror("rt_raster_to_gdal_mem: Could not get GDAL band for additional processing");
2030 if (allocBandNums) rtdealloc(bandNums);
2031 if (allocNodataValues) rtdealloc(excludeNodataValues);
2032 GDALClose(ds);
2033 return 0;
2034 }
2035
2036 /* PT_8BSI requires manual setting of pixels */
2037 if (pt == PT_8BSI) {
2038 uint32_t nXBlocks, nYBlocks;
2039 int nXBlockSize, nYBlockSize;
2040 uint32_t iXBlock, iYBlock;
2041 int nXValid, nYValid;
2042 int iX, iY;
2043 int iXMax, iYMax;
2044
2045 int x, y, z;
2046 uint32_t valueslen = 0;
2047 int16_t *values = NULL;
2048 double value = 0.;
2049
2050 /* this makes use of GDAL's "natural" blocks */
2051 GDALGetBlockSize(band, &nXBlockSize, &nYBlockSize);
2052 nXBlocks = (width + nXBlockSize - 1) / nXBlockSize;
2053 nYBlocks = (height + nYBlockSize - 1) / nYBlockSize;
2054 RASTER_DEBUGF(4, "(nXBlockSize, nYBlockSize) = (%d, %d)", nXBlockSize, nYBlockSize);
2055 RASTER_DEBUGF(4, "(nXBlocks, nYBlocks) = (%d, %d)", nXBlocks, nYBlocks);
2056
2057 /* length is for the desired pixel type */
2058 valueslen = rt_pixtype_size(PT_16BSI) * nXBlockSize * nYBlockSize;
2059 values = rtalloc(valueslen);
2060 if (NULL == values) {
2061 rterror("rt_raster_to_gdal_mem: Could not allocate memory for GDAL band pixel values");
2062 if (allocBandNums) rtdealloc(bandNums);
2063 if (allocNodataValues) rtdealloc(excludeNodataValues);
2064 GDALClose(ds);
2065 return 0;
2066 }
2067
2068 for (iYBlock = 0; iYBlock < nYBlocks; iYBlock++) {
2069 for (iXBlock = 0; iXBlock < nXBlocks; iXBlock++) {
2070 memset(values, 0, valueslen);
2071
2072 x = iXBlock * nXBlockSize;
2073 y = iYBlock * nYBlockSize;
2074 RASTER_DEBUGF(4, "(iXBlock, iYBlock) = (%d, %d)", iXBlock, iYBlock);
2075 RASTER_DEBUGF(4, "(x, y) = (%d, %d)", x, y);
2076
2077 /* valid block width */
2078 if ((iXBlock + 1) * nXBlockSize > width)
2079 nXValid = width - (iXBlock * nXBlockSize);
2080 else
2081 nXValid = nXBlockSize;
2082
2083 /* valid block height */
2084 if ((iYBlock + 1) * nYBlockSize > height)
2085 nYValid = height - (iYBlock * nYBlockSize);
2086 else
2087 nYValid = nYBlockSize;
2088
2089 RASTER_DEBUGF(4, "(nXValid, nYValid) = (%d, %d)", nXValid, nYValid);
2090
2091 /* convert 8BSI values to 16BSI */
2092 z = 0;
2093 iYMax = y + nYValid;
2094 iXMax = x + nXValid;
2095 for (iY = y; iY < iYMax; iY++) {
2096 for (iX = x; iX < iXMax; iX++) {
2097 if (rt_band_get_pixel(rtband, iX, iY, &value, NULL) != ES_NONE) {
2098 rterror("rt_raster_to_gdal_mem: Could not get pixel value to convert from 8BSI to 16BSI");
2099 rtdealloc(values);
2100 if (allocBandNums) rtdealloc(bandNums);
2101 if (allocNodataValues) rtdealloc(excludeNodataValues);
2102 GDALClose(ds);
2103 return 0;
2104 }
2105
2106 values[z++] = rt_util_clamp_to_16BSI(value);
2107 }
2108 }
2109
2110 /* burn values */
2111 if (GDALRasterIO(
2112 band, GF_Write,
2113 x, y,
2114 nXValid, nYValid,
2115 values, nXValid, nYValid,
2116 gdal_pt,
2117 0, 0
2118 ) != CE_None) {
2119 rterror("rt_raster_to_gdal_mem: Could not write converted 8BSI to 16BSI values to GDAL band");
2120 rtdealloc(values);
2121 if (allocBandNums) rtdealloc(bandNums);
2122 if (allocNodataValues) rtdealloc(excludeNodataValues);
2123 GDALClose(ds);
2124 return 0;
2125 }
2126 }
2127 }
2128
2129 rtdealloc(values);
2130 }
2131
2132 /* Add nodata value for band */
2133 if (rt_band_get_hasnodata_flag(rtband) != FALSE && excludeNodataValues[i]) {
2134 rt_band_get_nodata(rtband, &nodata);
2135 if (GDALSetRasterNoDataValue(band, nodata) != CE_None)
2136 rtwarn("rt_raster_to_gdal_mem: Could not set nodata value for band");
2137 RASTER_DEBUGF(3, "nodata value set to %f", GDALGetRasterNoDataValue(band, NULL));
2138 }
2139
2140 #if POSTGIS_DEBUG_LEVEL > 3
2141 {
2142 GDALRasterBandH _grb = NULL;
2143 double _min;
2144 double _max;
2145 double _mean;
2146 double _stddev;
2147
2148 _grb = GDALGetRasterBand(ds, i + 1);
2149 GDALComputeRasterStatistics(_grb, FALSE, &_min, &_max, &_mean, &_stddev, NULL, NULL);
2150 RASTER_DEBUGF(4, "GDAL Band %d stats: %f, %f, %f, %f", i + 1, _min, _max, _mean, _stddev);
2151 }
2152 #endif
2153
2154 }
2155
2156 /* necessary??? */
2157 GDALFlushCache(ds);
2158
2159 if (allocBandNums) rtdealloc(bandNums);
2160 if (allocNodataValues) rtdealloc(excludeNodataValues);
2161
2162 return ds;
2163 }
2164
2165 /******************************************************************************
2166 * rt_raster_from_gdal_dataset()
2167 ******************************************************************************/
2168
2169 /**
2170 * Return a raster from a GDAL dataset
2171 *
2172 * @param ds : the GDAL dataset to convert to a raster
2173 *
2174 * @return raster or NULL
2175 */
2176 rt_raster
rt_raster_from_gdal_dataset(GDALDatasetH ds)2177 rt_raster_from_gdal_dataset(GDALDatasetH ds) {
2178 rt_raster rast = NULL;
2179 double gt[6] = {0};
2180 CPLErr cplerr;
2181 uint32_t width = 0;
2182 uint32_t height = 0;
2183 uint32_t numBands = 0;
2184 uint32_t i = 0;
2185 char *authname = NULL;
2186 char *authcode = NULL;
2187
2188 GDALRasterBandH gdband = NULL;
2189 GDALDataType gdpixtype = GDT_Unknown;
2190 rt_band band;
2191 int32_t idx;
2192 rt_pixtype pt = PT_END;
2193 uint32_t ptlen = 0;
2194 int hasnodata = 0;
2195 double nodataval;
2196
2197 int x;
2198 int y;
2199
2200 uint32_t nXBlocks, nYBlocks;
2201 int nXBlockSize, nYBlockSize;
2202 uint32_t iXBlock, iYBlock;
2203 uint32_t nXValid, nYValid;
2204 uint32_t iY;
2205
2206 uint8_t *values = NULL;
2207 uint32_t valueslen = 0;
2208 uint8_t *ptr = NULL;
2209
2210 assert(NULL != ds);
2211
2212 /* raster size */
2213 width = GDALGetRasterXSize(ds);
2214 height = GDALGetRasterYSize(ds);
2215 RASTER_DEBUGF(3, "Raster dimensions (width x height): %d x %d", width, height);
2216
2217 /* create new raster */
2218 RASTER_DEBUG(3, "Creating new raster");
2219 rast = rt_raster_new(width, height);
2220 if (NULL == rast) {
2221 rterror("rt_raster_from_gdal_dataset: Out of memory allocating new raster");
2222 return NULL;
2223 }
2224 RASTER_DEBUGF(3, "Created raster dimensions (width x height): %d x %d", rast->width, rast->height);
2225
2226 /* get raster attributes */
2227 cplerr = GDALGetGeoTransform(ds, gt);
2228 if (GDALGetGeoTransform(ds, gt) != CE_None) {
2229 RASTER_DEBUG(4, "Using default geotransform matrix (0, 1, 0, 0, 0, -1)");
2230 gt[0] = 0;
2231 gt[1] = 1;
2232 gt[2] = 0;
2233 gt[3] = 0;
2234 gt[4] = 0;
2235 gt[5] = -1;
2236 }
2237
2238 /* apply raster attributes */
2239 rt_raster_set_geotransform_matrix(rast, gt);
2240
2241 RASTER_DEBUGF(3, "Raster geotransform (%f, %f, %f, %f, %f, %f)",
2242 gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]);
2243
2244 /* srid */
2245 if (rt_util_gdal_sr_auth_info(ds, &authname, &authcode) == ES_NONE) {
2246 if (
2247 authname != NULL &&
2248 strcmp(authname, "EPSG") == 0 &&
2249 authcode != NULL
2250 ) {
2251 rt_raster_set_srid(rast, atoi(authcode));
2252 RASTER_DEBUGF(3, "New raster's SRID = %d", rast->srid);
2253 }
2254
2255 if (authname != NULL)
2256 rtdealloc(authname);
2257 if (authcode != NULL)
2258 rtdealloc(authcode);
2259 }
2260
2261 numBands = GDALGetRasterCount(ds);
2262
2263 #if POSTGIS_DEBUG_LEVEL > 3
2264 for (i = 1; i <= numBands; i++) {
2265 GDALRasterBandH _grb = NULL;
2266 double _min;
2267 double _max;
2268 double _mean;
2269 double _stddev;
2270
2271 _grb = GDALGetRasterBand(ds, i);
2272 GDALComputeRasterStatistics(_grb, FALSE, &_min, &_max, &_mean, &_stddev, NULL, NULL);
2273 RASTER_DEBUGF(4, "GDAL Band %d stats: %f, %f, %f, %f", i, _min, _max, _mean, _stddev);
2274 }
2275 #endif
2276
2277 /* copy bands */
2278 for (i = 1; i <= numBands; i++) {
2279 RASTER_DEBUGF(3, "Processing band %d of %d", i, numBands);
2280 gdband = NULL;
2281 gdband = GDALGetRasterBand(ds, i);
2282 if (NULL == gdband) {
2283 rterror("rt_raster_from_gdal_dataset: Could not get GDAL band");
2284 rt_raster_destroy(rast);
2285 return NULL;
2286 }
2287 RASTER_DEBUGF(4, "gdband @ %p", gdband);
2288
2289 /* pixtype */
2290 gdpixtype = GDALGetRasterDataType(gdband);
2291 RASTER_DEBUGF(4, "gdpixtype, size = %s, %d", GDALGetDataTypeName(gdpixtype), GDALGetDataTypeSize(gdpixtype) / 8);
2292 pt = rt_util_gdal_datatype_to_pixtype(gdpixtype);
2293 if (pt == PT_END) {
2294 rterror("rt_raster_from_gdal_dataset: Unknown pixel type for GDAL band");
2295 rt_raster_destroy(rast);
2296 return NULL;
2297 }
2298 ptlen = rt_pixtype_size(pt);
2299
2300 /* size: width and height */
2301 width = GDALGetRasterBandXSize(gdband);
2302 height = GDALGetRasterBandYSize(gdband);
2303 RASTER_DEBUGF(3, "GDAL band dimensions (width x height): %d x %d", width, height);
2304
2305 /* nodata */
2306 nodataval = GDALGetRasterNoDataValue(gdband, &hasnodata);
2307 RASTER_DEBUGF(3, "(hasnodata, nodataval) = (%d, %f)", hasnodata, nodataval);
2308
2309 /* create band object */
2310 idx = rt_raster_generate_new_band(
2311 rast, pt,
2312 (hasnodata ? nodataval : 0),
2313 hasnodata, nodataval, rt_raster_get_num_bands(rast)
2314 );
2315 if (idx < 0) {
2316 rterror("rt_raster_from_gdal_dataset: Could not allocate memory for raster band");
2317 rt_raster_destroy(rast);
2318 return NULL;
2319 }
2320 band = rt_raster_get_band(rast, idx);
2321 RASTER_DEBUGF(3, "Created band of dimension (width x height): %d x %d", band->width, band->height);
2322
2323 /* this makes use of GDAL's "natural" blocks */
2324 GDALGetBlockSize(gdband, &nXBlockSize, &nYBlockSize);
2325 nXBlocks = (width + nXBlockSize - 1) / nXBlockSize;
2326 nYBlocks = (height + nYBlockSize - 1) / nYBlockSize;
2327 RASTER_DEBUGF(4, "(nXBlockSize, nYBlockSize) = (%d, %d)", nXBlockSize, nYBlockSize);
2328 RASTER_DEBUGF(4, "(nXBlocks, nYBlocks) = (%d, %d)", nXBlocks, nYBlocks);
2329
2330 /* allocate memory for values */
2331 valueslen = ptlen * nXBlockSize * nYBlockSize;
2332 values = rtalloc(valueslen);
2333 if (values == NULL) {
2334 rterror("rt_raster_from_gdal_dataset: Could not allocate memory for GDAL band pixel values");
2335 rt_raster_destroy(rast);
2336 return NULL;
2337 }
2338 RASTER_DEBUGF(3, "values @ %p of length = %d", values, valueslen);
2339
2340 for (iYBlock = 0; iYBlock < nYBlocks; iYBlock++) {
2341 for (iXBlock = 0; iXBlock < nXBlocks; iXBlock++) {
2342 x = iXBlock * nXBlockSize;
2343 y = iYBlock * nYBlockSize;
2344 RASTER_DEBUGF(4, "(iXBlock, iYBlock) = (%d, %d)", iXBlock, iYBlock);
2345 RASTER_DEBUGF(4, "(x, y) = (%d, %d)", x, y);
2346
2347 memset(values, 0, valueslen);
2348
2349 /* valid block width */
2350 if ((iXBlock + 1) * nXBlockSize > width)
2351 nXValid = width - (iXBlock * nXBlockSize);
2352 else
2353 nXValid = nXBlockSize;
2354
2355 /* valid block height */
2356 if ((iYBlock + 1) * nYBlockSize > height)
2357 nYValid = height - (iYBlock * nYBlockSize);
2358 else
2359 nYValid = nYBlockSize;
2360
2361 RASTER_DEBUGF(4, "(nXValid, nYValid) = (%d, %d)", nXValid, nYValid);
2362
2363 cplerr = GDALRasterIO(
2364 gdband, GF_Read,
2365 x, y,
2366 nXValid, nYValid,
2367 values, nXValid, nYValid,
2368 gdpixtype,
2369 0, 0
2370 );
2371 if (cplerr != CE_None) {
2372 rterror("rt_raster_from_gdal_dataset: Could not get data from GDAL raster");
2373 rtdealloc(values);
2374 rt_raster_destroy(rast);
2375 return NULL;
2376 }
2377
2378 /* if block width is same as raster width, shortcut */
2379 if (nXBlocks == 1 && nYBlockSize > 1 && nXValid == width) {
2380 x = 0;
2381 y = nYBlockSize * iYBlock;
2382
2383 RASTER_DEBUGF(4, "Setting set of pixel lines at (%d, %d) for %d pixels", x, y, nXValid * nYValid);
2384 rt_band_set_pixel_line(band, x, y, values, nXValid * nYValid);
2385 }
2386 else {
2387 ptr = values;
2388 x = nXBlockSize * iXBlock;
2389 for (iY = 0; iY < nYValid; iY++) {
2390 y = iY + (nYBlockSize * iYBlock);
2391
2392 RASTER_DEBUGF(4, "Setting pixel line at (%d, %d) for %d pixels", x, y, nXValid);
2393 rt_band_set_pixel_line(band, x, y, ptr, nXValid);
2394 ptr += (nXValid * ptlen);
2395 }
2396 }
2397 }
2398 }
2399
2400 /* free memory */
2401 rtdealloc(values);
2402 }
2403
2404 return rast;
2405 }
2406
2407 /******************************************************************************
2408 * rt_raster_gdal_rasterize()
2409 ******************************************************************************/
2410
2411 typedef struct _rti_rasterize_arg_t* _rti_rasterize_arg;
2412 struct _rti_rasterize_arg_t {
2413 uint8_t noband;
2414
2415 uint32_t numbands;
2416
2417 OGRSpatialReferenceH src_sr;
2418
2419 rt_pixtype *pixtype;
2420 double *init;
2421 double *nodata;
2422 uint8_t *hasnodata;
2423 double *value;
2424 int *bandlist;
2425 };
2426
2427 static _rti_rasterize_arg
_rti_rasterize_arg_init()2428 _rti_rasterize_arg_init() {
2429 _rti_rasterize_arg arg = NULL;
2430
2431 arg = rtalloc(sizeof(struct _rti_rasterize_arg_t));
2432 if (arg == NULL) {
2433 rterror("_rti_rasterize_arg_init: Could not allocate memory for _rti_rasterize_arg");
2434 return NULL;
2435 }
2436
2437 arg->noband = 0;
2438
2439 arg->numbands = 0;
2440
2441 arg->src_sr = NULL;
2442
2443 arg->pixtype = NULL;
2444 arg->init = NULL;
2445 arg->nodata = NULL;
2446 arg->hasnodata = NULL;
2447 arg->value = NULL;
2448 arg->bandlist = NULL;
2449
2450 return arg;
2451 }
2452
2453 static void
_rti_rasterize_arg_destroy(_rti_rasterize_arg arg)2454 _rti_rasterize_arg_destroy(_rti_rasterize_arg arg) {
2455 if (arg->noband) {
2456 if (arg->pixtype != NULL)
2457 rtdealloc(arg->pixtype);
2458 if (arg->init != NULL)
2459 rtdealloc(arg->init);
2460 if (arg->nodata != NULL)
2461 rtdealloc(arg->nodata);
2462 if (arg->hasnodata != NULL)
2463 rtdealloc(arg->hasnodata);
2464 if (arg->value != NULL)
2465 rtdealloc(arg->value);
2466 }
2467
2468 if (arg->bandlist != NULL)
2469 rtdealloc(arg->bandlist);
2470
2471 if (arg->src_sr != NULL)
2472 OSRDestroySpatialReference(arg->src_sr);
2473
2474 rtdealloc(arg);
2475 }
2476
2477 /**
2478 * Return a raster of the provided geometry
2479 *
2480 * @param wkb : WKB representation of the geometry to convert
2481 * @param wkb_len : length of the WKB representation of the geometry
2482 * @param srs : the geometry's coordinate system in OGC WKT
2483 * @param num_bands : number of bands in the output raster
2484 * @param pixtype : data type of each band
2485 * @param init : array of values to initialize each band with
2486 * @param value : array of values for pixels of geometry
2487 * @param nodata : array of nodata values for each band
2488 * @param hasnodata : array flagging the presence of nodata for each band
2489 * @param width : the number of columns of the raster
2490 * @param height : the number of rows of the raster
2491 * @param scale_x : the pixel width of the raster
2492 * @param scale_y : the pixel height of the raster
2493 * @param ul_xw : the X value of upper-left corner of the raster
2494 * @param ul_yw : the Y value of upper-left corner of the raster
2495 * @param grid_xw : the X value of point on grid to align raster to
2496 * @param grid_yw : the Y value of point on grid to align raster to
2497 * @param skew_x : the X skew of the raster
2498 * @param skew_y : the Y skew of the raster
2499 * @param options : array of options. only option is "ALL_TOUCHED"
2500 *
2501 * @return the raster of the provided geometry or NULL
2502 */
2503 rt_raster
rt_raster_gdal_rasterize(const unsigned char * wkb,uint32_t wkb_len,const char * srs,uint32_t num_bands,rt_pixtype * pixtype,double * init,double * value,double * nodata,uint8_t * hasnodata,int * width,int * height,double * scale_x,double * scale_y,double * ul_xw,double * ul_yw,double * grid_xw,double * grid_yw,double * skew_x,double * skew_y,char ** options)2504 rt_raster_gdal_rasterize(
2505 const unsigned char *wkb, uint32_t wkb_len,
2506 const char *srs,
2507 uint32_t num_bands, rt_pixtype *pixtype,
2508 double *init, double *value,
2509 double *nodata, uint8_t *hasnodata,
2510 int *width, int *height,
2511 double *scale_x, double *scale_y,
2512 double *ul_xw, double *ul_yw,
2513 double *grid_xw, double *grid_yw,
2514 double *skew_x, double *skew_y,
2515 char **options
2516 ) {
2517 rt_raster rast = NULL;
2518 uint32_t i = 0;
2519 int err = 0;
2520
2521 _rti_rasterize_arg arg = NULL;
2522
2523 int _dim[2] = {0};
2524 double _scale[2] = {0};
2525 double _skew[2] = {0};
2526
2527 OGRErr ogrerr;
2528 OGRGeometryH src_geom;
2529 OGREnvelope src_env;
2530 rt_envelope extent;
2531 OGRwkbGeometryType wkbtype = wkbUnknown;
2532
2533 int ul_user = 0;
2534
2535 CPLErr cplerr;
2536 double _gt[6] = {0};
2537 GDALDriverH _drv = NULL;
2538 int unload_drv = 0;
2539 GDALDatasetH _ds = NULL;
2540 GDALRasterBandH _band = NULL;
2541
2542 uint16_t _width = 0;
2543 uint16_t _height = 0;
2544
2545 RASTER_DEBUG(3, "starting");
2546
2547 assert(NULL != wkb);
2548 assert(0 != wkb_len);
2549
2550 /* internal variables */
2551 arg = _rti_rasterize_arg_init();
2552 if (arg == NULL) {
2553 rterror("rt_raster_gdal_rasterize: Could not initialize internal variables");
2554 return NULL;
2555 }
2556
2557 /* no bands, raster is a mask */
2558 if (num_bands < 1) {
2559 arg->noband = 1;
2560 arg->numbands = 1;
2561
2562 arg->pixtype = (rt_pixtype *) rtalloc(sizeof(rt_pixtype));
2563 arg->pixtype[0] = PT_8BUI;
2564
2565 arg->init = (double *) rtalloc(sizeof(double));
2566 arg->init[0] = 0;
2567
2568 arg->nodata = (double *) rtalloc(sizeof(double));
2569 arg->nodata[0] = 0;
2570
2571 arg->hasnodata = (uint8_t *) rtalloc(sizeof(uint8_t));
2572 arg->hasnodata[0] = 1;
2573
2574 arg->value = (double *) rtalloc(sizeof(double));
2575 arg->value[0] = 1;
2576 }
2577 else {
2578 arg->noband = 0;
2579 arg->numbands = num_bands;
2580
2581 arg->pixtype = pixtype;
2582 arg->init = init;
2583 arg->nodata = nodata;
2584 arg->hasnodata = hasnodata;
2585 arg->value = value;
2586 }
2587
2588 /* OGR spatial reference */
2589 if (NULL != srs && strlen(srs)) {
2590 arg->src_sr = OSRNewSpatialReference(NULL);
2591 if (OSRSetFromUserInput(arg->src_sr, srs) != OGRERR_NONE) {
2592 rterror("rt_raster_gdal_rasterize: Could not create OSR spatial reference using the provided srs: %s", srs);
2593 _rti_rasterize_arg_destroy(arg);
2594 return NULL;
2595 }
2596 }
2597
2598 /* convert WKB to OGR Geometry */
2599 ogrerr = OGR_G_CreateFromWkb((unsigned char *) wkb, arg->src_sr, &src_geom, wkb_len);
2600 if (ogrerr != OGRERR_NONE) {
2601 rterror("rt_raster_gdal_rasterize: Could not create OGR Geometry from WKB");
2602
2603 _rti_rasterize_arg_destroy(arg);
2604 /* OGRCleanupAll(); */
2605
2606 return NULL;
2607 }
2608
2609 /* OGR Geometry is empty */
2610 if (OGR_G_IsEmpty(src_geom)) {
2611 rtinfo("Geometry provided is empty. Returning empty raster");
2612
2613 OGR_G_DestroyGeometry(src_geom);
2614 _rti_rasterize_arg_destroy(arg);
2615 /* OGRCleanupAll(); */
2616
2617 return rt_raster_new(0, 0);
2618 }
2619
2620 /* get envelope */
2621 OGR_G_GetEnvelope(src_geom, &src_env);
2622 rt_util_from_ogr_envelope(src_env, &extent);
2623
2624 RASTER_DEBUGF(3, "Suggested raster envelope: %f, %f, %f, %f",
2625 extent.MinX, extent.MinY, extent.MaxX, extent.MaxY);
2626
2627 /* user-defined scale */
2628 if (
2629 (NULL != scale_x) &&
2630 (NULL != scale_y) &&
2631 (FLT_NEQ(*scale_x, 0.0)) &&
2632 (FLT_NEQ(*scale_y, 0.0))
2633 ) {
2634 /* for now, force scale to be in left-right, top-down orientation */
2635 _scale[0] = fabs(*scale_x);
2636 _scale[1] = fabs(*scale_y);
2637 }
2638 /* user-defined width/height */
2639 else if ((NULL != width) && (NULL != height) && (*width != 0) && (*height != 0))
2640 {
2641 _dim[0] = abs(*width);
2642 _dim[1] = abs(*height);
2643
2644 if (FLT_NEQ(extent.MaxX, extent.MinX))
2645 _scale[0] = fabs((extent.MaxX - extent.MinX) / _dim[0]);
2646 else
2647 _scale[0] = 1.;
2648
2649 if (FLT_NEQ(extent.MaxY, extent.MinY))
2650 _scale[1] = fabs((extent.MaxY - extent.MinY) / _dim[1]);
2651 else
2652 _scale[1] = 1.;
2653 }
2654 else {
2655 rterror("rt_raster_gdal_rasterize: Values must be provided for width and height or X and Y of scale");
2656
2657 OGR_G_DestroyGeometry(src_geom);
2658 _rti_rasterize_arg_destroy(arg);
2659 /* OGRCleanupAll(); */
2660
2661 return NULL;
2662 }
2663 RASTER_DEBUGF(3, "scale (x, y) = %f, %f", _scale[0], -1 * _scale[1]);
2664 RASTER_DEBUGF(3, "dim (x, y) = %d, %d", _dim[0], _dim[1]);
2665
2666 /* user-defined skew */
2667 if (NULL != skew_x) {
2668 _skew[0] = *skew_x;
2669
2670 /*
2671 negative scale-x affects skew
2672 for now, force skew to be in left-right, top-down orientation
2673 */
2674 if (
2675 NULL != scale_x &&
2676 *scale_x < 0.
2677 ) {
2678 _skew[0] *= -1;
2679 }
2680 }
2681 if (NULL != skew_y) {
2682 _skew[1] = *skew_y;
2683
2684 /*
2685 positive scale-y affects skew
2686 for now, force skew to be in left-right, top-down orientation
2687 */
2688 if (
2689 NULL != scale_y &&
2690 *scale_y > 0.
2691 ) {
2692 _skew[1] *= -1;
2693 }
2694 }
2695
2696 /*
2697 if geometry is a point, a linestring or set of either and bounds not set,
2698 increase extent by a pixel to avoid missing points on border
2699
2700 a whole pixel is used instead of half-pixel due to backward
2701 compatibility with GDAL 1.6, 1.7 and 1.8. 1.9+ works fine with half-pixel.
2702 */
2703 wkbtype = wkbFlatten(OGR_G_GetGeometryType(src_geom));
2704 if ((
2705 (wkbtype == wkbPoint) ||
2706 (wkbtype == wkbMultiPoint) ||
2707 (wkbtype == wkbLineString) ||
2708 (wkbtype == wkbMultiLineString)
2709 ) &&
2710 _dim[0] == 0 &&
2711 _dim[1] == 0
2712 ) {
2713 int result;
2714 LWPOLY *epoly = NULL;
2715 LWGEOM *lwgeom = NULL;
2716 GEOSGeometry *egeom = NULL;
2717 GEOSGeometry *geom = NULL;
2718
2719 RASTER_DEBUG(3, "Testing geometry is properly contained by extent");
2720
2721 /*
2722 see if geometry is properly contained by extent
2723 all parts of geometry lies within extent
2724 */
2725
2726 /* initialize GEOS */
2727 initGEOS(rtinfo, lwgeom_geos_error);
2728
2729 /* convert envelope to geometry */
2730 RASTER_DEBUG(4, "Converting envelope to geometry");
2731 epoly = rt_util_envelope_to_lwpoly(extent);
2732 if (epoly == NULL) {
2733 rterror("rt_raster_gdal_rasterize: Could not create envelope's geometry to test if geometry is properly contained by extent");
2734
2735 OGR_G_DestroyGeometry(src_geom);
2736 _rti_rasterize_arg_destroy(arg);
2737 /* OGRCleanupAll(); */
2738
2739 return NULL;
2740 }
2741
2742 egeom = (GEOSGeometry *) LWGEOM2GEOS(lwpoly_as_lwgeom(epoly), 0);
2743 lwpoly_free(epoly);
2744
2745 /* convert WKB to geometry */
2746 RASTER_DEBUG(4, "Converting WKB to geometry");
2747 lwgeom = lwgeom_from_wkb(wkb, wkb_len, LW_PARSER_CHECK_NONE);
2748 geom = (GEOSGeometry *) LWGEOM2GEOS(lwgeom, 0);
2749 lwgeom_free(lwgeom);
2750
2751 result = GEOSRelatePattern(egeom, geom, "T**FF*FF*");
2752 GEOSGeom_destroy(geom);
2753 GEOSGeom_destroy(egeom);
2754
2755 if (result == 2) {
2756 rterror("rt_raster_gdal_rasterize: Could not test if geometry is properly contained by extent for geometry within extent");
2757
2758 OGR_G_DestroyGeometry(src_geom);
2759 _rti_rasterize_arg_destroy(arg);
2760 /* OGRCleanupAll(); */
2761
2762 return NULL;
2763 }
2764
2765 /* geometry NOT properly contained by extent */
2766 if (!result) {
2767
2768 #if POSTGIS_GDAL_VERSION > 18
2769
2770 /* check alignment flag: grid_xw */
2771 if (
2772 (NULL == ul_xw && NULL == ul_yw) &&
2773 (NULL != grid_xw && NULL != grid_yw) &&
2774 FLT_NEQ(*grid_xw, extent.MinX)
2775 ) {
2776 /* do nothing */
2777 RASTER_DEBUG(3, "Skipping extent adjustment on X-axis due to upcoming alignment");
2778 }
2779 else {
2780 RASTER_DEBUG(3, "Adjusting extent for GDAL > 1.8 by half the scale on X-axis");
2781 extent.MinX -= (_scale[0] / 2.);
2782 extent.MaxX += (_scale[0] / 2.);
2783 }
2784
2785 /* check alignment flag: grid_yw */
2786 if (
2787 (NULL == ul_xw && NULL == ul_yw) &&
2788 (NULL != grid_xw && NULL != grid_yw) &&
2789 FLT_NEQ(*grid_yw, extent.MaxY)
2790 ) {
2791 /* do nothing */
2792 RASTER_DEBUG(3, "Skipping extent adjustment on Y-axis due to upcoming alignment");
2793 }
2794 else {
2795 RASTER_DEBUG(3, "Adjusting extent for GDAL > 1.8 by half the scale on Y-axis");
2796 extent.MinY -= (_scale[1] / 2.);
2797 extent.MaxY += (_scale[1] / 2.);
2798 }
2799
2800 #else
2801
2802 /* check alignment flag: grid_xw */
2803 if (
2804 (NULL == ul_xw && NULL == ul_yw) &&
2805 (NULL != grid_xw && NULL != grid_yw) &&
2806 FLT_NEQ(*grid_xw, extent.MinX)
2807 ) {
2808 /* do nothing */
2809 RASTER_DEBUG(3, "Skipping extent adjustment on X-axis due to upcoming alignment");
2810 }
2811 else {
2812 RASTER_DEBUG(3, "Adjusting extent for GDAL <= 1.8 by the scale on X-axis");
2813 extent.MinX -= _scale[0];
2814 extent.MaxX += _scale[0];
2815 }
2816
2817
2818 /* check alignment flag: grid_yw */
2819 if (
2820 (NULL == ul_xw && NULL == ul_yw) &&
2821 (NULL != grid_xw && NULL != grid_yw) &&
2822 FLT_NEQ(*grid_yw, extent.MaxY)
2823 ) {
2824 /* do nothing */
2825 RASTER_DEBUG(3, "Skipping extent adjustment on Y-axis due to upcoming alignment");
2826 }
2827 else {
2828 RASTER_DEBUG(3, "Adjusting extent for GDAL <= 1.8 by the scale on Y-axis");
2829 extent.MinY -= _scale[1];
2830 extent.MaxY += _scale[1];
2831 }
2832
2833 #endif
2834
2835 }
2836
2837 RASTER_DEBUGF(3, "Adjusted extent: %f, %f, %f, %f",
2838 extent.MinX, extent.MinY, extent.MaxX, extent.MaxY);
2839
2840 extent.UpperLeftX = extent.MinX;
2841 extent.UpperLeftY = extent.MaxY;
2842 }
2843
2844 /* reprocess extent if skewed */
2845 if (FLT_NEQ(_skew[0], 0.0) || FLT_NEQ(_skew[1], 0.0))
2846 {
2847 rt_raster skewedrast;
2848
2849 RASTER_DEBUG(3, "Computing skewed extent's envelope");
2850
2851 skewedrast = rt_raster_compute_skewed_raster(
2852 extent,
2853 _skew,
2854 _scale,
2855 0.01
2856 );
2857 if (skewedrast == NULL) {
2858 rterror("rt_raster_gdal_rasterize: Could not compute skewed raster");
2859
2860 OGR_G_DestroyGeometry(src_geom);
2861 _rti_rasterize_arg_destroy(arg);
2862 /* OGRCleanupAll(); */
2863
2864 return NULL;
2865 }
2866
2867 _dim[0] = skewedrast->width;
2868 _dim[1] = skewedrast->height;
2869
2870 extent.UpperLeftX = skewedrast->ipX;
2871 extent.UpperLeftY = skewedrast->ipY;
2872
2873 rt_raster_destroy(skewedrast);
2874 }
2875
2876 /* raster dimensions */
2877 if (!_dim[0])
2878 _dim[0] = (int) fmax((fabs(extent.MaxX - extent.MinX) + (_scale[0] / 2.)) / _scale[0], 1);
2879 if (!_dim[1])
2880 _dim[1] = (int) fmax((fabs(extent.MaxY - extent.MinY) + (_scale[1] / 2.)) / _scale[1], 1);
2881
2882 /* temporary raster */
2883 rast = rt_raster_new(_dim[0], _dim[1]);
2884 if (rast == NULL) {
2885 rterror("rt_raster_gdal_rasterize: Out of memory allocating temporary raster");
2886
2887 OGR_G_DestroyGeometry(src_geom);
2888 _rti_rasterize_arg_destroy(arg);
2889 /* OGRCleanupAll(); */
2890
2891 return NULL;
2892 }
2893
2894 /* set raster's spatial attributes */
2895 rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
2896 rt_raster_set_scale(rast, _scale[0], -1 * _scale[1]);
2897 rt_raster_set_skews(rast, _skew[0], _skew[1]);
2898
2899 rt_raster_get_geotransform_matrix(rast, _gt);
2900 RASTER_DEBUGF(3, "Temp raster's geotransform: %f, %f, %f, %f, %f, %f",
2901 _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
2902 RASTER_DEBUGF(3, "Temp raster's dimensions (width x height): %d x %d",
2903 _dim[0], _dim[1]);
2904
2905 /* user-specified upper-left corner */
2906 if (
2907 NULL != ul_xw &&
2908 NULL != ul_yw
2909 ) {
2910 ul_user = 1;
2911
2912 RASTER_DEBUGF(4, "Using user-specified upper-left corner: %f, %f", *ul_xw, *ul_yw);
2913
2914 /* set upper-left corner */
2915 rt_raster_set_offsets(rast, *ul_xw, *ul_yw);
2916 extent.UpperLeftX = *ul_xw;
2917 extent.UpperLeftY = *ul_yw;
2918 }
2919 else if (
2920 ((NULL != ul_xw) && (NULL == ul_yw)) ||
2921 ((NULL == ul_xw) && (NULL != ul_yw))
2922 ) {
2923 rterror("rt_raster_gdal_rasterize: Both X and Y upper-left corner values must be provided");
2924
2925 rt_raster_destroy(rast);
2926 OGR_G_DestroyGeometry(src_geom);
2927 _rti_rasterize_arg_destroy(arg);
2928 /* OGRCleanupAll(); */
2929
2930 return NULL;
2931 }
2932
2933 /* alignment only considered if upper-left corner not provided */
2934 if (
2935 !ul_user && (
2936 (NULL != grid_xw) || (NULL != grid_yw)
2937 )
2938 ) {
2939
2940 if (
2941 ((NULL != grid_xw) && (NULL == grid_yw)) ||
2942 ((NULL == grid_xw) && (NULL != grid_yw))
2943 ) {
2944 rterror("rt_raster_gdal_rasterize: Both X and Y alignment values must be provided");
2945
2946 rt_raster_destroy(rast);
2947 OGR_G_DestroyGeometry(src_geom);
2948 _rti_rasterize_arg_destroy(arg);
2949 /* OGRCleanupAll(); */
2950
2951 return NULL;
2952 }
2953
2954 RASTER_DEBUGF(4, "Aligning extent to user-specified grid: %f, %f", *grid_xw, *grid_yw);
2955
2956 do {
2957 double _r[2] = {0};
2958 double _w[2] = {0};
2959
2960 /* raster is already aligned */
2961 if (FLT_EQ(*grid_xw, extent.UpperLeftX) && FLT_EQ(*grid_yw, extent.UpperLeftY)) {
2962 RASTER_DEBUG(3, "Skipping raster alignment as it is already aligned to grid");
2963 break;
2964 }
2965
2966 extent.UpperLeftX = rast->ipX;
2967 extent.UpperLeftY = rast->ipY;
2968 rt_raster_set_offsets(rast, *grid_xw, *grid_yw);
2969
2970 /* process upper-left corner */
2971 if (rt_raster_geopoint_to_cell(
2972 rast,
2973 extent.UpperLeftX, extent.UpperLeftY,
2974 &(_r[0]), &(_r[1]),
2975 NULL
2976 ) != ES_NONE) {
2977 rterror("rt_raster_gdal_rasterize: Could not compute raster pixel for spatial coordinates");
2978
2979 rt_raster_destroy(rast);
2980 OGR_G_DestroyGeometry(src_geom);
2981 _rti_rasterize_arg_destroy(arg);
2982 /* OGRCleanupAll(); */
2983
2984 return NULL;
2985 }
2986
2987 if (rt_raster_cell_to_geopoint(
2988 rast,
2989 _r[0], _r[1],
2990 &(_w[0]), &(_w[1]),
2991 NULL
2992 ) != ES_NONE) {
2993 rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
2994
2995 rt_raster_destroy(rast);
2996 OGR_G_DestroyGeometry(src_geom);
2997 _rti_rasterize_arg_destroy(arg);
2998 /* OGRCleanupAll(); */
2999
3000 return NULL;
3001 }
3002
3003 /* shift occurred */
3004 if (FLT_NEQ(_w[0], extent.UpperLeftX)) {
3005 if (NULL == width)
3006 rast->width++;
3007 else if (NULL == scale_x) {
3008 double _c[2] = {0};
3009
3010 rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
3011
3012 /* get upper-right corner */
3013 if (rt_raster_cell_to_geopoint(
3014 rast,
3015 rast->width, 0,
3016 &(_c[0]), &(_c[1]),
3017 NULL
3018 ) != ES_NONE) {
3019 rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
3020
3021 rt_raster_destroy(rast);
3022 OGR_G_DestroyGeometry(src_geom);
3023 _rti_rasterize_arg_destroy(arg);
3024 /* OGRCleanupAll(); */
3025
3026 return NULL;
3027 }
3028
3029 rast->scaleX = fabs((_c[0] - _w[0]) / ((double) rast->width));
3030 }
3031 }
3032 if (FLT_NEQ(_w[1], extent.UpperLeftY)) {
3033 if (NULL == height)
3034 rast->height++;
3035 else if (NULL == scale_y) {
3036 double _c[2] = {0};
3037
3038 rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
3039
3040 /* get upper-right corner */
3041 if (rt_raster_cell_to_geopoint(
3042 rast,
3043 0, rast->height,
3044 &(_c[0]), &(_c[1]),
3045 NULL
3046 ) != ES_NONE) {
3047 rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
3048
3049 rt_raster_destroy(rast);
3050 OGR_G_DestroyGeometry(src_geom);
3051 _rti_rasterize_arg_destroy(arg);
3052 /* OGRCleanupAll(); */
3053
3054 return NULL;
3055 }
3056
3057 rast->scaleY = -1 * fabs((_c[1] - _w[1]) / ((double) rast->height));
3058 }
3059 }
3060
3061 rt_raster_set_offsets(rast, _w[0], _w[1]);
3062 }
3063 while (0);
3064 }
3065
3066 /*
3067 after this point, rt_envelope extent is no longer used
3068 */
3069
3070 /* get key attributes from rast */
3071 _dim[0] = rast->width;
3072 _dim[1] = rast->height;
3073 rt_raster_get_geotransform_matrix(rast, _gt);
3074
3075 /* scale-x is negative or scale-y is positive */
3076 if ((
3077 (NULL != scale_x) && (*scale_x < 0.)
3078 ) || (
3079 (NULL != scale_y) && (*scale_y > 0)
3080 )) {
3081 double _w[2] = {0};
3082
3083 /* negative scale-x */
3084 if (
3085 (NULL != scale_x) &&
3086 (*scale_x < 0.)
3087 ) {
3088 RASTER_DEBUG(3, "Processing negative scale-x");
3089
3090 if (rt_raster_cell_to_geopoint(
3091 rast,
3092 _dim[0], 0,
3093 &(_w[0]), &(_w[1]),
3094 NULL
3095 ) != ES_NONE) {
3096 rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
3097
3098 rt_raster_destroy(rast);
3099 OGR_G_DestroyGeometry(src_geom);
3100 _rti_rasterize_arg_destroy(arg);
3101 /* OGRCleanupAll(); */
3102
3103 return NULL;
3104 }
3105
3106 _gt[0] = _w[0];
3107 _gt[1] = *scale_x;
3108
3109 /* check for skew */
3110 if (NULL != skew_x && FLT_NEQ(*skew_x, 0.0))
3111 _gt[2] = *skew_x;
3112 }
3113 /* positive scale-y */
3114 if (
3115 (NULL != scale_y) &&
3116 (*scale_y > 0)
3117 ) {
3118 RASTER_DEBUG(3, "Processing positive scale-y");
3119
3120 if (rt_raster_cell_to_geopoint(
3121 rast,
3122 0, _dim[1],
3123 &(_w[0]), &(_w[1]),
3124 NULL
3125 ) != ES_NONE) {
3126 rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
3127
3128 rt_raster_destroy(rast);
3129 OGR_G_DestroyGeometry(src_geom);
3130 _rti_rasterize_arg_destroy(arg);
3131 /* OGRCleanupAll(); */
3132
3133 return NULL;
3134 }
3135
3136 _gt[3] = _w[1];
3137 _gt[5] = *scale_y;
3138
3139 /* check for skew */
3140 if (NULL != skew_y && FLT_NEQ(*skew_y, 0.0))
3141 _gt[4] = *skew_y;
3142 }
3143 }
3144
3145 rt_raster_destroy(rast);
3146 rast = NULL;
3147
3148 RASTER_DEBUGF(3, "Applied geotransform: %f, %f, %f, %f, %f, %f",
3149 _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
3150 RASTER_DEBUGF(3, "Raster dimensions (width x height): %d x %d",
3151 _dim[0], _dim[1]);
3152
3153 /* load GDAL mem */
3154 if (!rt_util_gdal_driver_registered("MEM")) {
3155 RASTER_DEBUG(4, "Registering MEM driver");
3156 GDALRegister_MEM();
3157 unload_drv = 1;
3158 }
3159 _drv = GDALGetDriverByName("MEM");
3160 if (NULL == _drv) {
3161 rterror("rt_raster_gdal_rasterize: Could not load the MEM GDAL driver");
3162
3163 OGR_G_DestroyGeometry(src_geom);
3164 _rti_rasterize_arg_destroy(arg);
3165 /* OGRCleanupAll(); */
3166
3167 return NULL;
3168 }
3169
3170 /* unload driver from GDAL driver manager */
3171 if (unload_drv) {
3172 RASTER_DEBUG(4, "Deregistering MEM driver");
3173 GDALDeregisterDriver(_drv);
3174 }
3175
3176 _ds = GDALCreate(_drv, "", _dim[0], _dim[1], 0, GDT_Byte, NULL);
3177 if (NULL == _ds) {
3178 rterror("rt_raster_gdal_rasterize: Could not create a GDALDataset to rasterize the geometry into");
3179
3180 OGR_G_DestroyGeometry(src_geom);
3181 _rti_rasterize_arg_destroy(arg);
3182 /* OGRCleanupAll(); */
3183 if (unload_drv) GDALDestroyDriver(_drv);
3184
3185 return NULL;
3186 }
3187
3188 /* set geotransform */
3189 cplerr = GDALSetGeoTransform(_ds, _gt);
3190 if (cplerr != CE_None) {
3191 rterror("rt_raster_gdal_rasterize: Could not set geotransform on GDALDataset");
3192
3193 OGR_G_DestroyGeometry(src_geom);
3194 _rti_rasterize_arg_destroy(arg);
3195 /* OGRCleanupAll(); */
3196
3197 GDALClose(_ds);
3198 if (unload_drv) GDALDestroyDriver(_drv);
3199
3200 return NULL;
3201 }
3202
3203 /* set SRS */
3204 if (NULL != arg->src_sr) {
3205 char *_srs = NULL;
3206 OSRExportToWkt(arg->src_sr, &_srs);
3207
3208 cplerr = GDALSetProjection(_ds, _srs);
3209 CPLFree(_srs);
3210 if (cplerr != CE_None) {
3211 rterror("rt_raster_gdal_rasterize: Could not set projection on GDALDataset");
3212
3213 OGR_G_DestroyGeometry(src_geom);
3214 _rti_rasterize_arg_destroy(arg);
3215 /* OGRCleanupAll(); */
3216
3217 GDALClose(_ds);
3218 if (unload_drv) GDALDestroyDriver(_drv);
3219
3220 return NULL;
3221 }
3222 }
3223
3224 /* set bands */
3225 for (i = 0; i < arg->numbands; i++) {
3226 err = 0;
3227
3228 do {
3229 /* add band */
3230 cplerr = GDALAddBand(_ds, rt_util_pixtype_to_gdal_datatype(arg->pixtype[i]), NULL);
3231 if (cplerr != CE_None) {
3232 rterror("rt_raster_gdal_rasterize: Could not add band to GDALDataset");
3233 err = 1;
3234 break;
3235 }
3236
3237 _band = GDALGetRasterBand(_ds, i + 1);
3238 if (NULL == _band) {
3239 rterror("rt_raster_gdal_rasterize: Could not get band %d from GDALDataset", i + 1);
3240 err = 1;
3241 break;
3242 }
3243
3244 /* nodata value */
3245 if (arg->hasnodata[i]) {
3246 RASTER_DEBUGF(4, "Setting NODATA value of band %d to %f", i, arg->nodata[i]);
3247 cplerr = GDALSetRasterNoDataValue(_band, arg->nodata[i]);
3248 if (cplerr != CE_None) {
3249 rterror("rt_raster_gdal_rasterize: Could not set nodata value");
3250 err = 1;
3251 break;
3252 }
3253 RASTER_DEBUGF(4, "NODATA value set to %f", GDALGetRasterNoDataValue(_band, NULL));
3254 }
3255
3256 /* initial value */
3257 cplerr = GDALFillRaster(_band, arg->init[i], 0);
3258 if (cplerr != CE_None) {
3259 rterror("rt_raster_gdal_rasterize: Could not set initial value");
3260 err = 1;
3261 break;
3262 }
3263 }
3264 while (0);
3265
3266 if (err) {
3267
3268 OGR_G_DestroyGeometry(src_geom);
3269 _rti_rasterize_arg_destroy(arg);
3270
3271 /* OGRCleanupAll(); */
3272
3273 GDALClose(_ds);
3274 if (unload_drv) GDALDestroyDriver(_drv);
3275
3276 return NULL;
3277 }
3278 }
3279
3280 arg->bandlist = (int *) rtalloc(sizeof(int) * arg->numbands);
3281 for (i = 0; i < arg->numbands; i++) arg->bandlist[i] = i + 1;
3282
3283 /* burn geometry */
3284 cplerr = GDALRasterizeGeometries(
3285 _ds,
3286 arg->numbands, arg->bandlist,
3287 1, &src_geom,
3288 NULL, NULL,
3289 arg->value,
3290 options,
3291 NULL, NULL
3292 );
3293 if (cplerr != CE_None) {
3294 rterror("rt_raster_gdal_rasterize: Could not rasterize geometry");
3295
3296 OGR_G_DestroyGeometry(src_geom);
3297 _rti_rasterize_arg_destroy(arg);
3298 /* OGRCleanupAll(); */
3299
3300 GDALClose(_ds);
3301 if (unload_drv) GDALDestroyDriver(_drv);
3302
3303 return NULL;
3304 }
3305
3306 /* convert gdal dataset to raster */
3307 GDALFlushCache(_ds);
3308 RASTER_DEBUG(3, "Converting GDAL dataset to raster");
3309 rast = rt_raster_from_gdal_dataset(_ds);
3310
3311 OGR_G_DestroyGeometry(src_geom);
3312 /* OGRCleanupAll(); */
3313
3314 GDALClose(_ds);
3315 if (unload_drv) GDALDestroyDriver(_drv);
3316
3317 if (NULL == rast) {
3318 rterror("rt_raster_gdal_rasterize: Could not rasterize geometry");
3319 return NULL;
3320 }
3321
3322 /* width, height */
3323 _width = rt_raster_get_width(rast);
3324 _height = rt_raster_get_height(rast);
3325
3326 /* check each band for pixtype */
3327 for (i = 0; i < arg->numbands; i++) {
3328 uint8_t *data = NULL;
3329 rt_band band = NULL;
3330 rt_band oldband = NULL;
3331
3332 double val = 0;
3333 int nodata = 0;
3334 int hasnodata = 0;
3335 double nodataval = 0;
3336 int x = 0;
3337 int y = 0;
3338
3339 oldband = rt_raster_get_band(rast, i);
3340 if (oldband == NULL) {
3341 rterror("rt_raster_gdal_rasterize: Could not get band %d of output raster", i);
3342 _rti_rasterize_arg_destroy(arg);
3343 rt_raster_destroy(rast);
3344 return NULL;
3345 }
3346
3347 /* band is of user-specified type */
3348 if (rt_band_get_pixtype(oldband) == arg->pixtype[i])
3349 continue;
3350
3351 /* hasnodata, nodataval */
3352 hasnodata = rt_band_get_hasnodata_flag(oldband);
3353 if (hasnodata)
3354 rt_band_get_nodata(oldband, &nodataval);
3355
3356 /* allocate data */
3357 data = rtalloc(rt_pixtype_size(arg->pixtype[i]) * _width * _height);
3358 if (data == NULL) {
3359 rterror("rt_raster_gdal_rasterize: Could not allocate memory for band data");
3360 _rti_rasterize_arg_destroy(arg);
3361 rt_raster_destroy(rast);
3362 return NULL;
3363 }
3364 memset(data, 0, rt_pixtype_size(arg->pixtype[i]) * _width * _height);
3365
3366 /* create new band of correct type */
3367 band = rt_band_new_inline(
3368 _width, _height,
3369 arg->pixtype[i],
3370 hasnodata, nodataval,
3371 data
3372 );
3373 if (band == NULL) {
3374 rterror("rt_raster_gdal_rasterize: Could not create band");
3375 rtdealloc(data);
3376 _rti_rasterize_arg_destroy(arg);
3377 rt_raster_destroy(rast);
3378 return NULL;
3379 }
3380
3381 /* give ownership of data to band */
3382 rt_band_set_ownsdata_flag(band, 1);
3383
3384 /* copy pixel by pixel */
3385 for (x = 0; x < _width; x++) {
3386 for (y = 0; y < _height; y++) {
3387 err = rt_band_get_pixel(oldband, x, y, &val, &nodata);
3388 if (err != ES_NONE) {
3389 rterror("rt_raster_gdal_rasterize: Could not get pixel value");
3390 _rti_rasterize_arg_destroy(arg);
3391 rt_raster_destroy(rast);
3392 rt_band_destroy(band);
3393 return NULL;
3394 }
3395
3396 if (nodata)
3397 val = nodataval;
3398
3399 err = rt_band_set_pixel(band, x, y, val, NULL);
3400 if (err != ES_NONE) {
3401 rterror("rt_raster_gdal_rasterize: Could not set pixel value");
3402 _rti_rasterize_arg_destroy(arg);
3403 rt_raster_destroy(rast);
3404 rt_band_destroy(band);
3405 return NULL;
3406 }
3407 }
3408 }
3409
3410 /* replace band */
3411 oldband = rt_raster_replace_band(rast, band, i);
3412 if (oldband == NULL) {
3413 rterror("rt_raster_gdal_rasterize: Could not replace band %d of output raster", i);
3414 _rti_rasterize_arg_destroy(arg);
3415 rt_raster_destroy(rast);
3416 rt_band_destroy(band);
3417 return NULL;
3418 }
3419
3420 /* free oldband */
3421 rt_band_destroy(oldband);
3422 }
3423
3424 _rti_rasterize_arg_destroy(arg);
3425
3426 RASTER_DEBUG(3, "done");
3427
3428 return rast;
3429 }
3430
3431 /******************************************************************************
3432 * rt_raster_from_two_rasters()
3433 ******************************************************************************/
3434
3435 /*
3436 * Return raster of computed extent specified extenttype applied
3437 * on two input rasters. The raster returned should be freed by
3438 * the caller
3439 *
3440 * @param rast1 : the first raster
3441 * @param rast2 : the second raster
3442 * @param extenttype : type of extent for the output raster
3443 * @param *rtnraster : raster of computed extent
3444 * @param *offset : 4-element array indicating the X,Y offsets
3445 * for each raster. 0,1 for rast1 X,Y. 2,3 for rast2 X,Y.
3446 *
3447 * @return ES_NONE if success, ES_ERROR if error
3448 */
3449 rt_errorstate
rt_raster_from_two_rasters(rt_raster rast1,rt_raster rast2,rt_extenttype extenttype,rt_raster * rtnraster,double * offset)3450 rt_raster_from_two_rasters(
3451 rt_raster rast1, rt_raster rast2,
3452 rt_extenttype extenttype,
3453 rt_raster *rtnraster, double *offset
3454 ) {
3455 int i;
3456
3457 rt_raster _rast[2] = {rast1, rast2};
3458 double _offset[2][4] = {{0.}};
3459 uint16_t _dim[2][2] = {{0}};
3460
3461 rt_raster raster = NULL;
3462 int aligned = 0;
3463 int dim[2] = {0};
3464 double gt[6] = {0.};
3465
3466 assert(NULL != rast1);
3467 assert(NULL != rast2);
3468 assert(NULL != rtnraster);
3469
3470 /* set rtnraster to NULL */
3471 *rtnraster = NULL;
3472
3473 /* rasters must be aligned */
3474 if (rt_raster_same_alignment(rast1, rast2, &aligned, NULL) != ES_NONE) {
3475 rterror("rt_raster_from_two_rasters: Could not test for alignment on the two rasters");
3476 return ES_ERROR;
3477 }
3478 if (!aligned) {
3479 rterror("rt_raster_from_two_rasters: The two rasters provided do not have the same alignment");
3480 return ES_ERROR;
3481 }
3482
3483 /* dimensions */
3484 _dim[0][0] = rast1->width;
3485 _dim[0][1] = rast1->height;
3486 _dim[1][0] = rast2->width;
3487 _dim[1][1] = rast2->height;
3488
3489 /* get raster offsets */
3490 if (rt_raster_geopoint_to_cell(
3491 _rast[1],
3492 _rast[0]->ipX, _rast[0]->ipY,
3493 &(_offset[1][0]), &(_offset[1][1]),
3494 NULL
3495 ) != ES_NONE) {
3496 rterror("rt_raster_from_two_rasters: Could not compute offsets of the second raster relative to the first raster");
3497 return ES_ERROR;
3498 }
3499 _offset[1][0] = -1 * _offset[1][0];
3500 _offset[1][1] = -1 * _offset[1][1];
3501 _offset[1][2] = _offset[1][0] + _dim[1][0] - 1;
3502 _offset[1][3] = _offset[1][1] + _dim[1][1] - 1;
3503
3504 i = -1;
3505 switch (extenttype) {
3506 case ET_FIRST:
3507 i = 0;
3508 _offset[0][0] = 0.;
3509 _offset[0][1] = 0.;
3510 /* FALLTHROUGH */
3511 case ET_LAST:
3512 case ET_SECOND:
3513 if (i < 0) {
3514 i = 1;
3515 _offset[0][0] = -1 * _offset[1][0];
3516 _offset[0][1] = -1 * _offset[1][1];
3517 _offset[1][0] = 0.;
3518 _offset[1][1] = 0.;
3519 }
3520
3521 dim[0] = _dim[i][0];
3522 dim[1] = _dim[i][1];
3523 raster = rt_raster_new(
3524 dim[0],
3525 dim[1]
3526 );
3527 if (raster == NULL) {
3528 rterror("rt_raster_from_two_rasters: Could not create output raster");
3529 return ES_ERROR;
3530 }
3531 rt_raster_set_srid(raster, _rast[i]->srid);
3532 rt_raster_get_geotransform_matrix(_rast[i], gt);
3533 rt_raster_set_geotransform_matrix(raster, gt);
3534 break;
3535 case ET_UNION: {
3536 double off[4] = {0};
3537
3538 rt_raster_get_geotransform_matrix(_rast[0], gt);
3539 RASTER_DEBUGF(4, "gt = (%f, %f, %f, %f, %f, %f)",
3540 gt[0],
3541 gt[1],
3542 gt[2],
3543 gt[3],
3544 gt[4],
3545 gt[5]
3546 );
3547
3548 /* new raster upper left offset */
3549 off[0] = 0;
3550 if (_offset[1][0] < 0)
3551 off[0] = _offset[1][0];
3552 off[1] = 0;
3553 if (_offset[1][1] < 0)
3554 off[1] = _offset[1][1];
3555
3556 /* new raster lower right offset */
3557 off[2] = _dim[0][0] - 1;
3558 if ((int) _offset[1][2] >= _dim[0][0])
3559 off[2] = _offset[1][2];
3560 off[3] = _dim[0][1] - 1;
3561 if ((int) _offset[1][3] >= _dim[0][1])
3562 off[3] = _offset[1][3];
3563
3564 /* upper left corner */
3565 if (rt_raster_cell_to_geopoint(
3566 _rast[0],
3567 off[0], off[1],
3568 &(gt[0]), &(gt[3]),
3569 NULL
3570 ) != ES_NONE) {
3571 rterror("rt_raster_from_two_rasters: Could not get spatial coordinates of upper-left pixel of output raster");
3572 return ES_ERROR;
3573 }
3574
3575 dim[0] = off[2] - off[0] + 1;
3576 dim[1] = off[3] - off[1] + 1;
3577 RASTER_DEBUGF(4, "off = (%f, %f, %f, %f)",
3578 off[0],
3579 off[1],
3580 off[2],
3581 off[3]
3582 );
3583 RASTER_DEBUGF(4, "dim = (%d, %d)", dim[0], dim[1]);
3584
3585 raster = rt_raster_new(
3586 dim[0],
3587 dim[1]
3588 );
3589 if (raster == NULL) {
3590 rterror("rt_raster_from_two_rasters: Could not create output raster");
3591 return ES_ERROR;
3592 }
3593 rt_raster_set_srid(raster, _rast[0]->srid);
3594 rt_raster_set_geotransform_matrix(raster, gt);
3595 RASTER_DEBUGF(4, "gt = (%f, %f, %f, %f, %f, %f)",
3596 gt[0],
3597 gt[1],
3598 gt[2],
3599 gt[3],
3600 gt[4],
3601 gt[5]
3602 );
3603
3604 /* get offsets */
3605 if (rt_raster_geopoint_to_cell(
3606 _rast[0],
3607 gt[0], gt[3],
3608 &(_offset[0][0]), &(_offset[0][1]),
3609 NULL
3610 ) != ES_NONE) {
3611 rterror("rt_raster_from_two_rasters: Could not get offsets of the FIRST raster relative to the output raster");
3612 rt_raster_destroy(raster);
3613 return ES_ERROR;
3614 }
3615 _offset[0][0] *= -1;
3616 _offset[0][1] *= -1;
3617
3618 if (rt_raster_geopoint_to_cell(
3619 _rast[1],
3620 gt[0], gt[3],
3621 &(_offset[1][0]), &(_offset[1][1]),
3622 NULL
3623 ) != ES_NONE) {
3624 rterror("rt_raster_from_two_rasters: Could not get offsets of the SECOND raster relative to the output raster");
3625 rt_raster_destroy(raster);
3626 return ES_ERROR;
3627 }
3628 _offset[1][0] *= -1;
3629 _offset[1][1] *= -1;
3630 break;
3631 }
3632 case ET_INTERSECTION: {
3633 double off[4] = {0};
3634
3635 /* no intersection */
3636 if (
3637 (_offset[1][2] < 0 || _offset[1][0] > (_dim[0][0] - 1)) ||
3638 (_offset[1][3] < 0 || _offset[1][1] > (_dim[0][1] - 1))
3639 ) {
3640 RASTER_DEBUG(3, "The two rasters provided have no intersection. Returning no band raster");
3641
3642 raster = rt_raster_new(0, 0);
3643 if (raster == NULL) {
3644 rterror("rt_raster_from_two_rasters: Could not create output raster");
3645 return ES_ERROR;
3646 }
3647 rt_raster_set_srid(raster, _rast[0]->srid);
3648 rt_raster_set_scale(raster, 0, 0);
3649
3650 /* set offsets if provided */
3651 if (NULL != offset) {
3652 for (i = 0; i < 4; i++)
3653 offset[i] = _offset[i / 2][i % 2];
3654 }
3655
3656 *rtnraster = raster;
3657 return ES_NONE;
3658 }
3659
3660 if (_offset[1][0] > 0)
3661 off[0] = _offset[1][0];
3662 if (_offset[1][1] > 0)
3663 off[1] = _offset[1][1];
3664
3665 off[2] = _dim[0][0] - 1;
3666 if (_offset[1][2] < _dim[0][0])
3667 off[2] = _offset[1][2];
3668 off[3] = _dim[0][1] - 1;
3669 if (_offset[1][3] < _dim[0][1])
3670 off[3] = _offset[1][3];
3671
3672 dim[0] = off[2] - off[0] + 1;
3673 dim[1] = off[3] - off[1] + 1;
3674 raster = rt_raster_new(
3675 dim[0],
3676 dim[1]
3677 );
3678 if (raster == NULL) {
3679 rterror("rt_raster_from_two_rasters: Could not create output raster");
3680 return ES_ERROR;
3681 }
3682 rt_raster_set_srid(raster, _rast[0]->srid);
3683
3684 /* get upper-left corner */
3685 rt_raster_get_geotransform_matrix(_rast[0], gt);
3686 if (rt_raster_cell_to_geopoint(
3687 _rast[0],
3688 off[0], off[1],
3689 &(gt[0]), &(gt[3]),
3690 gt
3691 ) != ES_NONE) {
3692 rterror("rt_raster_from_two_rasters: Could not get spatial coordinates of upper-left pixel of output raster");
3693 rt_raster_destroy(raster);
3694 return ES_ERROR;
3695 }
3696
3697 rt_raster_set_geotransform_matrix(raster, gt);
3698
3699 /* get offsets */
3700 if (rt_raster_geopoint_to_cell(
3701 _rast[0],
3702 gt[0], gt[3],
3703 &(_offset[0][0]), &(_offset[0][1]),
3704 NULL
3705 ) != ES_NONE) {
3706 rterror("rt_raster_from_two_rasters: Could not get pixel coordinates to compute the offsets of the FIRST raster relative to the output raster");
3707 rt_raster_destroy(raster);
3708 return ES_ERROR;
3709 }
3710 _offset[0][0] *= -1;
3711 _offset[0][1] *= -1;
3712
3713 if (rt_raster_geopoint_to_cell(
3714 _rast[1],
3715 gt[0], gt[3],
3716 &(_offset[1][0]), &(_offset[1][1]),
3717 NULL
3718 ) != ES_NONE) {
3719 rterror("rt_raster_from_two_rasters: Could not get pixel coordinates to compute the offsets of the SECOND raster relative to the output raster");
3720 rt_raster_destroy(raster);
3721 return ES_ERROR;
3722 }
3723 _offset[1][0] *= -1;
3724 _offset[1][1] *= -1;
3725 break;
3726 }
3727 case ET_CUSTOM:
3728 rterror("rt_raster_from_two_rasters: Extent type ET_CUSTOM is not supported by this function");
3729 break;
3730 }
3731
3732 /* set offsets if provided */
3733 if (NULL != offset) {
3734 for (i = 0; i < 4; i++)
3735 offset[i] = _offset[i / 2][i % 2];
3736 }
3737
3738 *rtnraster = raster;
3739 return ES_NONE;
3740 }
3741