PostGIS  2.4.9dev-r@@SVN_REVISION@@
rt_raster.c
Go to the documentation of this file.
1 /*
2  *
3  * WKTRaster - Raster Types for PostGIS
4  * http://trac.osgeo.org/postgis/wiki/WKTRaster
5  *
6  * Copyright (C) 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 
47 rt_raster
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
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
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
121 rt_raster_get_width(rt_raster raster) {
122 
123  assert(NULL != raster);
124 
125  return raster->width;
126 }
127 
128 uint16_t
129 rt_raster_get_height(rt_raster raster) {
130 
131  assert(NULL != raster);
132 
133  return raster->height;
134 }
135 
136 void
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
150 rt_raster_get_x_scale(rt_raster raster) {
151 
152 
153  assert(NULL != raster);
154 
155  return raster->scaleX;
156 }
157 
158 double
159 rt_raster_get_y_scale(rt_raster raster) {
160 
161 
162  assert(NULL != raster);
163 
164  return raster->scaleY;
165 }
166 
167 void
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
181 rt_raster_get_x_skew(rt_raster raster) {
182 
183 
184  assert(NULL != raster);
185 
186  return raster->skewX;
187 }
188 
189 double
190 rt_raster_get_y_skew(rt_raster raster) {
191 
192 
193  assert(NULL != raster);
194 
195  return raster->skewY;
196 }
197 
198 void
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
213 rt_raster_get_x_offset(rt_raster raster) {
214 
215 
216  assert(NULL != raster);
217 
218  return raster->ipX;
219 }
220 
221 double
222 rt_raster_get_y_offset(rt_raster raster) {
223 
224 
225  assert(NULL != raster);
226 
227  return raster->ipY;
228 }
229 
230 void
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
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
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
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
356 rt_raster_get_srid(rt_raster raster) {
357  assert(NULL != raster);
358 
359  return clamp_srid(raster->srid);
360 }
361 
362 void
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 int
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
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 
404 int
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 
484 int
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 
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 
705 void
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 
726 void
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 
754 rt_errorstate
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 (
771  FLT_EQ(_gt[1], 0) ||
772  FLT_EQ(_gt[5], 0)
773  ) {
774  rt_raster_get_geotransform_matrix(raster, _gt);
775  }
776 
777  RASTER_DEBUGF(4, "gt = (%f, %f, %f, %f, %f, %f)",
778  _gt[0],
779  _gt[1],
780  _gt[2],
781  _gt[3],
782  _gt[4],
783  _gt[5]
784  );
785 
786  GDALApplyGeoTransform(_gt, xr, yr, xw, yw);
787  RASTER_DEBUGF(4, "GDALApplyGeoTransform (c -> g) for (%f, %f) = (%f, %f)",
788  xr, yr, *xw, *yw);
789 
790  return ES_NONE;
791 }
792 
805 rt_errorstate
806 rt_raster_geopoint_to_cell(
807  rt_raster raster,
808  double xw, double yw,
809  double *xr, double *yr,
810  double *igt
811 ) {
812  double _igt[6] = {0};
813  double rnd = 0;
814 
815  assert(NULL != raster);
816  assert(NULL != xr && NULL != yr);
817 
818  if (igt != NULL)
819  memcpy(_igt, igt, sizeof(double) * 6);
820 
821  /* matrix is not set */
822  if (
823  FLT_EQ(_igt[0], 0.) &&
824  FLT_EQ(_igt[1], 0.) &&
825  FLT_EQ(_igt[2], 0.) &&
826  FLT_EQ(_igt[3], 0.) &&
827  FLT_EQ(_igt[4], 0.) &&
828  FLT_EQ(_igt[5], 0.)
829  ) {
830  if (rt_raster_get_inverse_geotransform_matrix(raster, NULL, _igt) != ES_NONE) {
831  rterror("rt_raster_geopoint_to_cell: Could not get inverse geotransform matrix");
832  return ES_ERROR;
833  }
834  }
835 
836  GDALApplyGeoTransform(_igt, xw, yw, xr, yr);
837  RASTER_DEBUGF(4, "GDALApplyGeoTransform (g -> c) for (%f, %f) = (%f, %f)",
838  xw, yw, *xr, *yr);
839 
840  rnd = ROUND(*xr, 0);
841  if (FLT_EQ(rnd, *xr))
842  *xr = rnd;
843  else
844  *xr = floor(*xr);
845 
846  rnd = ROUND(*yr, 0);
847  if (FLT_EQ(rnd, *yr))
848  *yr = rnd;
849  else
850  *yr = floor(*yr);
851 
852  RASTER_DEBUGF(4, "Corrected GDALApplyGeoTransform (g -> c) for (%f, %f) = (%f, %f)",
853  xw, yw, *xr, *yr);
854 
855  return ES_NONE;
856 }
857 
858 /******************************************************************************
859 * rt_raster_get_envelope()
860 ******************************************************************************/
861 
872 rt_errorstate
873 rt_raster_get_envelope(
874  rt_raster raster,
875  rt_envelope *env
876 ) {
877  int i;
878  int rtn;
879  int set = 0;
880  double _r[2] = {0.};
881  double _w[2] = {0.};
882  double _gt[6] = {0.};
883 
884  assert(raster != NULL);
885  assert(env != NULL);
886 
887  rt_raster_get_geotransform_matrix(raster, _gt);
888 
889  for (i = 0; i < 4; i++) {
890  switch (i) {
891  case 0:
892  _r[0] = 0;
893  _r[1] = 0;
894  break;
895  case 1:
896  _r[0] = 0;
897  _r[1] = raster->height;
898  break;
899  case 2:
900  _r[0] = raster->width;
901  _r[1] = raster->height;
902  break;
903  case 3:
904  _r[0] = raster->width;
905  _r[1] = 0;
906  break;
907  }
908 
909  rtn = rt_raster_cell_to_geopoint(
910  raster,
911  _r[0], _r[1],
912  &(_w[0]), &(_w[1]),
913  _gt
914  );
915  if (rtn != ES_NONE) {
916  rterror("rt_raster_get_envelope: Could not compute spatial coordinates for raster pixel");
917  return ES_ERROR;
918  }
919 
920  if (!set) {
921  set = 1;
922  env->MinX = _w[0];
923  env->MaxX = _w[0];
924  env->MinY = _w[1];
925  env->MaxY = _w[1];
926  }
927  else {
928  if (_w[0] < env->MinX)
929  env->MinX = _w[0];
930  else if (_w[0] > env->MaxX)
931  env->MaxX = _w[0];
932 
933  if (_w[1] < env->MinY)
934  env->MinY = _w[1];
935  else if (_w[1] > env->MaxY)
936  env->MaxY = _w[1];
937  }
938  }
939 
940  return ES_NONE;
941 }
942 
943 /******************************************************************************
944 * rt_raster_compute_skewed_raster()
945 ******************************************************************************/
946 
947 /*
948  * Compute skewed extent that covers unskewed extent.
949  *
950  * @param envelope : unskewed extent of type rt_envelope
951  * @param skew : pointer to 2-element array (x, y) of skew
952  * @param scale : pointer to 2-element array (x, y) of scale
953  * @param tolerance : value between 0 and 1 where the smaller the tolerance
954  * results in an extent approaching the "minimum" skewed extent.
955  * If value <= 0, tolerance = 0.1. If value > 1, tolerance = 1.
956  *
957  * @return skewed raster who's extent covers unskewed extent, NULL on error
958  */
959 rt_raster
960 rt_raster_compute_skewed_raster(
961  rt_envelope extent,
962  double *skew,
963  double *scale,
964  double tolerance
965 ) {
966  uint32_t run = 0;
967  uint32_t max_run = 1;
968  double dbl_run = 0;
969 
970  int rtn;
971  int covers = 0;
972  rt_raster raster;
973  double _gt[6] = {0};
974  double _igt[6] = {0};
975  int _d[2] = {1, -1};
976  int _dlast = 0;
977  int _dlastpos = 0;
978  double _w[2] = {0};
979  double _r[2] = {0};
980  double _xy[2] = {0};
981  int i;
982  int j;
983  int x;
984  int y;
985 
986  LWGEOM *geom = NULL;
987  GEOSGeometry *sgeom = NULL;
988  GEOSGeometry *ngeom = NULL;
989 
990  if (
991  (tolerance < 0.) ||
992  FLT_EQ(tolerance, 0.)
993  ) {
994  tolerance = 0.1;
995  }
996  else if (tolerance > 1.)
997  tolerance = 1;
998 
999  dbl_run = tolerance;
1000  while (dbl_run < 10) {
1001  dbl_run *= 10.;
1002  max_run *= 10;
1003  }
1004 
1005  /* scale must be provided */
1006  if (scale == NULL)
1007  return NULL;
1008  for (i = 0; i < 2; i++) {
1009  if (FLT_EQ(scale[i], 0)) {
1010  rterror("rt_raster_compute_skewed_raster: Scale cannot be zero");
1011  return 0;
1012  }
1013 
1014  if (i < 1)
1015  _gt[1] = fabs(scale[i] * tolerance);
1016  else
1017  _gt[5] = fabs(scale[i] * tolerance);
1018  }
1019  /* conform scale-y to be negative */
1020  _gt[5] *= -1;
1021 
1022  /* skew not provided or skew is zero, return raster of correct dim and spatial attributes */
1023  if (
1024  (skew == NULL) || (
1025  FLT_EQ(skew[0], 0) &&
1026  FLT_EQ(skew[1], 0)
1027  )
1028  ) {
1029  int _dim[2] = {
1030  (int) fmax((fabs(extent.MaxX - extent.MinX) + (fabs(scale[0]) / 2.)) / fabs(scale[0]), 1),
1031  (int) fmax((fabs(extent.MaxY - extent.MinY) + (fabs(scale[1]) / 2.)) / fabs(scale[1]), 1)
1032  };
1033 
1034  raster = rt_raster_new(_dim[0], _dim[1]);
1035  if (raster == NULL) {
1036  rterror("rt_raster_compute_skewed_raster: Could not create output raster");
1037  return NULL;
1038  }
1039 
1040  rt_raster_set_offsets(raster, extent.MinX, extent.MaxY);
1041  rt_raster_set_scale(raster, fabs(scale[0]), -1 * fabs(scale[1]));
1042  rt_raster_set_skews(raster, skew[0], skew[1]);
1043 
1044  return raster;
1045  }
1046 
1047  /* direction to shift upper-left corner */
1048  if (skew[0] > 0.)
1049  _d[0] = -1;
1050  if (skew[1] < 0.)
1051  _d[1] = 1;
1052 
1053  /* geotransform */
1054  _gt[0] = extent.UpperLeftX;
1055  _gt[2] = skew[0] * tolerance;
1056  _gt[3] = extent.UpperLeftY;
1057  _gt[4] = skew[1] * tolerance;
1058 
1059  RASTER_DEBUGF(4, "Initial geotransform: %f, %f, %f, %f, %f, %f",
1060  _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]
1061  );
1062  RASTER_DEBUGF(4, "Delta: %d, %d", _d[0], _d[1]);
1063 
1064  /* simple raster */
1065  if ((raster = rt_raster_new(1, 1)) == NULL) {
1066  rterror("rt_raster_compute_skewed_raster: Out of memory allocating extent raster");
1067  return NULL;
1068  }
1069  rt_raster_set_geotransform_matrix(raster, _gt);
1070 
1071  /* get inverse geotransform matrix */
1072  if (!GDALInvGeoTransform(_gt, _igt)) {
1073  rterror("rt_raster_compute_skewed_raster: Could not compute inverse geotransform matrix");
1074  rt_raster_destroy(raster);
1075  return NULL;
1076  }
1077  RASTER_DEBUGF(4, "Inverse geotransform: %f, %f, %f, %f, %f, %f",
1078  _igt[0], _igt[1], _igt[2], _igt[3], _igt[4], _igt[5]
1079  );
1080 
1081  /* shift along axis */
1082  for (i = 0; i < 2; i++) {
1083  covers = 0;
1084  run = 0;
1085 
1086  RASTER_DEBUGF(3, "Shifting along %s axis", i < 1 ? "X" : "Y");
1087 
1088  do {
1089 
1090  /* prevent possible infinite loop */
1091  if (run > max_run) {
1092  rterror("rt_raster_compute_skewed_raster: Could not compute skewed extent due to check preventing infinite loop");
1093  rt_raster_destroy(raster);
1094  return NULL;
1095  }
1096 
1097  /*
1098  check the four corners that they are covered along the specific axis
1099  pixel column should be >= 0
1100  */
1101  for (j = 0; j < 4; j++) {
1102  switch (j) {
1103  /* upper-left */
1104  case 0:
1105  _xy[0] = extent.MinX;
1106  _xy[1] = extent.MaxY;
1107  break;
1108  /* lower-left */
1109  case 1:
1110  _xy[0] = extent.MinX;
1111  _xy[1] = extent.MinY;
1112  break;
1113  /* lower-right */
1114  case 2:
1115  _xy[0] = extent.MaxX;
1116  _xy[1] = extent.MinY;
1117  break;
1118  /* upper-right */
1119  case 3:
1120  _xy[0] = extent.MaxX;
1121  _xy[1] = extent.MaxY;
1122  break;
1123  }
1124 
1125  rtn = rt_raster_geopoint_to_cell(
1126  raster,
1127  _xy[0], _xy[1],
1128  &(_r[0]), &(_r[1]),
1129  _igt
1130  );
1131  if (rtn != ES_NONE) {
1132  rterror("rt_raster_compute_skewed_raster: Could not compute raster pixel for spatial coordinates");
1133  rt_raster_destroy(raster);
1134  return NULL;
1135  }
1136 
1137  RASTER_DEBUGF(4, "Point %d at cell %d x %d", j, (int) _r[0], (int) _r[1]);
1138 
1139  /* raster doesn't cover point */
1140  if ((int) _r[i] < 0) {
1141  RASTER_DEBUGF(4, "Point outside of skewed extent: %d", j);
1142  covers = 0;
1143 
1144  if (_dlastpos != j) {
1145  _dlast = (int) _r[i];
1146  _dlastpos = j;
1147  }
1148  else if ((int) _r[i] < _dlast) {
1149  RASTER_DEBUG(4, "Point going in wrong direction. Reversing direction");
1150  _d[i] *= -1;
1151  _dlastpos = -1;
1152  run = 0;
1153  }
1154 
1155  break;
1156  }
1157 
1158  covers++;
1159  }
1160 
1161  if (!covers) {
1162  x = 0;
1163  y = 0;
1164  if (i < 1)
1165  x = _d[i] * fabs(_r[i]);
1166  else
1167  y = _d[i] * fabs(_r[i]);
1168 
1169  rtn = rt_raster_cell_to_geopoint(
1170  raster,
1171  x, y,
1172  &(_w[0]), &(_w[1]),
1173  _gt
1174  );
1175  if (rtn != ES_NONE) {
1176  rterror("rt_raster_compute_skewed_raster: Could not compute spatial coordinates for raster pixel");
1177  rt_raster_destroy(raster);
1178  return NULL;
1179  }
1180 
1181  /* adjust ul */
1182  if (i < 1)
1183  _gt[0] = _w[i];
1184  else
1185  _gt[3] = _w[i];
1186  rt_raster_set_geotransform_matrix(raster, _gt);
1187  RASTER_DEBUGF(4, "Shifted geotransform: %f, %f, %f, %f, %f, %f",
1188  _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]
1189  );
1190 
1191  /* get inverse geotransform matrix */
1192  if (!GDALInvGeoTransform(_gt, _igt)) {
1193  rterror("rt_raster_compute_skewed_raster: Could not compute inverse geotransform matrix");
1194  rt_raster_destroy(raster);
1195  return NULL;
1196  }
1197  RASTER_DEBUGF(4, "Inverse geotransform: %f, %f, %f, %f, %f, %f",
1198  _igt[0], _igt[1], _igt[2], _igt[3], _igt[4], _igt[5]
1199  );
1200  }
1201 
1202  run++;
1203  }
1204  while (!covers);
1205  }
1206 
1207  /* covers test */
1208  rtn = rt_raster_geopoint_to_cell(
1209  raster,
1210  extent.MaxX, extent.MinY,
1211  &(_r[0]), &(_r[1]),
1212  _igt
1213  );
1214  if (rtn != ES_NONE) {
1215  rterror("rt_raster_compute_skewed_raster: Could not compute raster pixel for spatial coordinates");
1216  rt_raster_destroy(raster);
1217  return NULL;
1218  }
1219 
1220  RASTER_DEBUGF(4, "geopoint %f x %f at cell %d x %d", extent.MaxX, extent.MinY, (int) _r[0], (int) _r[1]);
1221 
1222  raster->width = _r[0];
1223  raster->height = _r[1];
1224 
1225  /* initialize GEOS */
1226  initGEOS(rtinfo, lwgeom_geos_error);
1227 
1228  /* create reference LWPOLY */
1229  {
1230  LWPOLY *npoly = rt_util_envelope_to_lwpoly(extent);
1231  if (npoly == NULL) {
1232  rterror("rt_raster_compute_skewed_raster: Could not build extent's geometry for covers test");
1233  rt_raster_destroy(raster);
1234  return NULL;
1235  }
1236 
1237  ngeom = (GEOSGeometry *) LWGEOM2GEOS(lwpoly_as_lwgeom(npoly), 0);
1238  lwpoly_free(npoly);
1239  }
1240 
1241  do {
1242  covers = 0;
1243 
1244  /* construct sgeom from raster */
1245  if ((rt_raster_get_convex_hull(raster, &geom) != ES_NONE) || geom == NULL) {
1246  rterror("rt_raster_compute_skewed_raster: Could not build skewed extent's geometry for covers test");
1247  GEOSGeom_destroy(ngeom);
1248  rt_raster_destroy(raster);
1249  return NULL;
1250  }
1251 
1252  sgeom = (GEOSGeometry *) LWGEOM2GEOS(geom, 0);
1253  lwgeom_free(geom);
1254 
1255  covers = GEOSRelatePattern(sgeom, ngeom, "******FF*");
1256  GEOSGeom_destroy(sgeom);
1257 
1258  if (covers == 2) {
1259  rterror("rt_raster_compute_skewed_raster: Could not run covers test");
1260  GEOSGeom_destroy(ngeom);
1261  rt_raster_destroy(raster);
1262  return NULL;
1263  }
1264 
1265  if (covers)
1266  break;
1267 
1268  raster->width++;
1269  raster->height++;
1270  }
1271  while (!covers);
1272 
1273  RASTER_DEBUGF(4, "Skewed extent does cover normal extent with dimensions %d x %d", raster->width, raster->height);
1274 
1275  raster->width = (int) ((((double) raster->width) * fabs(_gt[1]) + fabs(scale[0] / 2.)) / fabs(scale[0]));
1276  raster->height = (int) ((((double) raster->height) * fabs(_gt[5]) + fabs(scale[1] / 2.)) / fabs(scale[1]));
1277  _gt[1] = fabs(scale[0]);
1278  _gt[5] = -1 * fabs(scale[1]);
1279  _gt[2] = skew[0];
1280  _gt[4] = skew[1];
1281  rt_raster_set_geotransform_matrix(raster, _gt);
1282 
1283  /* minimize width/height */
1284  for (i = 0; i < 2; i++) {
1285  covers = 1;
1286  do {
1287  if (i < 1)
1288  raster->width--;
1289  else
1290  raster->height--;
1291 
1292  /* construct sgeom from raster */
1293  if ((rt_raster_get_convex_hull(raster, &geom) != ES_NONE) || geom == NULL) {
1294  rterror("rt_raster_compute_skewed_raster: Could not build skewed extent's geometry for minimizing dimensions");
1295  GEOSGeom_destroy(ngeom);
1296  rt_raster_destroy(raster);
1297  return NULL;
1298  }
1299 
1300  sgeom = (GEOSGeometry *) LWGEOM2GEOS(geom, 0);
1301  lwgeom_free(geom);
1302 
1303  covers = GEOSRelatePattern(sgeom, ngeom, "******FF*");
1304  GEOSGeom_destroy(sgeom);
1305 
1306  if (covers == 2) {
1307  rterror("rt_raster_compute_skewed_raster: Could not run covers test for minimizing dimensions");
1308  GEOSGeom_destroy(ngeom);
1309  rt_raster_destroy(raster);
1310  return NULL;
1311  }
1312 
1313  if (!covers) {
1314  if (i < 1)
1315  raster->width++;
1316  else
1317  raster->height++;
1318 
1319  break;
1320  }
1321  }
1322  while (covers);
1323  }
1324 
1325  GEOSGeom_destroy(ngeom);
1326 
1327  return raster;
1328 }
1329 
1337 int
1338 rt_raster_is_empty(rt_raster raster) {
1339  return (NULL == raster || raster->height <= 0 || raster->width <= 0);
1340 }
1341 
1350 int
1351 rt_raster_has_band(rt_raster raster, int nband) {
1352  return !(NULL == raster || nband >= raster->numBands || nband < 0);
1353 }
1354 
1355 /******************************************************************************
1356 * rt_raster_copy_band()
1357 ******************************************************************************/
1358 
1373 int
1374 rt_raster_copy_band(
1375  rt_raster torast, rt_raster fromrast,
1376  int fromindex, int toindex
1377 ) {
1378  rt_band srcband = NULL;
1379  rt_band dstband = NULL;
1380 
1381  assert(NULL != torast);
1382  assert(NULL != fromrast);
1383 
1384  /* Check raster dimensions */
1385  if (torast->height != fromrast->height || torast->width != fromrast->width) {
1386  rtwarn("rt_raster_copy_band: Attempting to add a band with different width or height");
1387  return -1;
1388  }
1389 
1390  /* Check bands limits */
1391  if (fromrast->numBands < 1) {
1392  rtwarn("rt_raster_copy_band: Second raster has no band");
1393  return -1;
1394  }
1395  else if (fromindex < 0) {
1396  rtwarn("rt_raster_copy_band: Band index for second raster < 0. Defaulted to 0");
1397  fromindex = 0;
1398  }
1399  else if (fromindex >= fromrast->numBands) {
1400  rtwarn("rt_raster_copy_band: Band index for second raster > number of bands, truncated from %u to %u", fromindex, fromrast->numBands - 1);
1401  fromindex = fromrast->numBands - 1;
1402  }
1403 
1404  if (toindex < 0) {
1405  rtwarn("rt_raster_copy_band: Band index for first raster < 0. Defaulted to 0");
1406  toindex = 0;
1407  }
1408  else if (toindex > torast->numBands) {
1409  rtwarn("rt_raster_copy_band: Band index for first raster > number of bands, truncated from %u to %u", toindex, torast->numBands);
1410  toindex = torast->numBands;
1411  }
1412 
1413  /* Get band from source raster */
1414  srcband = rt_raster_get_band(fromrast, fromindex);
1415 
1416  /* duplicate band */
1417  dstband = rt_band_duplicate(srcband);
1418 
1419  /* Add band to the second raster */
1420  return rt_raster_add_band(torast, dstband, toindex);
1421 }
1422 
1423 /******************************************************************************
1424 * rt_raster_from_band()
1425 ******************************************************************************/
1426 
1438 rt_raster
1439 rt_raster_from_band(rt_raster raster, uint32_t *bandNums, int count) {
1440  rt_raster rast = NULL;
1441  int i = 0;
1442  int j = 0;
1443  int idx;
1444  int32_t flag;
1445  double gt[6] = {0.};
1446 
1447  assert(NULL != raster);
1448  assert(NULL != bandNums);
1449 
1450  RASTER_DEBUGF(3, "rt_raster_from_band: source raster has %d bands",
1451  rt_raster_get_num_bands(raster));
1452 
1453  /* create new raster */
1454  rast = rt_raster_new(raster->width, raster->height);
1455  if (NULL == rast) {
1456  rterror("rt_raster_from_band: Out of memory allocating new raster");
1457  return NULL;
1458  }
1459 
1460  /* copy raster attributes */
1461  rt_raster_get_geotransform_matrix(raster, gt);
1462  rt_raster_set_geotransform_matrix(rast, gt);
1463 
1464  /* srid */
1465  rt_raster_set_srid(rast, raster->srid);
1466 
1467  /* copy bands */
1468  for (i = 0; i < count; i++) {
1469  idx = bandNums[i];
1470  flag = rt_raster_copy_band(rast, raster, idx, i);
1471 
1472  if (flag < 0) {
1473  rterror("rt_raster_from_band: Could not copy band");
1474  for (j = 0; j < i; j++) rt_band_destroy(rast->bands[j]);
1475  rt_raster_destroy(rast);
1476  return NULL;
1477  }
1478 
1479  RASTER_DEBUGF(3, "rt_raster_from_band: band created at index %d",
1480  flag);
1481  }
1482 
1483  RASTER_DEBUGF(3, "rt_raster_from_band: new raster has %d bands",
1484  rt_raster_get_num_bands(rast));
1485  return rast;
1486 }
1487 
1488 /******************************************************************************
1489 * rt_raster_replace_band()
1490 ******************************************************************************/
1491 
1501 rt_band
1502 rt_raster_replace_band(rt_raster raster, rt_band band, int index) {
1503  rt_band oldband = NULL;
1504  assert(NULL != raster);
1505  assert(NULL != band);
1506 
1507  if (band->width != raster->width || band->height != raster->height) {
1508  rterror("rt_raster_replace_band: Band does not match raster's dimensions: %dx%d band to %dx%d raster",
1509  band->width, band->height, raster->width, raster->height);
1510  return 0;
1511  }
1512 
1513  if (index >= raster->numBands || index < 0) {
1514  rterror("rt_raster_replace_band: Band index is not valid");
1515  return 0;
1516  }
1517 
1518  oldband = rt_raster_get_band(raster, index);
1519  RASTER_DEBUGF(3, "rt_raster_replace_band: old band at %p", oldband);
1520  RASTER_DEBUGF(3, "rt_raster_replace_band: new band at %p", band);
1521 
1522  raster->bands[index] = band;
1523  RASTER_DEBUGF(3, "rt_raster_replace_band: new band at %p", raster->bands[index]);
1524 
1525  band->raster = raster;
1526  oldband->raster = NULL;
1527 
1528  return oldband;
1529 }
1530 
1531 /******************************************************************************
1532 * rt_raster_clone()
1533 ******************************************************************************/
1534 
1543 rt_raster
1544 rt_raster_clone(rt_raster raster, uint8_t deep) {
1545  rt_raster rtn = NULL;
1546  double gt[6] = {0};
1547 
1548  assert(NULL != raster);
1549 
1550  if (deep) {
1551  int numband = rt_raster_get_num_bands(raster);
1552  uint32_t *nband = NULL;
1553  int i = 0;
1554 
1555  nband = rtalloc(sizeof(uint32_t) * numband);
1556  if (nband == NULL) {
1557  rterror("rt_raster_clone: Could not allocate memory for deep clone");
1558  return NULL;
1559  }
1560  for (i = 0; i < numband; i++)
1561  nband[i] = i;
1562 
1563  rtn = rt_raster_from_band(raster, nband, numband);
1564  rtdealloc(nband);
1565 
1566  return rtn;
1567  }
1568 
1569  rtn = rt_raster_new(
1570  rt_raster_get_width(raster),
1571  rt_raster_get_height(raster)
1572  );
1573  if (rtn == NULL) {
1574  rterror("rt_raster_clone: Could not create cloned raster");
1575  return NULL;
1576  }
1577 
1578  rt_raster_get_geotransform_matrix(raster, gt);
1579  rt_raster_set_geotransform_matrix(rtn, gt);
1580  rt_raster_set_srid(rtn, rt_raster_get_srid(raster));
1581 
1582  return rtn;
1583 }
1584 
1585 /******************************************************************************
1586 * rt_raster_to_gdal()
1587 ******************************************************************************/
1588 
1601 uint8_t*
1602 rt_raster_to_gdal(
1603  rt_raster raster, const char *srs,
1604  char *format, char **options, uint64_t *gdalsize
1605 ) {
1606  GDALDriverH src_drv = NULL;
1607  int destroy_src_drv = 0;
1608  GDALDatasetH src_ds = NULL;
1609 
1610  vsi_l_offset rtn_lenvsi;
1611  uint64_t rtn_len = 0;
1612 
1613  GDALDriverH rtn_drv = NULL;
1614  GDALDatasetH rtn_ds = NULL;
1615  uint8_t *rtn = NULL;
1616 
1617  assert(NULL != raster);
1618  assert(NULL != gdalsize);
1619 
1620  /* any supported format is possible */
1621  rt_util_gdal_register_all(0);
1622  RASTER_DEBUG(3, "loaded all supported GDAL formats");
1623 
1624  /* output format not specified */
1625  if (format == NULL || !strlen(format))
1626  format = "GTiff";
1627  RASTER_DEBUGF(3, "output format is %s", format);
1628 
1629  /* load raster into a GDAL MEM raster */
1630  src_ds = rt_raster_to_gdal_mem(raster, srs, NULL, NULL, 0, &src_drv, &destroy_src_drv);
1631  if (NULL == src_ds) {
1632  rterror("rt_raster_to_gdal: Could not convert raster to GDAL MEM format");
1633  return 0;
1634  }
1635 
1636  /* load driver */
1637  rtn_drv = GDALGetDriverByName(format);
1638  if (NULL == rtn_drv) {
1639  rterror("rt_raster_to_gdal: Could not load the output GDAL driver");
1640  GDALClose(src_ds);
1641  if (destroy_src_drv) GDALDestroyDriver(src_drv);
1642  return 0;
1643  }
1644  RASTER_DEBUG(3, "Output driver loaded");
1645 
1646  /* convert GDAL MEM raster to output format */
1647  RASTER_DEBUG(3, "Copying GDAL MEM raster to memory file in output format");
1648  rtn_ds = GDALCreateCopy(
1649  rtn_drv,
1650  "/vsimem/out.dat", /* should be fine assuming this is in a process */
1651  src_ds,
1652  FALSE, /* should copy be strictly equivelent? */
1653  options, /* format options */
1654  NULL, /* progress function */
1655  NULL /* progress data */
1656  );
1657 
1658  /* close source dataset */
1659  GDALClose(src_ds);
1660  if (destroy_src_drv) GDALDestroyDriver(src_drv);
1661  RASTER_DEBUG(3, "Closed GDAL MEM raster");
1662 
1663  if (NULL == rtn_ds) {
1664  rterror("rt_raster_to_gdal: Could not create the output GDAL dataset");
1665  return 0;
1666  }
1667 
1668  RASTER_DEBUGF(4, "dataset SRS: %s", GDALGetProjectionRef(rtn_ds));
1669 
1670  /* close dataset, this also flushes any pending writes */
1671  GDALClose(rtn_ds);
1672  RASTER_DEBUG(3, "Closed GDAL output raster");
1673 
1674  RASTER_DEBUG(3, "Done copying GDAL MEM raster to memory file in output format");
1675 
1676  /* from memory file to buffer */
1677  RASTER_DEBUG(3, "Copying GDAL memory file to buffer");
1678  rtn = VSIGetMemFileBuffer("/vsimem/out.dat", &rtn_lenvsi, TRUE);
1679  RASTER_DEBUG(3, "Done copying GDAL memory file to buffer");
1680  if (NULL == rtn) {
1681  rterror("rt_raster_to_gdal: Could not create the output GDAL raster");
1682  return 0;
1683  }
1684 
1685  rtn_len = (uint64_t) rtn_lenvsi;
1686  *gdalsize = rtn_len;
1687 
1688  return rtn;
1689 }
1690 
1691 /******************************************************************************
1692 * rt_raster_gdal_drivers()
1693 ******************************************************************************/
1694 
1704 rt_gdaldriver
1705 rt_raster_gdal_drivers(uint32_t *drv_count, uint8_t cancc) {
1706  const char *state;
1707  const char *txt;
1708  int txt_len;
1709  GDALDriverH *drv = NULL;
1710  rt_gdaldriver rtn = NULL;
1711  int count;
1712  int i;
1713  uint32_t j;
1714 
1715  assert(drv_count != NULL);
1716 
1717  rt_util_gdal_register_all(0);
1718  count = GDALGetDriverCount();
1719  RASTER_DEBUGF(3, "%d drivers found", count);
1720 
1721  rtn = (rt_gdaldriver) rtalloc(count * sizeof(struct rt_gdaldriver_t));
1722  if (NULL == rtn) {
1723  rterror("rt_raster_gdal_drivers: Could not allocate memory for gdaldriver structure");
1724  return 0;
1725  }
1726 
1727  for (i = 0, j = 0; i < count; i++) {
1728  drv = GDALGetDriver(i);
1729 
1730 #ifdef GDAL_DCAP_RASTER
1731  /* Starting with GDAL 2.0, vector drivers can also be returned */
1732  /* Only keep raster drivers */
1733  state = GDALGetMetadataItem(drv, GDAL_DCAP_RASTER, NULL);
1734  if (state == NULL || !EQUAL(state, "YES"))
1735  continue;
1736 #endif
1737 
1738  if (cancc) {
1739  /* CreateCopy support */
1740  state = GDALGetMetadataItem(drv, GDAL_DCAP_CREATECOPY, NULL);
1741  if (state == NULL) continue;
1742 
1743  /* VirtualIO support */
1744  state = GDALGetMetadataItem(drv, GDAL_DCAP_VIRTUALIO, NULL);
1745  if (state == NULL) continue;
1746  }
1747 
1748  /* index of driver */
1749  rtn[j].idx = i;
1750 
1751  /* short name */
1752  txt = GDALGetDriverShortName(drv);
1753  txt_len = strlen(txt);
1754 
1755  if (cancc) {
1756  RASTER_DEBUGF(3, "driver %s (%d) supports CreateCopy() and VirtualIO()", txt, i);
1757  }
1758 
1759  txt_len = (txt_len + 1) * sizeof(char);
1760  rtn[j].short_name = (char *) rtalloc(txt_len);
1761  memcpy(rtn[j].short_name, txt, txt_len);
1762 
1763  /* long name */
1764  txt = GDALGetDriverLongName(drv);
1765  txt_len = strlen(txt);
1766 
1767  txt_len = (txt_len + 1) * sizeof(char);
1768  rtn[j].long_name = (char *) rtalloc(txt_len);
1769  memcpy(rtn[j].long_name, txt, txt_len);
1770 
1771  /* creation options */
1772  txt = GDALGetDriverCreationOptionList(drv);
1773  txt_len = strlen(txt);
1774 
1775  txt_len = (txt_len + 1) * sizeof(char);
1776  rtn[j].create_options = (char *) rtalloc(txt_len);
1777  memcpy(rtn[j].create_options, txt, txt_len);
1778 
1779  j++;
1780  }
1781 
1782  /* free unused memory */
1783  rtn = rtrealloc(rtn, j * sizeof(struct rt_gdaldriver_t));
1784  *drv_count = j;
1785 
1786  return rtn;
1787 }
1788 
1789 /******************************************************************************
1790 * rt_raster_to_gdal_mem()
1791 ******************************************************************************/
1792 
1808 GDALDatasetH
1809 rt_raster_to_gdal_mem(
1810  rt_raster raster,
1811  const char *srs,
1812  uint32_t *bandNums,
1813  int *excludeNodataValues,
1814  int count,
1815  GDALDriverH *rtn_drv, int *destroy_rtn_drv
1816 ) {
1817  GDALDriverH drv = NULL;
1818  GDALDatasetH ds = NULL;
1819  double gt[6] = {0.0};
1820  CPLErr cplerr;
1821  GDALDataType gdal_pt = GDT_Unknown;
1822  GDALRasterBandH band;
1823  void *pVoid;
1824  char *pszDataPointer;
1825  char szGDALOption[50];
1826  char *apszOptions[4];
1827  double nodata = 0.0;
1828  int allocBandNums = 0;
1829  int allocNodataValues = 0;
1830 
1831  int i;
1832  int numBands;
1833  uint32_t width = 0;
1834  uint32_t height = 0;
1835  rt_band rtband = NULL;
1836  rt_pixtype pt = PT_END;
1837 
1838  assert(NULL != raster);
1839  assert(NULL != rtn_drv);
1840  assert(NULL != destroy_rtn_drv);
1841 
1842  *destroy_rtn_drv = 0;
1843 
1844  /* store raster in GDAL MEM raster */
1845  if (!rt_util_gdal_driver_registered("MEM")) {
1846  RASTER_DEBUG(4, "Registering MEM driver");
1847  GDALRegister_MEM();
1848  *destroy_rtn_drv = 1;
1849  }
1850  drv = GDALGetDriverByName("MEM");
1851  if (NULL == drv) {
1852  rterror("rt_raster_to_gdal_mem: Could not load the MEM GDAL driver");
1853  return 0;
1854  }
1855  *rtn_drv = drv;
1856 
1857  /* unload driver from GDAL driver manager */
1858  if (*destroy_rtn_drv) {
1859  RASTER_DEBUG(4, "Deregistering MEM driver");
1860  GDALDeregisterDriver(drv);
1861  }
1862 
1863  width = rt_raster_get_width(raster);
1864  height = rt_raster_get_height(raster);
1865  ds = GDALCreate(
1866  drv, "",
1867  width, height,
1868  0, GDT_Byte, NULL
1869  );
1870  if (NULL == ds) {
1871  rterror("rt_raster_to_gdal_mem: Could not create a GDALDataset to convert into");
1872  return 0;
1873  }
1874 
1875  /* add geotransform */
1876  rt_raster_get_geotransform_matrix(raster, gt);
1877  cplerr = GDALSetGeoTransform(ds, gt);
1878  if (cplerr != CE_None) {
1879  rterror("rt_raster_to_gdal_mem: Could not set geotransformation");
1880  GDALClose(ds);
1881  return 0;
1882  }
1883 
1884  /* set spatial reference */
1885  if (NULL != srs && strlen(srs)) {
1886  char *_srs = rt_util_gdal_convert_sr(srs, 0);
1887  if (_srs == NULL) {
1888  rterror("rt_raster_to_gdal_mem: Could not convert srs to GDAL accepted format");
1889  GDALClose(ds);
1890  return 0;
1891  }
1892 
1893  cplerr = GDALSetProjection(ds, _srs);
1894  CPLFree(_srs);
1895  if (cplerr != CE_None) {
1896  rterror("rt_raster_to_gdal_mem: Could not set projection");
1897  GDALClose(ds);
1898  return 0;
1899  }
1900  RASTER_DEBUGF(3, "Projection set to: %s", GDALGetProjectionRef(ds));
1901  }
1902 
1903  /* process bandNums */
1904  numBands = rt_raster_get_num_bands(raster);
1905  if (NULL != bandNums && count > 0) {
1906  for (i = 0; i < count; i++) {
1907  if (bandNums[i] >= numBands) {
1908  rterror("rt_raster_to_gdal_mem: The band index %d is invalid", bandNums[i]);
1909  GDALClose(ds);
1910  return 0;
1911  }
1912  }
1913  }
1914  else {
1915  count = numBands;
1916  bandNums = (uint32_t *) rtalloc(sizeof(uint32_t) * count);
1917  if (NULL == bandNums) {
1918  rterror("rt_raster_to_gdal_mem: Could not allocate memory for band indices");
1919  GDALClose(ds);
1920  return 0;
1921  }
1922  allocBandNums = 1;
1923  for (i = 0; i < count; i++) bandNums[i] = i;
1924  }
1925 
1926  /* process exclude_nodata_values */
1927  if (NULL == excludeNodataValues) {
1928  excludeNodataValues = (int *) rtalloc(sizeof(int) * count);
1929  if (NULL == excludeNodataValues) {
1930  rterror("rt_raster_to_gdal_mem: Could not allocate memory for NODATA flags");
1931  GDALClose(ds);
1932  return 0;
1933  }
1934  allocNodataValues = 1;
1935  for (i = 0; i < count; i++) excludeNodataValues[i] = 1;
1936  }
1937 
1938  /* add band(s) */
1939  for (i = 0; i < count; i++) {
1940  rtband = rt_raster_get_band(raster, bandNums[i]);
1941  if (NULL == rtband) {
1942  rterror("rt_raster_to_gdal_mem: Could not get requested band index %d", bandNums[i]);
1943  if (allocBandNums) rtdealloc(bandNums);
1944  if (allocNodataValues) rtdealloc(excludeNodataValues);
1945  GDALClose(ds);
1946  return 0;
1947  }
1948 
1949  pt = rt_band_get_pixtype(rtband);
1950  gdal_pt = rt_util_pixtype_to_gdal_datatype(pt);
1951  if (gdal_pt == GDT_Unknown)
1952  rtwarn("rt_raster_to_gdal_mem: Unknown pixel type for band");
1953 
1954  /*
1955  For all pixel types other than PT_8BSI, set pointer to start of data
1956  */
1957  if (pt != PT_8BSI) {
1958  pVoid = rt_band_get_data(rtband);
1959  RASTER_DEBUGF(4, "Band data is at pos %p", pVoid);
1960 
1961  pszDataPointer = (char *) rtalloc(20 * sizeof (char));
1962  sprintf(pszDataPointer, "%p", pVoid);
1963  RASTER_DEBUGF(4, "rt_raster_to_gdal_mem: szDatapointer is %p",
1964  pszDataPointer);
1965 
1966  if (strncasecmp(pszDataPointer, "0x", 2) == 0)
1967  sprintf(szGDALOption, "DATAPOINTER=%s", pszDataPointer);
1968  else
1969  sprintf(szGDALOption, "DATAPOINTER=0x%s", pszDataPointer);
1970 
1971  RASTER_DEBUG(3, "Storing info for GDAL MEM raster band");
1972 
1973  apszOptions[0] = szGDALOption;
1974  apszOptions[1] = NULL; /* pixel offset, not needed */
1975  apszOptions[2] = NULL; /* line offset, not needed */
1976  apszOptions[3] = NULL;
1977 
1978  /* free */
1979  rtdealloc(pszDataPointer);
1980 
1981  /* add band */
1982  if (GDALAddBand(ds, gdal_pt, apszOptions) == CE_Failure) {
1983  rterror("rt_raster_to_gdal_mem: Could not add GDAL raster band");
1984  if (allocBandNums) rtdealloc(bandNums);
1985  GDALClose(ds);
1986  return 0;
1987  }
1988  }
1989  /*
1990  PT_8BSI is special as GDAL has no equivalent pixel type.
1991  Must convert 8BSI to 16BSI so create basic band
1992  */
1993  else {
1994  /* add band */
1995  if (GDALAddBand(ds, gdal_pt, NULL) == CE_Failure) {
1996  rterror("rt_raster_to_gdal_mem: Could not add GDAL raster band");
1997  if (allocBandNums) rtdealloc(bandNums);
1998  if (allocNodataValues) rtdealloc(excludeNodataValues);
1999  GDALClose(ds);
2000  return 0;
2001  }
2002  }
2003 
2004  /* check band count */
2005  if (GDALGetRasterCount(ds) != i + 1) {
2006  rterror("rt_raster_to_gdal_mem: Error creating GDAL MEM raster band");
2007  if (allocBandNums) rtdealloc(bandNums);
2008  if (allocNodataValues) rtdealloc(excludeNodataValues);
2009  GDALClose(ds);
2010  return 0;
2011  }
2012 
2013  /* get new band */
2014  band = NULL;
2015  band = GDALGetRasterBand(ds, i + 1);
2016  if (NULL == band) {
2017  rterror("rt_raster_to_gdal_mem: Could not get GDAL band for additional processing");
2018  if (allocBandNums) rtdealloc(bandNums);
2019  if (allocNodataValues) rtdealloc(excludeNodataValues);
2020  GDALClose(ds);
2021  return 0;
2022  }
2023 
2024  /* PT_8BSI requires manual setting of pixels */
2025  if (pt == PT_8BSI) {
2026  int nXBlocks, nYBlocks;
2027  int nXBlockSize, nYBlockSize;
2028  int iXBlock, iYBlock;
2029  int nXValid, nYValid;
2030  int iX, iY;
2031  int iXMax, iYMax;
2032 
2033  int x, y, z;
2034  uint32_t valueslen = 0;
2035  int16_t *values = NULL;
2036  double value = 0.;
2037 
2038  /* this makes use of GDAL's "natural" blocks */
2039  GDALGetBlockSize(band, &nXBlockSize, &nYBlockSize);
2040  nXBlocks = (width + nXBlockSize - 1) / nXBlockSize;
2041  nYBlocks = (height + nYBlockSize - 1) / nYBlockSize;
2042  RASTER_DEBUGF(4, "(nXBlockSize, nYBlockSize) = (%d, %d)", nXBlockSize, nYBlockSize);
2043  RASTER_DEBUGF(4, "(nXBlocks, nYBlocks) = (%d, %d)", nXBlocks, nYBlocks);
2044 
2045  /* length is for the desired pixel type */
2046  valueslen = rt_pixtype_size(PT_16BSI) * nXBlockSize * nYBlockSize;
2047  values = rtalloc(valueslen);
2048  if (NULL == values) {
2049  rterror("rt_raster_to_gdal_mem: Could not allocate memory for GDAL band pixel values");
2050  if (allocBandNums) rtdealloc(bandNums);
2051  if (allocNodataValues) rtdealloc(excludeNodataValues);
2052  GDALClose(ds);
2053  return 0;
2054  }
2055 
2056  for (iYBlock = 0; iYBlock < nYBlocks; iYBlock++) {
2057  for (iXBlock = 0; iXBlock < nXBlocks; iXBlock++) {
2058  memset(values, 0, valueslen);
2059 
2060  x = iXBlock * nXBlockSize;
2061  y = iYBlock * nYBlockSize;
2062  RASTER_DEBUGF(4, "(iXBlock, iYBlock) = (%d, %d)", iXBlock, iYBlock);
2063  RASTER_DEBUGF(4, "(x, y) = (%d, %d)", x, y);
2064 
2065  /* valid block width */
2066  if ((iXBlock + 1) * nXBlockSize > width)
2067  nXValid = width - (iXBlock * nXBlockSize);
2068  else
2069  nXValid = nXBlockSize;
2070 
2071  /* valid block height */
2072  if ((iYBlock + 1) * nYBlockSize > height)
2073  nYValid = height - (iYBlock * nYBlockSize);
2074  else
2075  nYValid = nYBlockSize;
2076 
2077  RASTER_DEBUGF(4, "(nXValid, nYValid) = (%d, %d)", nXValid, nYValid);
2078 
2079  /* convert 8BSI values to 16BSI */
2080  z = 0;
2081  iYMax = y + nYValid;
2082  iXMax = x + nXValid;
2083  for (iY = y; iY < iYMax; iY++) {
2084  for (iX = x; iX < iXMax; iX++) {
2085  if (rt_band_get_pixel(rtband, iX, iY, &value, NULL) != ES_NONE) {
2086  rterror("rt_raster_to_gdal_mem: Could not get pixel value to convert from 8BSI to 16BSI");
2087  rtdealloc(values);
2088  if (allocBandNums) rtdealloc(bandNums);
2089  if (allocNodataValues) rtdealloc(excludeNodataValues);
2090  GDALClose(ds);
2091  return 0;
2092  }
2093 
2094  values[z++] = rt_util_clamp_to_16BSI(value);
2095  }
2096  }
2097 
2098  /* burn values */
2099  if (GDALRasterIO(
2100  band, GF_Write,
2101  x, y,
2102  nXValid, nYValid,
2103  values, nXValid, nYValid,
2104  gdal_pt,
2105  0, 0
2106  ) != CE_None) {
2107  rterror("rt_raster_to_gdal_mem: Could not write converted 8BSI to 16BSI values to GDAL band");
2108  rtdealloc(values);
2109  if (allocBandNums) rtdealloc(bandNums);
2110  if (allocNodataValues) rtdealloc(excludeNodataValues);
2111  GDALClose(ds);
2112  return 0;
2113  }
2114  }
2115  }
2116 
2117  rtdealloc(values);
2118  }
2119 
2120  /* Add nodata value for band */
2121  if (rt_band_get_hasnodata_flag(rtband) != FALSE && excludeNodataValues[i]) {
2122  rt_band_get_nodata(rtband, &nodata);
2123  if (GDALSetRasterNoDataValue(band, nodata) != CE_None)
2124  rtwarn("rt_raster_to_gdal_mem: Could not set nodata value for band");
2125  RASTER_DEBUGF(3, "nodata value set to %f", GDALGetRasterNoDataValue(band, NULL));
2126  }
2127 
2128 #if POSTGIS_DEBUG_LEVEL > 3
2129  {
2130  GDALRasterBandH _grb = NULL;
2131  double _min;
2132  double _max;
2133  double _mean;
2134  double _stddev;
2135 
2136  _grb = GDALGetRasterBand(ds, i + 1);
2137  GDALComputeRasterStatistics(_grb, FALSE, &_min, &_max, &_mean, &_stddev, NULL, NULL);
2138  RASTER_DEBUGF(4, "GDAL Band %d stats: %f, %f, %f, %f", i + 1, _min, _max, _mean, _stddev);
2139  }
2140 #endif
2141 
2142  }
2143 
2144  /* necessary??? */
2145  GDALFlushCache(ds);
2146 
2147  if (allocBandNums) rtdealloc(bandNums);
2148  if (allocNodataValues) rtdealloc(excludeNodataValues);
2149 
2150  return ds;
2151 }
2152 
2153 /******************************************************************************
2154 * rt_raster_from_gdal_dataset()
2155 ******************************************************************************/
2156 
2164 rt_raster
2165 rt_raster_from_gdal_dataset(GDALDatasetH ds) {
2166  rt_raster rast = NULL;
2167  double gt[6] = {0};
2168  CPLErr cplerr;
2169  uint32_t width = 0;
2170  uint32_t height = 0;
2171  uint32_t numBands = 0;
2172  int i = 0;
2173  char *authname = NULL;
2174  char *authcode = NULL;
2175 
2176  GDALRasterBandH gdband = NULL;
2177  GDALDataType gdpixtype = GDT_Unknown;
2178  rt_band band;
2179  int32_t idx;
2180  rt_pixtype pt = PT_END;
2181  uint32_t ptlen = 0;
2182  int hasnodata = 0;
2183  double nodataval;
2184 
2185  int x;
2186  int y;
2187 
2188  int nXBlocks, nYBlocks;
2189  int nXBlockSize, nYBlockSize;
2190  int iXBlock, iYBlock;
2191  int nXValid, nYValid;
2192  int iY;
2193 
2194  uint8_t *values = NULL;
2195  uint32_t valueslen = 0;
2196  uint8_t *ptr = NULL;
2197 
2198  assert(NULL != ds);
2199 
2200  /* raster size */
2201  width = GDALGetRasterXSize(ds);
2202  height = GDALGetRasterYSize(ds);
2203  RASTER_DEBUGF(3, "Raster dimensions (width x height): %d x %d", width, height);
2204 
2205  /* create new raster */
2206  RASTER_DEBUG(3, "Creating new raster");
2207  rast = rt_raster_new(width, height);
2208  if (NULL == rast) {
2209  rterror("rt_raster_from_gdal_dataset: Out of memory allocating new raster");
2210  return NULL;
2211  }
2212  RASTER_DEBUGF(3, "Created raster dimensions (width x height): %d x %d", rast->width, rast->height);
2213 
2214  /* get raster attributes */
2215  cplerr = GDALGetGeoTransform(ds, gt);
2216  if (GDALGetGeoTransform(ds, gt) != CE_None) {
2217  RASTER_DEBUG(4, "Using default geotransform matrix (0, 1, 0, 0, 0, -1)");
2218  gt[0] = 0;
2219  gt[1] = 1;
2220  gt[2] = 0;
2221  gt[3] = 0;
2222  gt[4] = 0;
2223  gt[5] = -1;
2224  }
2225 
2226  /* apply raster attributes */
2227  rt_raster_set_geotransform_matrix(rast, gt);
2228 
2229  RASTER_DEBUGF(3, "Raster geotransform (%f, %f, %f, %f, %f, %f)",
2230  gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]);
2231 
2232  /* srid */
2233  if (rt_util_gdal_sr_auth_info(ds, &authname, &authcode) == ES_NONE) {
2234  if (
2235  authname != NULL &&
2236  strcmp(authname, "EPSG") == 0 &&
2237  authcode != NULL
2238  ) {
2239  rt_raster_set_srid(rast, atoi(authcode));
2240  RASTER_DEBUGF(3, "New raster's SRID = %d", rast->srid);
2241  }
2242 
2243  if (authname != NULL)
2244  rtdealloc(authname);
2245  if (authcode != NULL)
2246  rtdealloc(authcode);
2247  }
2248 
2249  numBands = GDALGetRasterCount(ds);
2250 
2251 #if POSTGIS_DEBUG_LEVEL > 3
2252  for (i = 1; i <= numBands; i++) {
2253  GDALRasterBandH _grb = NULL;
2254  double _min;
2255  double _max;
2256  double _mean;
2257  double _stddev;
2258 
2259  _grb = GDALGetRasterBand(ds, i);
2260  GDALComputeRasterStatistics(_grb, FALSE, &_min, &_max, &_mean, &_stddev, NULL, NULL);
2261  RASTER_DEBUGF(4, "GDAL Band %d stats: %f, %f, %f, %f", i, _min, _max, _mean, _stddev);
2262  }
2263 #endif
2264 
2265  /* copy bands */
2266  for (i = 1; i <= numBands; i++) {
2267  RASTER_DEBUGF(3, "Processing band %d of %d", i, numBands);
2268  gdband = NULL;
2269  gdband = GDALGetRasterBand(ds, i);
2270  if (NULL == gdband) {
2271  rterror("rt_raster_from_gdal_dataset: Could not get GDAL band");
2272  rt_raster_destroy(rast);
2273  return NULL;
2274  }
2275  RASTER_DEBUGF(4, "gdband @ %p", gdband);
2276 
2277  /* pixtype */
2278  gdpixtype = GDALGetRasterDataType(gdband);
2279  RASTER_DEBUGF(4, "gdpixtype, size = %s, %d", GDALGetDataTypeName(gdpixtype), GDALGetDataTypeSize(gdpixtype) / 8);
2280  pt = rt_util_gdal_datatype_to_pixtype(gdpixtype);
2281  if (pt == PT_END) {
2282  rterror("rt_raster_from_gdal_dataset: Unknown pixel type for GDAL band");
2283  rt_raster_destroy(rast);
2284  return NULL;
2285  }
2286  ptlen = rt_pixtype_size(pt);
2287 
2288  /* size: width and height */
2289  width = GDALGetRasterBandXSize(gdband);
2290  height = GDALGetRasterBandYSize(gdband);
2291  RASTER_DEBUGF(3, "GDAL band dimensions (width x height): %d x %d", width, height);
2292 
2293  /* nodata */
2294  nodataval = GDALGetRasterNoDataValue(gdband, &hasnodata);
2295  RASTER_DEBUGF(3, "(hasnodata, nodataval) = (%d, %f)", hasnodata, nodataval);
2296 
2297  /* create band object */
2298  idx = rt_raster_generate_new_band(
2299  rast, pt,
2300  (hasnodata ? nodataval : 0),
2301  hasnodata, nodataval, rt_raster_get_num_bands(rast)
2302  );
2303  if (idx < 0) {
2304  rterror("rt_raster_from_gdal_dataset: Could not allocate memory for raster band");
2305  rt_raster_destroy(rast);
2306  return NULL;
2307  }
2308  band = rt_raster_get_band(rast, idx);
2309  RASTER_DEBUGF(3, "Created band of dimension (width x height): %d x %d", band->width, band->height);
2310 
2311  /* this makes use of GDAL's "natural" blocks */
2312  GDALGetBlockSize(gdband, &nXBlockSize, &nYBlockSize);
2313  nXBlocks = (width + nXBlockSize - 1) / nXBlockSize;
2314  nYBlocks = (height + nYBlockSize - 1) / nYBlockSize;
2315  RASTER_DEBUGF(4, "(nXBlockSize, nYBlockSize) = (%d, %d)", nXBlockSize, nYBlockSize);
2316  RASTER_DEBUGF(4, "(nXBlocks, nYBlocks) = (%d, %d)", nXBlocks, nYBlocks);
2317 
2318  /* allocate memory for values */
2319  valueslen = ptlen * nXBlockSize * nYBlockSize;
2320  values = rtalloc(valueslen);
2321  if (values == NULL) {
2322  rterror("rt_raster_from_gdal_dataset: Could not allocate memory for GDAL band pixel values");
2323  rt_raster_destroy(rast);
2324  return NULL;
2325  }
2326  RASTER_DEBUGF(3, "values @ %p of length = %d", values, valueslen);
2327 
2328  for (iYBlock = 0; iYBlock < nYBlocks; iYBlock++) {
2329  for (iXBlock = 0; iXBlock < nXBlocks; iXBlock++) {
2330  x = iXBlock * nXBlockSize;
2331  y = iYBlock * nYBlockSize;
2332  RASTER_DEBUGF(4, "(iXBlock, iYBlock) = (%d, %d)", iXBlock, iYBlock);
2333  RASTER_DEBUGF(4, "(x, y) = (%d, %d)", x, y);
2334 
2335  memset(values, 0, valueslen);
2336 
2337  /* valid block width */
2338  if ((iXBlock + 1) * nXBlockSize > width)
2339  nXValid = width - (iXBlock * nXBlockSize);
2340  else
2341  nXValid = nXBlockSize;
2342 
2343  /* valid block height */
2344  if ((iYBlock + 1) * nYBlockSize > height)
2345  nYValid = height - (iYBlock * nYBlockSize);
2346  else
2347  nYValid = nYBlockSize;
2348 
2349  RASTER_DEBUGF(4, "(nXValid, nYValid) = (%d, %d)", nXValid, nYValid);
2350 
2351  cplerr = GDALRasterIO(
2352  gdband, GF_Read,
2353  x, y,
2354  nXValid, nYValid,
2355  values, nXValid, nYValid,
2356  gdpixtype,
2357  0, 0
2358  );
2359  if (cplerr != CE_None) {
2360  rterror("rt_raster_from_gdal_dataset: Could not get data from GDAL raster");
2361  rtdealloc(values);
2362  rt_raster_destroy(rast);
2363  return NULL;
2364  }
2365 
2366  /* if block width is same as raster width, shortcut */
2367  if (nXBlocks == 1 && nYBlockSize > 1 && nXValid == width) {
2368  x = 0;
2369  y = nYBlockSize * iYBlock;
2370 
2371  RASTER_DEBUGF(4, "Setting set of pixel lines at (%d, %d) for %d pixels", x, y, nXValid * nYValid);
2372  rt_band_set_pixel_line(band, x, y, values, nXValid * nYValid);
2373  }
2374  else {
2375  ptr = values;
2376  x = nXBlockSize * iXBlock;
2377  for (iY = 0; iY < nYValid; iY++) {
2378  y = iY + (nYBlockSize * iYBlock);
2379 
2380  RASTER_DEBUGF(4, "Setting pixel line at (%d, %d) for %d pixels", x, y, nXValid);
2381  rt_band_set_pixel_line(band, x, y, ptr, nXValid);
2382  ptr += (nXValid * ptlen);
2383  }
2384  }
2385  }
2386  }
2387 
2388  /* free memory */
2389  rtdealloc(values);
2390  }
2391 
2392  return rast;
2393 }
2394 
2395 /******************************************************************************
2396 * rt_raster_gdal_rasterize()
2397 ******************************************************************************/
2398 
2399 typedef struct _rti_rasterize_arg_t* _rti_rasterize_arg;
2400 struct _rti_rasterize_arg_t {
2401  uint8_t noband;
2402 
2403  uint32_t numbands;
2404 
2405  OGRSpatialReferenceH src_sr;
2406 
2407  rt_pixtype *pixtype;
2408  double *init;
2409  double *nodata;
2410  uint8_t *hasnodata;
2411  double *value;
2412  int *bandlist;
2413 };
2414 
2415 static _rti_rasterize_arg
2416 _rti_rasterize_arg_init() {
2417  _rti_rasterize_arg arg = NULL;
2418 
2419  arg = rtalloc(sizeof(struct _rti_rasterize_arg_t));
2420  if (arg == NULL) {
2421  rterror("_rti_rasterize_arg_init: Could not allocate memory for _rti_rasterize_arg");
2422  return NULL;
2423  }
2424 
2425  arg->noband = 0;
2426 
2427  arg->numbands = 0;
2428 
2429  arg->src_sr = NULL;
2430 
2431  arg->pixtype = NULL;
2432  arg->init = NULL;
2433  arg->nodata = NULL;
2434  arg->hasnodata = NULL;
2435  arg->value = NULL;
2436  arg->bandlist = NULL;
2437 
2438  return arg;
2439 }
2440 
2441 static void
2442 _rti_rasterize_arg_destroy(_rti_rasterize_arg arg) {
2443  if (arg->noband) {
2444  if (arg->pixtype != NULL)
2445  rtdealloc(arg->pixtype);
2446  if (arg->init != NULL)
2447  rtdealloc(arg->init);
2448  if (arg->nodata != NULL)
2449  rtdealloc(arg->nodata);
2450  if (arg->hasnodata != NULL)
2451  rtdealloc(arg->hasnodata);
2452  if (arg->value != NULL)
2453  rtdealloc(arg->value);
2454  }
2455 
2456  if (arg->bandlist != NULL)
2457  rtdealloc(arg->bandlist);
2458 
2459  if (arg->src_sr != NULL)
2460  OSRDestroySpatialReference(arg->src_sr);
2461 
2462  rtdealloc(arg);
2463 }
2464 
2491 rt_raster
2492 rt_raster_gdal_rasterize(
2493  const unsigned char *wkb, uint32_t wkb_len,
2494  const char *srs,
2495  uint32_t num_bands, rt_pixtype *pixtype,
2496  double *init, double *value,
2497  double *nodata, uint8_t *hasnodata,
2498  int *width, int *height,
2499  double *scale_x, double *scale_y,
2500  double *ul_xw, double *ul_yw,
2501  double *grid_xw, double *grid_yw,
2502  double *skew_x, double *skew_y,
2503  char **options
2504 ) {
2505  rt_raster rast = NULL;
2506  int i = 0;
2507  int err = 0;
2508 
2509  _rti_rasterize_arg arg = NULL;
2510 
2511  int _dim[2] = {0};
2512  double _scale[2] = {0};
2513  double _skew[2] = {0};
2514 
2515  OGRErr ogrerr;
2516  OGRGeometryH src_geom;
2517  OGREnvelope src_env;
2518  rt_envelope extent;
2519  OGRwkbGeometryType wkbtype = wkbUnknown;
2520 
2521  int ul_user = 0;
2522 
2523  CPLErr cplerr;
2524  double _gt[6] = {0};
2525  GDALDriverH _drv = NULL;
2526  int unload_drv = 0;
2527  GDALDatasetH _ds = NULL;
2528  GDALRasterBandH _band = NULL;
2529 
2530  uint16_t _width = 0;
2531  uint16_t _height = 0;
2532 
2533  RASTER_DEBUG(3, "starting");
2534 
2535  assert(NULL != wkb);
2536  assert(0 != wkb_len);
2537 
2538  /* internal variables */
2539  arg = _rti_rasterize_arg_init();
2540  if (arg == NULL) {
2541  rterror("rt_raster_gdal_rasterize: Could not initialize internal variables");
2542  return NULL;
2543  }
2544 
2545  /* no bands, raster is a mask */
2546  if (num_bands < 1) {
2547  arg->noband = 1;
2548  arg->numbands = 1;
2549 
2550  arg->pixtype = (rt_pixtype *) rtalloc(sizeof(rt_pixtype));
2551  arg->pixtype[0] = PT_8BUI;
2552 
2553  arg->init = (double *) rtalloc(sizeof(double));
2554  arg->init[0] = 0;
2555 
2556  arg->nodata = (double *) rtalloc(sizeof(double));
2557  arg->nodata[0] = 0;
2558 
2559  arg->hasnodata = (uint8_t *) rtalloc(sizeof(uint8_t));
2560  arg->hasnodata[0] = 1;
2561 
2562  arg->value = (double *) rtalloc(sizeof(double));
2563  arg->value[0] = 1;
2564  }
2565  else {
2566  arg->noband = 0;
2567  arg->numbands = num_bands;
2568 
2569  arg->pixtype = pixtype;
2570  arg->init = init;
2571  arg->nodata = nodata;
2572  arg->hasnodata = hasnodata;
2573  arg->value = value;
2574  }
2575 
2576  /* OGR spatial reference */
2577  if (NULL != srs && strlen(srs)) {
2578  arg->src_sr = OSRNewSpatialReference(NULL);
2579  if (OSRSetFromUserInput(arg->src_sr, srs) != OGRERR_NONE) {
2580  rterror("rt_raster_gdal_rasterize: Could not create OSR spatial reference using the provided srs: %s", srs);
2581  _rti_rasterize_arg_destroy(arg);
2582  return NULL;
2583  }
2584  }
2585 
2586  /* convert WKB to OGR Geometry */
2587  ogrerr = OGR_G_CreateFromWkb((unsigned char *) wkb, arg->src_sr, &src_geom, wkb_len);
2588  if (ogrerr != OGRERR_NONE) {
2589  rterror("rt_raster_gdal_rasterize: Could not create OGR Geometry from WKB");
2590 
2591  _rti_rasterize_arg_destroy(arg);
2592  /* OGRCleanupAll(); */
2593 
2594  return NULL;
2595  }
2596 
2597  /* OGR Geometry is empty */
2598  if (OGR_G_IsEmpty(src_geom)) {
2599  rtinfo("Geometry provided is empty. Returning empty raster");
2600 
2601  OGR_G_DestroyGeometry(src_geom);
2602  _rti_rasterize_arg_destroy(arg);
2603  /* OGRCleanupAll(); */
2604 
2605  return rt_raster_new(0, 0);
2606  }
2607 
2608  /* get envelope */
2609  OGR_G_GetEnvelope(src_geom, &src_env);
2610  rt_util_from_ogr_envelope(src_env, &extent);
2611 
2612  RASTER_DEBUGF(3, "Suggested raster envelope: %f, %f, %f, %f",
2613  extent.MinX, extent.MinY, extent.MaxX, extent.MaxY);
2614 
2615  /* user-defined scale */
2616  if (
2617  (NULL != scale_x) &&
2618  (NULL != scale_y) &&
2619  (FLT_NEQ(*scale_x, 0.0)) &&
2620  (FLT_NEQ(*scale_y, 0.0))
2621  ) {
2622  /* for now, force scale to be in left-right, top-down orientation */
2623  _scale[0] = fabs(*scale_x);
2624  _scale[1] = fabs(*scale_y);
2625  }
2626  /* user-defined width/height */
2627  else if (
2628  (NULL != width) &&
2629  (NULL != height) &&
2630  (FLT_NEQ(*width, 0.0)) &&
2631  (FLT_NEQ(*height, 0.0))
2632  ) {
2633  _dim[0] = abs(*width);
2634  _dim[1] = abs(*height);
2635 
2636  if (FLT_NEQ(extent.MaxX, extent.MinX))
2637  _scale[0] = fabs((extent.MaxX - extent.MinX) / _dim[0]);
2638  else
2639  _scale[0] = 1.;
2640 
2641  if (FLT_NEQ(extent.MaxY, extent.MinY))
2642  _scale[1] = fabs((extent.MaxY - extent.MinY) / _dim[1]);
2643  else
2644  _scale[1] = 1.;
2645  }
2646  else {
2647  rterror("rt_raster_gdal_rasterize: Values must be provided for width and height or X and Y of scale");
2648 
2649  OGR_G_DestroyGeometry(src_geom);
2650  _rti_rasterize_arg_destroy(arg);
2651  /* OGRCleanupAll(); */
2652 
2653  return NULL;
2654  }
2655  RASTER_DEBUGF(3, "scale (x, y) = %f, %f", _scale[0], -1 * _scale[1]);
2656  RASTER_DEBUGF(3, "dim (x, y) = %d, %d", _dim[0], _dim[1]);
2657 
2658  /* user-defined skew */
2659  if (NULL != skew_x) {
2660  _skew[0] = *skew_x;
2661 
2662  /*
2663  negative scale-x affects skew
2664  for now, force skew to be in left-right, top-down orientation
2665  */
2666  if (
2667  NULL != scale_x &&
2668  *scale_x < 0.
2669  ) {
2670  _skew[0] *= -1;
2671  }
2672  }
2673  if (NULL != skew_y) {
2674  _skew[1] = *skew_y;
2675 
2676  /*
2677  positive scale-y affects skew
2678  for now, force skew to be in left-right, top-down orientation
2679  */
2680  if (
2681  NULL != scale_y &&
2682  *scale_y > 0.
2683  ) {
2684  _skew[1] *= -1;
2685  }
2686  }
2687 
2688  /*
2689  if geometry is a point, a linestring or set of either and bounds not set,
2690  increase extent by a pixel to avoid missing points on border
2691 
2692  a whole pixel is used instead of half-pixel due to backward
2693  compatibility with GDAL 1.6, 1.7 and 1.8. 1.9+ works fine with half-pixel.
2694  */
2695  wkbtype = wkbFlatten(OGR_G_GetGeometryType(src_geom));
2696  if ((
2697  (wkbtype == wkbPoint) ||
2698  (wkbtype == wkbMultiPoint) ||
2699  (wkbtype == wkbLineString) ||
2700  (wkbtype == wkbMultiLineString)
2701  ) &&
2702  _dim[0] == 0 &&
2703  _dim[1] == 0
2704  ) {
2705  int result;
2706  LWPOLY *epoly = NULL;
2707  LWGEOM *lwgeom = NULL;
2708  GEOSGeometry *egeom = NULL;
2709  GEOSGeometry *geom = NULL;
2710 
2711  RASTER_DEBUG(3, "Testing geometry is properly contained by extent");
2712 
2713  /*
2714  see if geometry is properly contained by extent
2715  all parts of geometry lies within extent
2716  */
2717 
2718  /* initialize GEOS */
2719  initGEOS(rtinfo, lwgeom_geos_error);
2720 
2721  /* convert envelope to geometry */
2722  RASTER_DEBUG(4, "Converting envelope to geometry");
2723  epoly = rt_util_envelope_to_lwpoly(extent);
2724  if (epoly == NULL) {
2725  rterror("rt_raster_gdal_rasterize: Could not create envelope's geometry to test if geometry is properly contained by extent");
2726 
2727  OGR_G_DestroyGeometry(src_geom);
2728  _rti_rasterize_arg_destroy(arg);
2729  /* OGRCleanupAll(); */
2730 
2731  return NULL;
2732  }
2733 
2734  egeom = (GEOSGeometry *) LWGEOM2GEOS(lwpoly_as_lwgeom(epoly), 0);
2735  lwpoly_free(epoly);
2736 
2737  /* convert WKB to geometry */
2738  RASTER_DEBUG(4, "Converting WKB to geometry");
2739  lwgeom = lwgeom_from_wkb(wkb, wkb_len, LW_PARSER_CHECK_NONE);
2740  geom = (GEOSGeometry *) LWGEOM2GEOS(lwgeom, 0);
2741  lwgeom_free(lwgeom);
2742 
2743  result = GEOSRelatePattern(egeom, geom, "T**FF*FF*");
2744  GEOSGeom_destroy(geom);
2745  GEOSGeom_destroy(egeom);
2746 
2747  if (result == 2) {
2748  rterror("rt_raster_gdal_rasterize: Could not test if geometry is properly contained by extent for geometry within extent");
2749 
2750  OGR_G_DestroyGeometry(src_geom);
2751  _rti_rasterize_arg_destroy(arg);
2752  /* OGRCleanupAll(); */
2753 
2754  return NULL;
2755  }
2756 
2757  /* geometry NOT properly contained by extent */
2758  if (!result) {
2759 
2760 #if POSTGIS_GDAL_VERSION > 18
2761 
2762  /* check alignment flag: grid_xw */
2763  if (
2764  (NULL == ul_xw && NULL == ul_yw) &&
2765  (NULL != grid_xw && NULL != grid_xw) &&
2766  FLT_NEQ(*grid_xw, extent.MinX)
2767  ) {
2768  /* do nothing */
2769  RASTER_DEBUG(3, "Skipping extent adjustment on X-axis due to upcoming alignment");
2770  }
2771  else {
2772  RASTER_DEBUG(3, "Adjusting extent for GDAL > 1.8 by half the scale on X-axis");
2773  extent.MinX -= (_scale[0] / 2.);
2774  extent.MaxX += (_scale[0] / 2.);
2775  }
2776 
2777  /* check alignment flag: grid_yw */
2778  if (
2779  (NULL == ul_xw && NULL == ul_yw) &&
2780  (NULL != grid_xw && NULL != grid_xw) &&
2781  FLT_NEQ(*grid_yw, extent.MaxY)
2782  ) {
2783  /* do nothing */
2784  RASTER_DEBUG(3, "Skipping extent adjustment on Y-axis due to upcoming alignment");
2785  }
2786  else {
2787  RASTER_DEBUG(3, "Adjusting extent for GDAL > 1.8 by half the scale on Y-axis");
2788  extent.MinY -= (_scale[1] / 2.);
2789  extent.MaxY += (_scale[1] / 2.);
2790  }
2791 
2792 #else
2793 
2794  /* check alignment flag: grid_xw */
2795  if (
2796  (NULL == ul_xw && NULL == ul_yw) &&
2797  (NULL != grid_xw && NULL != grid_xw) &&
2798  FLT_NEQ(*grid_xw, extent.MinX)
2799  ) {
2800  /* do nothing */
2801  RASTER_DEBUG(3, "Skipping extent adjustment on X-axis due to upcoming alignment");
2802  }
2803  else {
2804  RASTER_DEBUG(3, "Adjusting extent for GDAL <= 1.8 by the scale on X-axis");
2805  extent.MinX -= _scale[0];
2806  extent.MaxX += _scale[0];
2807  }
2808 
2809 
2810  /* check alignment flag: grid_yw */
2811  if (
2812  (NULL == ul_xw && NULL == ul_yw) &&
2813  (NULL != grid_xw && NULL != grid_xw) &&
2814  FLT_NEQ(*grid_yw, extent.MaxY)
2815  ) {
2816  /* do nothing */
2817  RASTER_DEBUG(3, "Skipping extent adjustment on Y-axis due to upcoming alignment");
2818  }
2819  else {
2820  RASTER_DEBUG(3, "Adjusting extent for GDAL <= 1.8 by the scale on Y-axis");
2821  extent.MinY -= _scale[1];
2822  extent.MaxY += _scale[1];
2823  }
2824 
2825 #endif
2826 
2827  }
2828 
2829  RASTER_DEBUGF(3, "Adjusted extent: %f, %f, %f, %f",
2830  extent.MinX, extent.MinY, extent.MaxX, extent.MaxY);
2831 
2832  extent.UpperLeftX = extent.MinX;
2833  extent.UpperLeftY = extent.MaxY;
2834  }
2835 
2836  /* reprocess extent if skewed */
2837  if (
2838  FLT_NEQ(_skew[0], 0) ||
2839  FLT_NEQ(_skew[1], 0)
2840  ) {
2841  rt_raster skewedrast;
2842 
2843  RASTER_DEBUG(3, "Computing skewed extent's envelope");
2844 
2845  skewedrast = rt_raster_compute_skewed_raster(
2846  extent,
2847  _skew,
2848  _scale,
2849  0.01
2850  );
2851  if (skewedrast == NULL) {
2852  rterror("rt_raster_gdal_rasterize: Could not compute skewed raster");
2853 
2854  OGR_G_DestroyGeometry(src_geom);
2855  _rti_rasterize_arg_destroy(arg);
2856  /* OGRCleanupAll(); */
2857 
2858  return NULL;
2859  }
2860 
2861  _dim[0] = skewedrast->width;
2862  _dim[1] = skewedrast->height;
2863 
2864  extent.UpperLeftX = skewedrast->ipX;
2865  extent.UpperLeftY = skewedrast->ipY;
2866 
2867  rt_raster_destroy(skewedrast);
2868  }
2869 
2870  /* raster dimensions */
2871  if (!_dim[0])
2872  _dim[0] = (int) fmax((fabs(extent.MaxX - extent.MinX) + (_scale[0] / 2.)) / _scale[0], 1);
2873  if (!_dim[1])
2874  _dim[1] = (int) fmax((fabs(extent.MaxY - extent.MinY) + (_scale[1] / 2.)) / _scale[1], 1);
2875 
2876  /* temporary raster */
2877  rast = rt_raster_new(_dim[0], _dim[1]);
2878  if (rast == NULL) {
2879  rterror("rt_raster_gdal_rasterize: Out of memory allocating temporary raster");
2880 
2881  OGR_G_DestroyGeometry(src_geom);
2882  _rti_rasterize_arg_destroy(arg);
2883  /* OGRCleanupAll(); */
2884 
2885  return NULL;
2886  }
2887 
2888  /* set raster's spatial attributes */
2889  rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
2890  rt_raster_set_scale(rast, _scale[0], -1 * _scale[1]);
2891  rt_raster_set_skews(rast, _skew[0], _skew[1]);
2892 
2893  rt_raster_get_geotransform_matrix(rast, _gt);
2894  RASTER_DEBUGF(3, "Temp raster's geotransform: %f, %f, %f, %f, %f, %f",
2895  _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
2896  RASTER_DEBUGF(3, "Temp raster's dimensions (width x height): %d x %d",
2897  _dim[0], _dim[1]);
2898 
2899  /* user-specified upper-left corner */
2900  if (
2901  NULL != ul_xw &&
2902  NULL != ul_yw
2903  ) {
2904  ul_user = 1;
2905 
2906  RASTER_DEBUGF(4, "Using user-specified upper-left corner: %f, %f", *ul_xw, *ul_yw);
2907 
2908  /* set upper-left corner */
2909  rt_raster_set_offsets(rast, *ul_xw, *ul_yw);
2910  extent.UpperLeftX = *ul_xw;
2911  extent.UpperLeftY = *ul_yw;
2912  }
2913  else if (
2914  ((NULL != ul_xw) && (NULL == ul_yw)) ||
2915  ((NULL == ul_xw) && (NULL != ul_yw))
2916  ) {
2917  rterror("rt_raster_gdal_rasterize: Both X and Y upper-left corner values must be provided");
2918 
2919  rt_raster_destroy(rast);
2920  OGR_G_DestroyGeometry(src_geom);
2921  _rti_rasterize_arg_destroy(arg);
2922  /* OGRCleanupAll(); */
2923 
2924  return NULL;
2925  }
2926 
2927  /* alignment only considered if upper-left corner not provided */
2928  if (
2929  !ul_user && (
2930  (NULL != grid_xw) || (NULL != grid_yw)
2931  )
2932  ) {
2933 
2934  if (
2935  ((NULL != grid_xw) && (NULL == grid_yw)) ||
2936  ((NULL == grid_xw) && (NULL != grid_yw))
2937  ) {
2938  rterror("rt_raster_gdal_rasterize: Both X and Y alignment values must be provided");
2939 
2940  rt_raster_destroy(rast);
2941  OGR_G_DestroyGeometry(src_geom);
2942  _rti_rasterize_arg_destroy(arg);
2943  /* OGRCleanupAll(); */
2944 
2945  return NULL;
2946  }
2947 
2948  RASTER_DEBUGF(4, "Aligning extent to user-specified grid: %f, %f", *grid_xw, *grid_yw);
2949 
2950  do {
2951  double _r[2] = {0};
2952  double _w[2] = {0};
2953 
2954  /* raster is already aligned */
2955  if (FLT_EQ(*grid_xw, extent.UpperLeftX) && FLT_EQ(*grid_yw, extent.UpperLeftY)) {
2956  RASTER_DEBUG(3, "Skipping raster alignment as it is already aligned to grid");
2957  break;
2958  }
2959 
2960  extent.UpperLeftX = rast->ipX;
2961  extent.UpperLeftY = rast->ipY;
2962  rt_raster_set_offsets(rast, *grid_xw, *grid_yw);
2963 
2964  /* process upper-left corner */
2965  if (rt_raster_geopoint_to_cell(
2966  rast,
2967  extent.UpperLeftX, extent.UpperLeftY,
2968  &(_r[0]), &(_r[1]),
2969  NULL
2970  ) != ES_NONE) {
2971  rterror("rt_raster_gdal_rasterize: Could not compute raster pixel for spatial coordinates");
2972 
2973  rt_raster_destroy(rast);
2974  OGR_G_DestroyGeometry(src_geom);
2975  _rti_rasterize_arg_destroy(arg);
2976  /* OGRCleanupAll(); */
2977 
2978  return NULL;
2979  }
2980 
2981  if (rt_raster_cell_to_geopoint(
2982  rast,
2983  _r[0], _r[1],
2984  &(_w[0]), &(_w[1]),
2985  NULL
2986  ) != ES_NONE) {
2987  rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
2988 
2989  rt_raster_destroy(rast);
2990  OGR_G_DestroyGeometry(src_geom);
2991  _rti_rasterize_arg_destroy(arg);
2992  /* OGRCleanupAll(); */
2993 
2994  return NULL;
2995  }
2996 
2997  /* shift occurred */
2998  if (FLT_NEQ(_w[0], extent.UpperLeftX)) {
2999  if (NULL == width)
3000  rast->width++;
3001  else if (NULL == scale_x) {
3002  double _c[2] = {0};
3003 
3004  rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
3005 
3006  /* get upper-right corner */
3007  if (rt_raster_cell_to_geopoint(
3008  rast,
3009  rast->width, 0,
3010  &(_c[0]), &(_c[1]),
3011  NULL
3012  ) != ES_NONE) {
3013  rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
3014 
3015  rt_raster_destroy(rast);
3016  OGR_G_DestroyGeometry(src_geom);
3017  _rti_rasterize_arg_destroy(arg);
3018  /* OGRCleanupAll(); */
3019 
3020  return NULL;
3021  }
3022 
3023  rast->scaleX = fabs((_c[0] - _w[0]) / ((double) rast->width));
3024  }
3025  }
3026  if (FLT_NEQ(_w[1], extent.UpperLeftY)) {
3027  if (NULL == height)
3028  rast->height++;
3029  else if (NULL == scale_y) {
3030  double _c[2] = {0};
3031 
3032  rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
3033 
3034  /* get upper-right corner */
3035  if (rt_raster_cell_to_geopoint(
3036  rast,
3037  0, rast->height,
3038  &(_c[0]), &(_c[1]),
3039  NULL
3040  ) != ES_NONE) {
3041  rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
3042 
3043  rt_raster_destroy(rast);
3044  OGR_G_DestroyGeometry(src_geom);
3045  _rti_rasterize_arg_destroy(arg);
3046  /* OGRCleanupAll(); */
3047 
3048  return NULL;
3049  }
3050 
3051  rast->scaleY = -1 * fabs((_c[1] - _w[1]) / ((double) rast->height));
3052  }
3053  }
3054 
3055  rt_raster_set_offsets(rast, _w[0], _w[1]);
3056  }
3057  while (0);
3058  }
3059 
3060  /*
3061  after this point, rt_envelope extent is no longer used
3062  */
3063 
3064  /* get key attributes from rast */
3065  _dim[0] = rast->width;
3066  _dim[1] = rast->height;
3067  rt_raster_get_geotransform_matrix(rast, _gt);
3068 
3069  /* scale-x is negative or scale-y is positive */
3070  if ((
3071  (NULL != scale_x) && (*scale_x < 0.)
3072  ) || (
3073  (NULL != scale_y) && (*scale_y > 0)
3074  )) {
3075  double _w[2] = {0};
3076 
3077  /* negative scale-x */
3078  if (
3079  (NULL != scale_x) &&
3080  (*scale_x < 0.)
3081  ) {
3082  RASTER_DEBUG(3, "Processing negative scale-x");
3083 
3084  if (rt_raster_cell_to_geopoint(
3085  rast,
3086  _dim[0], 0,
3087  &(_w[0]), &(_w[1]),
3088  NULL
3089  ) != ES_NONE) {
3090  rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
3091 
3092  rt_raster_destroy(rast);
3093  OGR_G_DestroyGeometry(src_geom);
3094  _rti_rasterize_arg_destroy(arg);
3095  /* OGRCleanupAll(); */
3096 
3097  return NULL;
3098  }
3099 
3100  _gt[0] = _w[0];
3101  _gt[1] = *scale_x;
3102 
3103  /* check for skew */
3104  if (NULL != skew_x && FLT_NEQ(*skew_x, 0))
3105  _gt[2] = *skew_x;
3106  }
3107  /* positive scale-y */
3108  if (
3109  (NULL != scale_y) &&
3110  (*scale_y > 0)
3111  ) {
3112  RASTER_DEBUG(3, "Processing positive scale-y");
3113 
3114  if (rt_raster_cell_to_geopoint(
3115  rast,
3116  0, _dim[1],
3117  &(_w[0]), &(_w[1]),
3118  NULL
3119  ) != ES_NONE) {
3120  rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
3121 
3122  rt_raster_destroy(rast);
3123  OGR_G_DestroyGeometry(src_geom);
3124  _rti_rasterize_arg_destroy(arg);
3125  /* OGRCleanupAll(); */
3126 
3127  return NULL;
3128  }
3129 
3130  _gt[3] = _w[1];
3131  _gt[5] = *scale_y;
3132 
3133  /* check for skew */
3134  if (NULL != skew_y && FLT_NEQ(*skew_y, 0))
3135  _gt[4] = *skew_y;
3136  }
3137  }
3138 
3139  rt_raster_destroy(rast);
3140  rast = NULL;
3141 
3142  RASTER_DEBUGF(3, "Applied geotransform: %f, %f, %f, %f, %f, %f",
3143  _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
3144  RASTER_DEBUGF(3, "Raster dimensions (width x height): %d x %d",
3145  _dim[0], _dim[1]);
3146 
3147  /* load GDAL mem */
3148  if (!rt_util_gdal_driver_registered("MEM")) {
3149  RASTER_DEBUG(4, "Registering MEM driver");
3150  GDALRegister_MEM();
3151  unload_drv = 1;
3152  }
3153  _drv = GDALGetDriverByName("MEM");
3154  if (NULL == _drv) {
3155  rterror("rt_raster_gdal_rasterize: Could not load the MEM GDAL driver");
3156 
3157  OGR_G_DestroyGeometry(src_geom);
3158  _rti_rasterize_arg_destroy(arg);
3159  /* OGRCleanupAll(); */
3160 
3161  return NULL;
3162  }
3163 
3164  /* unload driver from GDAL driver manager */
3165  if (unload_drv) {
3166  RASTER_DEBUG(4, "Deregistering MEM driver");
3167  GDALDeregisterDriver(_drv);
3168  }
3169 
3170  _ds = GDALCreate(_drv, "", _dim[0], _dim[1], 0, GDT_Byte, NULL);
3171  if (NULL == _ds) {
3172  rterror("rt_raster_gdal_rasterize: Could not create a GDALDataset to rasterize the geometry into");
3173 
3174  OGR_G_DestroyGeometry(src_geom);
3175  _rti_rasterize_arg_destroy(arg);
3176  /* OGRCleanupAll(); */
3177  if (unload_drv) GDALDestroyDriver(_drv);
3178 
3179  return NULL;
3180  }
3181 
3182  /* set geotransform */
3183  cplerr = GDALSetGeoTransform(_ds, _gt);
3184  if (cplerr != CE_None) {
3185  rterror("rt_raster_gdal_rasterize: Could not set geotransform on GDALDataset");
3186 
3187  OGR_G_DestroyGeometry(src_geom);
3188  _rti_rasterize_arg_destroy(arg);
3189  /* OGRCleanupAll(); */
3190 
3191  GDALClose(_ds);
3192  if (unload_drv) GDALDestroyDriver(_drv);
3193 
3194  return NULL;
3195  }
3196 
3197  /* set SRS */
3198  if (NULL != arg->src_sr) {
3199  char *_srs = NULL;
3200  OSRExportToWkt(arg->src_sr, &_srs);
3201 
3202  cplerr = GDALSetProjection(_ds, _srs);
3203  CPLFree(_srs);
3204  if (cplerr != CE_None) {
3205  rterror("rt_raster_gdal_rasterize: Could not set projection on GDALDataset");
3206 
3207  OGR_G_DestroyGeometry(src_geom);
3208  _rti_rasterize_arg_destroy(arg);
3209  /* OGRCleanupAll(); */
3210 
3211  GDALClose(_ds);
3212  if (unload_drv) GDALDestroyDriver(_drv);
3213 
3214  return NULL;
3215  }
3216  }
3217 
3218  /* set bands */
3219  for (i = 0; i < arg->numbands; i++) {
3220  err = 0;
3221 
3222  do {
3223  /* add band */
3224  cplerr = GDALAddBand(_ds, rt_util_pixtype_to_gdal_datatype(arg->pixtype[i]), NULL);
3225  if (cplerr != CE_None) {
3226  rterror("rt_raster_gdal_rasterize: Could not add band to GDALDataset");
3227  err = 1;
3228  break;
3229  }
3230 
3231  _band = GDALGetRasterBand(_ds, i + 1);
3232  if (NULL == _band) {
3233  rterror("rt_raster_gdal_rasterize: Could not get band %d from GDALDataset", i + 1);
3234  err = 1;
3235  break;
3236  }
3237 
3238  /* nodata value */
3239  if (arg->hasnodata[i]) {
3240  RASTER_DEBUGF(4, "Setting NODATA value of band %d to %f", i, arg->nodata[i]);
3241  cplerr = GDALSetRasterNoDataValue(_band, arg->nodata[i]);
3242  if (cplerr != CE_None) {
3243  rterror("rt_raster_gdal_rasterize: Could not set nodata value");
3244  err = 1;
3245  break;
3246  }
3247  RASTER_DEBUGF(4, "NODATA value set to %f", GDALGetRasterNoDataValue(_band, NULL));
3248  }
3249 
3250  /* initial value */
3251  cplerr = GDALFillRaster(_band, arg->init[i], 0);
3252  if (cplerr != CE_None) {
3253  rterror("rt_raster_gdal_rasterize: Could not set initial value");
3254  err = 1;
3255  break;
3256  }
3257  }
3258  while (0);
3259 
3260  if (err) {
3261 
3262  OGR_G_DestroyGeometry(src_geom);
3263  _rti_rasterize_arg_destroy(arg);
3264 
3265  /* OGRCleanupAll(); */
3266 
3267  GDALClose(_ds);
3268  if (unload_drv) GDALDestroyDriver(_drv);
3269 
3270  return NULL;
3271  }
3272  }
3273 
3274  arg->bandlist = (int *) rtalloc(sizeof(int) * arg->numbands);
3275  for (i = 0; i < arg->numbands; i++) arg->bandlist[i] = i + 1;
3276 
3277  /* burn geometry */
3278  cplerr = GDALRasterizeGeometries(
3279  _ds,
3280  arg->numbands, arg->bandlist,
3281  1, &src_geom,
3282  NULL, NULL,
3283  arg->value,
3284  options,
3285  NULL, NULL
3286  );
3287  if (cplerr != CE_None) {
3288  rterror("rt_raster_gdal_rasterize: Could not rasterize geometry");
3289 
3290  OGR_G_DestroyGeometry(src_geom);
3291  _rti_rasterize_arg_destroy(arg);
3292  /* OGRCleanupAll(); */
3293 
3294  GDALClose(_ds);
3295  if (unload_drv) GDALDestroyDriver(_drv);
3296 
3297  return NULL;
3298  }
3299 
3300  /* convert gdal dataset to raster */
3301  GDALFlushCache(_ds);
3302  RASTER_DEBUG(3, "Converting GDAL dataset to raster");
3303  rast = rt_raster_from_gdal_dataset(_ds);
3304 
3305  OGR_G_DestroyGeometry(src_geom);
3306  /* OGRCleanupAll(); */
3307 
3308  GDALClose(_ds);
3309  if (unload_drv) GDALDestroyDriver(_drv);
3310 
3311  if (NULL == rast) {
3312  rterror("rt_raster_gdal_rasterize: Could not rasterize geometry");
3313  return NULL;
3314  }
3315 
3316  /* width, height */
3317  _width = rt_raster_get_width(rast);
3318  _height = rt_raster_get_height(rast);
3319 
3320  /* check each band for pixtype */
3321  for (i = 0; i < arg->numbands; i++) {
3322  uint8_t *data = NULL;
3323  rt_band band = NULL;
3324  rt_band oldband = NULL;
3325 
3326  double val = 0;
3327  int nodata = 0;
3328  int hasnodata = 0;
3329  double nodataval = 0;
3330  int x = 0;
3331  int y = 0;
3332 
3333  oldband = rt_raster_get_band(rast, i);
3334  if (oldband == NULL) {
3335  rterror("rt_raster_gdal_rasterize: Could not get band %d of output raster", i);
3336  _rti_rasterize_arg_destroy(arg);
3337  rt_raster_destroy(rast);
3338  return NULL;
3339  }
3340 
3341  /* band is of user-specified type */
3342  if (rt_band_get_pixtype(oldband) == arg->pixtype[i])
3343  continue;
3344 
3345  /* hasnodata, nodataval */
3346  hasnodata = rt_band_get_hasnodata_flag(oldband);
3347  if (hasnodata)
3348  rt_band_get_nodata(oldband, &nodataval);
3349 
3350  /* allocate data */
3351  data = rtalloc(rt_pixtype_size(arg->pixtype[i]) * _width * _height);
3352  if (data == NULL) {
3353  rterror("rt_raster_gdal_rasterize: Could not allocate memory for band data");
3354  _rti_rasterize_arg_destroy(arg);
3355  rt_raster_destroy(rast);
3356  return NULL;
3357  }
3358  memset(data, 0, rt_pixtype_size(arg->pixtype[i]) * _width * _height);
3359 
3360  /* create new band of correct type */
3361  band = rt_band_new_inline(
3362  _width, _height,
3363  arg->pixtype[i],
3364  hasnodata, nodataval,
3365  data
3366  );
3367  if (band == NULL) {
3368  rterror("rt_raster_gdal_rasterize: Could not create band");
3369  rtdealloc(data);
3370  _rti_rasterize_arg_destroy(arg);
3371  rt_raster_destroy(rast);
3372  return NULL;
3373  }
3374 
3375  /* give ownership of data to band */
3376  rt_band_set_ownsdata_flag(band, 1);
3377 
3378  /* copy pixel by pixel */
3379  for (x = 0; x < _width; x++) {
3380  for (y = 0; y < _height; y++) {
3381  err = rt_band_get_pixel(oldband, x, y, &val, &nodata);
3382  if (err != ES_NONE) {
3383  rterror("rt_raster_gdal_rasterize: Could not get pixel value");
3384  _rti_rasterize_arg_destroy(arg);
3385  rt_raster_destroy(rast);
3386  rt_band_destroy(band);
3387  return NULL;
3388  }
3389 
3390  if (nodata)
3391  val = nodataval;
3392 
3393  err = rt_band_set_pixel(band, x, y, val, NULL);
3394  if (err != ES_NONE) {
3395  rterror("rt_raster_gdal_rasterize: Could not set pixel value");
3396  _rti_rasterize_arg_destroy(arg);
3397  rt_raster_destroy(rast);
3398  rt_band_destroy(band);
3399  return NULL;
3400  }
3401  }
3402  }
3403 
3404  /* replace band */
3405  oldband = rt_raster_replace_band(rast, band, i);
3406  if (oldband == NULL) {
3407  rterror("rt_raster_gdal_rasterize: Could not replace band %d of output raster", i);
3408  _rti_rasterize_arg_destroy(arg);
3409  rt_raster_destroy(rast);
3410  rt_band_destroy(band);
3411  return NULL;
3412  }
3413 
3414  /* free oldband */
3415  rt_band_destroy(oldband);
3416  }
3417 
3418  _rti_rasterize_arg_destroy(arg);
3419 
3420  RASTER_DEBUG(3, "done");
3421 
3422  return rast;
3423 }
3424 
3425 /******************************************************************************
3426 * rt_raster_from_two_rasters()
3427 ******************************************************************************/
3428 
3429 /*
3430  * Return raster of computed extent specified extenttype applied
3431  * on two input rasters. The raster returned should be freed by
3432  * the caller
3433  *
3434  * @param rast1 : the first raster
3435  * @param rast2 : the second raster
3436  * @param extenttype : type of extent for the output raster
3437  * @param *rtnraster : raster of computed extent
3438  * @param *offset : 4-element array indicating the X,Y offsets
3439  * for each raster. 0,1 for rast1 X,Y. 2,3 for rast2 X,Y.
3440  *
3441  * @return ES_NONE if success, ES_ERROR if error
3442  */
3443 rt_errorstate
3444 rt_raster_from_two_rasters(
3445  rt_raster rast1, rt_raster rast2,
3446  rt_extenttype extenttype,
3447  rt_raster *rtnraster, double *offset
3448 ) {
3449  int i;
3450 
3451  rt_raster _rast[2] = {rast1, rast2};
3452  double _offset[2][4] = {{0.}};
3453  uint16_t _dim[2][2] = {{0}};
3454 
3455  rt_raster raster = NULL;
3456  int aligned = 0;
3457  int dim[2] = {0};
3458  double gt[6] = {0.};
3459 
3460  assert(NULL != rast1);
3461  assert(NULL != rast2);
3462  assert(NULL != rtnraster);
3463 
3464  /* set rtnraster to NULL */
3465  *rtnraster = NULL;
3466 
3467  /* rasters must be aligned */
3468  if (rt_raster_same_alignment(rast1, rast2, &aligned, NULL) != ES_NONE) {
3469  rterror("rt_raster_from_two_rasters: Could not test for alignment on the two rasters");
3470  return ES_ERROR;
3471  }
3472  if (!aligned) {
3473  rterror("rt_raster_from_two_rasters: The two rasters provided do not have the same alignment");
3474  return ES_ERROR;
3475  }
3476 
3477  /* dimensions */
3478  _dim[0][0] = rast1->width;
3479  _dim[0][1] = rast1->height;
3480  _dim[1][0] = rast2->width;
3481  _dim[1][1] = rast2->height;
3482 
3483  /* get raster offsets */
3484  if (rt_raster_geopoint_to_cell(
3485  _rast[1],
3486  _rast[0]->ipX, _rast[0]->ipY,
3487  &(_offset[1][0]), &(_offset[1][1]),
3488  NULL
3489  ) != ES_NONE) {
3490  rterror("rt_raster_from_two_rasters: Could not compute offsets of the second raster relative to the first raster");
3491  return ES_ERROR;
3492  }
3493  _offset[1][0] = -1 * _offset[1][0];
3494  _offset[1][1] = -1 * _offset[1][1];
3495  _offset[1][2] = _offset[1][0] + _dim[1][0] - 1;
3496  _offset[1][3] = _offset[1][1] + _dim[1][1] - 1;
3497 
3498  i = -1;
3499  switch (extenttype) {
3500  case ET_FIRST:
3501  i = 0;
3502  _offset[0][0] = 0.;
3503  _offset[0][1] = 0.;
3504  case ET_LAST:
3505  case ET_SECOND:
3506  if (i < 0) {
3507  i = 1;
3508  _offset[0][0] = -1 * _offset[1][0];
3509  _offset[0][1] = -1 * _offset[1][1];
3510  _offset[1][0] = 0.;
3511  _offset[1][1] = 0.;
3512  }
3513 
3514  dim[0] = _dim[i][0];
3515  dim[1] = _dim[i][1];
3516  raster = rt_raster_new(
3517  dim[0],
3518  dim[1]
3519  );
3520  if (raster == NULL) {
3521  rterror("rt_raster_from_two_rasters: Could not create output raster");
3522  return ES_ERROR;
3523  }
3524  rt_raster_set_srid(raster, _rast[i]->srid);
3525  rt_raster_get_geotransform_matrix(_rast[i], gt);
3526  rt_raster_set_geotransform_matrix(raster, gt);
3527  break;
3528  case ET_UNION: {
3529  double off[4] = {0};
3530 
3531  rt_raster_get_geotransform_matrix(_rast[0], gt);
3532  RASTER_DEBUGF(4, "gt = (%f, %f, %f, %f, %f, %f)",
3533  gt[0],
3534  gt[1],
3535  gt[2],
3536  gt[3],
3537  gt[4],
3538  gt[5]
3539  );
3540 
3541  /* new raster upper left offset */
3542  off[0] = 0;
3543  if (_offset[1][0] < 0)
3544  off[0] = _offset[1][0];
3545  off[1] = 0;
3546  if (_offset[1][1] < 0)
3547  off[1] = _offset[1][1];
3548 
3549  /* new raster lower right offset */
3550  off[2] = _dim[0][0] - 1;
3551  if ((int) _offset[1][2] >= _dim[0][0])
3552  off[2] = _offset[1][2];
3553  off[3] = _dim[0][1] - 1;
3554  if ((int) _offset[1][3] >= _dim[0][1])
3555  off[3] = _offset[1][3];
3556 
3557  /* upper left corner */
3558  if (rt_raster_cell_to_geopoint(
3559  _rast[0],
3560  off[0], off[1],
3561  &(gt[0]), &(gt[3]),
3562  NULL
3563  ) != ES_NONE) {
3564  rterror("rt_raster_from_two_rasters: Could not get spatial coordinates of upper-left pixel of output raster");
3565  return ES_ERROR;
3566  }
3567 
3568  dim[0] = off[2] - off[0] + 1;
3569  dim[1] = off[3] - off[1] + 1;
3570  RASTER_DEBUGF(4, "off = (%f, %f, %f, %f)",
3571  off[0],
3572  off[1],
3573  off[2],
3574  off[3]
3575  );
3576  RASTER_DEBUGF(4, "dim = (%d, %d)", dim[0], dim[1]);
3577 
3578  raster = rt_raster_new(
3579  dim[0],
3580  dim[1]
3581  );
3582  if (raster == NULL) {
3583  rterror("rt_raster_from_two_rasters: Could not create output raster");
3584  return ES_ERROR;
3585  }
3586  rt_raster_set_srid(raster, _rast[0]->srid);
3587  rt_raster_set_geotransform_matrix(raster, gt);
3588  RASTER_DEBUGF(4, "gt = (%f, %f, %f, %f, %f, %f)",
3589  gt[0],
3590  gt[1],
3591  gt[2],
3592  gt[3],
3593  gt[4],
3594  gt[5]
3595  );
3596 
3597  /* get offsets */
3598  if (rt_raster_geopoint_to_cell(
3599  _rast[0],
3600  gt[0], gt[3],
3601  &(_offset[0][0]), &(_offset[0][1]),
3602  NULL
3603  ) != ES_NONE) {
3604  rterror("rt_raster_from_two_rasters: Could not get offsets of the FIRST raster relative to the output raster");
3605  rt_raster_destroy(raster);
3606  return ES_ERROR;
3607  }
3608  _offset[0][0] *= -1;
3609  _offset[0][1] *= -1;
3610 
3611  if (rt_raster_geopoint_to_cell(
3612  _rast[1],
3613  gt[0], gt[3],
3614  &(_offset[1][0]), &(_offset[1][1]),
3615  NULL
3616  ) != ES_NONE) {
3617  rterror("rt_raster_from_two_rasters: Could not get offsets of the SECOND raster relative to the output raster");
3618  rt_raster_destroy(raster);
3619  return ES_ERROR;
3620  }
3621  _offset[1][0] *= -1;
3622  _offset[1][1] *= -1;
3623  break;
3624  }
3625  case ET_INTERSECTION: {
3626  double off[4] = {0};
3627 
3628  /* no intersection */
3629  if (
3630  (_offset[1][2] < 0 || _offset[1][0] > (_dim[0][0] - 1)) ||
3631  (_offset[1][3] < 0 || _offset[1][1] > (_dim[0][1] - 1))
3632  ) {
3633  RASTER_DEBUG(3, "The two rasters provided have no intersection. Returning no band raster");
3634 
3635  raster = rt_raster_new(0, 0);
3636  if (raster == NULL) {
3637  rterror("rt_raster_from_two_rasters: Could not create output raster");
3638  return ES_ERROR;
3639  }
3640  rt_raster_set_srid(raster, _rast[0]->srid);
3641  rt_raster_set_scale(raster, 0, 0);
3642 
3643  /* set offsets if provided */
3644  if (NULL != offset) {
3645  for (i = 0; i < 4; i++)
3646  offset[i] = _offset[i / 2][i % 2];
3647  }
3648 
3649  *rtnraster = raster;
3650  return ES_NONE;
3651  }
3652 
3653  if (_offset[1][0] > 0)
3654  off[0] = _offset[1][0];
3655  if (_offset[1][1] > 0)
3656  off[1] = _offset[1][1];
3657 
3658  off[2] = _dim[0][0] - 1;
3659  if (_offset[1][2] < _dim[0][0])
3660  off[2] = _offset[1][2];
3661  off[3] = _dim[0][1] - 1;
3662  if (_offset[1][3] < _dim[0][1])
3663  off[3] = _offset[1][3];
3664 
3665  dim[0] = off[2] - off[0] + 1;
3666  dim[1] = off[3] - off[1] + 1;
3667  raster = rt_raster_new(
3668  dim[0],
3669  dim[1]
3670  );
3671  if (raster == NULL) {
3672  rterror("rt_raster_from_two_rasters: Could not create output raster");
3673  return ES_ERROR;
3674  }
3675  rt_raster_set_srid(raster, _rast[0]->srid);
3676 
3677  /* get upper-left corner */
3678  rt_raster_get_geotransform_matrix(_rast[0], gt);
3679  if (rt_raster_cell_to_geopoint(
3680  _rast[0],
3681  off[0], off[1],
3682  &(gt[0]), &(gt[3]),
3683  gt
3684  ) != ES_NONE) {
3685  rterror("rt_raster_from_two_rasters: Could not get spatial coordinates of upper-left pixel of output raster");
3686  rt_raster_destroy(raster);
3687  return ES_ERROR;
3688  }
3689 
3690  rt_raster_set_geotransform_matrix(raster, gt);
3691 
3692  /* get offsets */
3693  if (rt_raster_geopoint_to_cell(
3694  _rast[0],
3695  gt[0], gt[3],
3696  &(_offset[0][0]), &(_offset[0][1]),
3697  NULL
3698  ) != ES_NONE) {
3699  rterror("rt_raster_from_two_rasters: Could not get pixel coordinates to compute the offsets of the FIRST raster relative to the output raster");
3700  rt_raster_destroy(raster);
3701  return ES_ERROR;
3702  }
3703  _offset[0][0] *= -1;
3704  _offset[0][1] *= -1;
3705 
3706  if (rt_raster_geopoint_to_cell(
3707  _rast[1],
3708  gt[0], gt[3],
3709  &(_offset[1][0]), &(_offset[1][1]),
3710  NULL
3711  ) != ES_NONE) {
3712  rterror("rt_raster_from_two_rasters: Could not get pixel coordinates to compute the offsets of the SECOND raster relative to the output raster");
3713  rt_raster_destroy(raster);
3714  return ES_ERROR;
3715  }
3716  _offset[1][0] *= -1;
3717  _offset[1][1] *= -1;
3718  break;
3719  }
3720  case ET_CUSTOM:
3721  rterror("rt_raster_from_two_rasters: Extent type ET_CUSTOM is not supported by this function");
3722  break;
3723  }
3724 
3725  /* set offsets if provided */
3726  if (NULL != offset) {
3727  for (i = 0; i < 4; i++)
3728  offset[i] = _offset[i / 2][i % 2];
3729  }
3730 
3731  *rtnraster = raster;
3732  return ES_NONE;
3733 }
rt_raster_get_phys_params
void rt_raster_get_phys_params(rt_raster rast, double *i_mag, double *j_mag, double *theta_i, double *theta_ij)
Calculates and returns the physically significant descriptors embodied in the geotransform attached t...
Definition: rt_raster.c:231
rt_util_clamp_to_32BSI
int32_t rt_util_clamp_to_32BSI(double value)
Definition: rt_util.c:69
rt_raster_from_two_rasters
rt_errorstate rt_raster_from_two_rasters(rt_raster rast1, rt_raster rast2, rt_extenttype extenttype, rt_raster *rtnraster, double *offset)
Definition: rt_raster.c:3444
clamp_srid
int clamp_srid(int srid)
Return a valid SRID from an arbitrary integer Raises a notice if what comes out is different from wha...
Definition: lwutil.c:380
rt_raster_set_phys_params
void rt_raster_set_phys_params(rt_raster rast, double i_mag, double j_mag, double theta_i, double theta_ij)
Calculates the geotransform coefficients and applies them to the supplied raster. ...
Definition: rt_raster.c:297
rt_raster_t::skewY
double skewY
Definition: librtcore.h:2249
rt_envelope::MinY
double MinY
Definition: librtcore.h:167
rt_raster_gdal_rasterize
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)
Return a raster of the provided geometry.
Definition: rt_raster.c:2492
rt_raster_t::srid
int32_t srid
Definition: librtcore.h:2251
rt_raster_get_inverse_geotransform_matrix
rt_errorstate rt_raster_get_inverse_geotransform_matrix(rt_raster raster, double *gt, double *igt)
Get 6-element array of raster inverse geotransform matrix.
Definition: rt_raster.c:676
rt_raster_is_empty
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
Definition: rt_raster.c:1338
rt_envelope::UpperLeftX
double UpperLeftX
Definition: librtcore.h:170
rt_util_clamp_to_32BUI
uint32_t rt_util_clamp_to_32BUI(double value)
Definition: rt_util.c:74
rt_band_t
Definition: librtcore.h:2264
rt_util_gdal_register_all
int rt_util_gdal_register_all(int force_register_all)
Definition: rt_util.c:334
_rti_rasterize_arg
struct _rti_rasterize_arg_t * _rti_rasterize_arg
Definition: rt_raster.c:2399
rt_raster_get_x_offset
double rt_raster_get_x_offset(rt_raster raster)
Get raster x offset, in projection units.
Definition: rt_raster.c:213
covers
Datum covers(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:2023
rt_raster_same_alignment
rt_errorstate rt_raster_same_alignment(rt_raster rast1, rt_raster rast2, int *aligned, char **reason)
Definition: rt_spatial_relationship.c:45
rt_raster_add_band
int rt_raster_add_band(rt_raster raster, rt_band band, int index)
Add band data to a raster.
Definition: rt_raster.c:405
rtrowdump.raster
raster
Be careful!! Zeros function&#39;s input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
rt_raster_set_scale
void rt_raster_set_scale(rt_raster raster, double scaleX, double scaleY)
Set scale in projection units.
Definition: rt_raster.c:137
rt_raster_t::numBands
uint16_t numBands
Definition: librtcore.h:2241
rt_raster_replace_band
rt_band rt_raster_replace_band(rt_raster raster, rt_band band, int index)
Replace band at provided index with new band.
Definition: rt_raster.c:1502
_rti_rasterize_arg_t::numbands
uint32_t numbands
Definition: rt_raster.c:2403
rt_raster
struct rt_raster_t * rt_raster
Types definitions.
Definition: librtcore.h:145
rt_envelope::MaxY
double MaxY
Definition: librtcore.h:168
rt_band_t::raster
rt_raster raster
Definition: librtcore.h:2275
rt_envelope::MaxX
double MaxX
Definition: librtcore.h:166
lwgeom_free
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1099
rt_gdaldriver
struct rt_gdaldriver_t * rt_gdaldriver
Definition: librtcore.h:154
ovdump.band
band
Definition: ovdump.py:57
librtcore_internal.h
rt_raster_cell_to_geopoint
rt_errorstate rt_raster_cell_to_geopoint(rt_raster raster, double xr, double yr, double *xw, double *yw, double *gt)
Convert an xr, yr raster point to an xw, yw point on map.
Definition: rt_raster.c:755
rt_gdaldriver_t::short_name
char * short_name
Definition: librtcore.h:2424
rt_extenttype
rt_extenttype
Definition: librtcore.h:200
rt_raster_get_num_bands
int rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:372
rt_band_is_offline
int rt_band_is_offline(rt_band band)
Return non-zero if the given band data is on the filesystem.
Definition: rt_band.c:228
FLT_EQ
#define FLT_EQ(x, y)
Definition: librtcore.h:2185
rterror
void rterror(const char *fmt,...)
Wrappers used for reporting errors and info.
Definition: rt_context.c:199
rtalloc
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition: rt_context.c:171
ovdump.data
data
Definition: ovdump.py:103
rt_gdaldriver_t::idx
int idx
Definition: librtcore.h:2423
rtrealloc
void * rtrealloc(void *mem, size_t size)
Definition: rt_context.c:179
rt_band_t::height
uint16_t height
Definition: librtcore.h:2268
rt_util_envelope_to_lwpoly
LWPOLY * rt_util_envelope_to_lwpoly(rt_envelope ext)
Definition: rt_util.c:435
_rti_rasterize_arg_t::src_sr
OGRSpatialReferenceH src_sr
Definition: rt_raster.c:2405
rt_errorstate
rt_errorstate
Enum definitions.
Definition: librtcore.h:179
rt_raster_set_srid
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster&#39;s SRID.
Definition: rt_raster.c:363
rt_raster_t::ipY
double ipY
Definition: librtcore.h:2247
window.gt
gt
Definition: window.py:77
rt_raster_geopoint_to_cell
rt_errorstate rt_raster_geopoint_to_cell(rt_raster raster, double xw, double yw, double *xr, double *yr, double *igt)
Convert an xw,yw map point to a xr,yr raster point.
Definition: rt_raster.c:806
rt_raster_compute_skewed_raster
rt_raster rt_raster_compute_skewed_raster(rt_envelope extent, double *skew, double *scale, double tolerance)
Definition: rt_raster.c:960
_rti_rasterize_arg_t::value
double * value
Definition: rt_raster.c:2411
rt_pixtype
rt_pixtype
Definition: librtcore.h:185
rt_band_destroy
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition: rt_band.c:242
_rti_rasterize_arg_init
static _rti_rasterize_arg _rti_rasterize_arg_init()
Definition: rt_raster.c:2416
rt_raster_set_geotransform_matrix
void rt_raster_set_geotransform_matrix(rt_raster raster, double *gt)
Set raster&#39;s geotransform using 6-element array.
Definition: rt_raster.c:727
_rti_rasterize_arg_t::noband
uint8_t noband
Definition: rt_raster.c:2401
lwpoly_as_lwgeom
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:288
rt_gdaldriver_t
Definition: librtcore.h:2422
rt_raster_clone
rt_raster rt_raster_clone(rt_raster raster, uint8_t deep)
Clone an existing raster.
Definition: rt_raster.c:1544
rt_raster_get_x_scale
double rt_raster_get_x_scale(rt_raster raster)
Get scale X in projection units.
Definition: rt_raster.c:150
rt_raster_t::height
uint16_t height
Definition: librtcore.h:2253
_rti_rasterize_arg_t::init
double * init
Definition: rt_raster.c:2408
rt_raster_to_gdal_mem
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)
Return GDAL dataset using GDAL MEM driver from raster.
Definition: rt_raster.c:1809
rt_band_get_nodata
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1597
rt_raster_get_convex_hull
rt_errorstate rt_raster_get_convex_hull(rt_raster raster, LWGEOM **hull)
Get raster&#39;s convex hull.
Definition: rt_geometry.c:803
pixval.ds
ds
Definition: pixval.py:66
uint32_t
unsigned int uint32_t
Definition: uthash.h:78
LW_PARSER_CHECK_NONE
#define LW_PARSER_CHECK_NONE
Definition: liblwgeom.h:2013
rt_raster_get_y_scale
double rt_raster_get_y_scale(rt_raster raster)
Get scale Y in projection units.
Definition: rt_raster.c:159
rt_raster_from_gdal_dataset
rt_raster rt_raster_from_gdal_dataset(GDALDatasetH ds)
Return a raster from a GDAL dataset.
Definition: rt_raster.c:2165
rt_util_clamp_to_16BUI
uint16_t rt_util_clamp_to_16BUI(double value)
Definition: rt_util.c:64
_rti_rasterize_arg_destroy
static void _rti_rasterize_arg_destroy(_rti_rasterize_arg arg)
Definition: rt_raster.c:2442
_rti_rasterize_arg_t::nodata
double * nodata
Definition: rt_raster.c:2409
rt_band_set_ownsdata_flag
void rt_band_set_ownsdata_flag(rt_band band, int flag)
Definition: rt_band.c:534
pixval.nband
nband
Definition: pixval.py:52
rt_util_pixtype_to_gdal_datatype
GDALDataType rt_util_pixtype_to_gdal_datatype(rt_pixtype pt)
Convert rt_pixtype to GDALDataType.
Definition: rt_util.c:118
rt_util_gdal_driver_registered
int rt_util_gdal_driver_registered(const char *drv)
Definition: rt_util.c:353
lwgeom_geos_error
void lwgeom_geos_error(const char *fmt,...)
Definition: liblwgeom/lwgeom_geos.c:42
rtwarn
void rtwarn(const char *fmt,...)
Definition: rt_context.c:224
lwpoly_free
void lwpoly_free(LWPOLY *poly)
Definition: lwpoly.c:174
rt_envelope
Definition: librtcore.h:164
rt_band_get_pixel
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_band.c:1088
rt_raster_copy_band
int rt_raster_copy_band(rt_raster torast, rt_raster fromrast, int fromindex, int toindex)
Copy one band from one raster to another.
Definition: rt_raster.c:1374
rt_raster_t::width
uint16_t width
Definition: librtcore.h:2252
SRID_UNKNOWN
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:188
rt_raster_get_geotransform_matrix
void rt_raster_get_geotransform_matrix(rt_raster raster, double *gt)
Get 6-element array of raster geotransform matrix.
Definition: rt_raster.c:706
rt_util_clamp_to_8BSI
int8_t rt_util_clamp_to_8BSI(double value)
Definition: rt_util.c:49
rt_util_dbl_trunc_warning
int rt_util_dbl_trunc_warning(double initialvalue, int32_t checkvalint, uint32_t checkvaluint, float checkvalfloat, double checkvaldouble, rt_pixtype pixtype)
Definition: rt_util.c:625
rt_raster_get_y_offset
double rt_raster_get_y_offset(rt_raster raster)
Get raster y offset, in projection units.
Definition: rt_raster.c:222
FLT_NEQ
#define FLT_NEQ(x, y)
Definition: librtcore.h:2184
rt_util_gdal_datatype_to_pixtype
rt_pixtype rt_util_gdal_datatype_to_pixtype(GDALDataType gdt)
Convert GDALDataType to rt_pixtype.
Definition: rt_util.c:153
rtinfo
void rtinfo(const char *fmt,...)
Definition: rt_context.c:211
rt_band_get_hasnodata_flag
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:541
genraster.count
int count
Definition: genraster.py:56
rt_raster_t::skewX
double skewX
Definition: librtcore.h:2248
rt_band_get_data
void * rt_band_get_data(rt_band band)
Get pointer to raster band data.
Definition: rt_band.c:302
_rti_rasterize_arg_t::pixtype
rt_pixtype * pixtype
Definition: rt_raster.c:2407
RASTER_DEBUGF
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:299
rt_raster_get_x_skew
double rt_raster_get_x_skew(rt_raster raster)
Get skew about the X axis.
Definition: rt_raster.c:181
rt_raster_gdal_drivers
rt_gdaldriver rt_raster_gdal_drivers(uint32_t *drv_count, uint8_t cancc)
Returns a set of available GDAL drivers.
Definition: rt_raster.c:1705
LWGEOM2GEOS
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, int autofix)
Definition: liblwgeom/lwgeom_geos.c:357
rt_raster_set_offsets
void rt_raster_set_offsets(rt_raster raster, double x, double y)
Set insertion points in projection units.
Definition: rt_raster.c:199
rt_band_t::width
uint16_t width
Definition: librtcore.h:2267
_rt_raster_geotransform_warn_offline_band
static void _rt_raster_geotransform_warn_offline_band(rt_raster raster)
Definition: rt_raster.c:95
rt_band_new_inline
rt_band rt_band_new_inline(uint16_t width, uint16_t height, rt_pixtype pixtype, uint32_t hasnodata, double nodataval, uint8_t *data)
Create an in-db rt_band with no data.
Definition: rt_band.c:58
rt_raster_get_srid
int32_t rt_raster_get_srid(rt_raster raster)
Get raster&#39;s SRID.
Definition: rt_raster.c:356
rt_util_clamp_to_4BUI
uint8_t rt_util_clamp_to_4BUI(double value)
Definition: rt_util.c:44
rt_util_clamp_to_8BUI
uint8_t rt_util_clamp_to_8BUI(double value)
Definition: rt_util.c:54
rt_pixtype_size
int rt_pixtype_size(rt_pixtype pixtype)
Return size in bytes of a value in the given pixtype.
Definition: rt_pixel.c:39
rt_util_clamp_to_1BB
uint8_t rt_util_clamp_to_1BB(double value)
Definition: rt_util.c:34
rt_raster_new
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_raster.c:48
_rti_rasterize_arg_t::bandlist
int * bandlist
Definition: rt_raster.c:2412
rt_band_set_isnodata_flag
rt_errorstate rt_band_set_isnodata_flag(rt_band band, int flag)
Set isnodata flag value.
Definition: rt_band.c:562
rt_raster_get_envelope
rt_errorstate rt_raster_get_envelope(rt_raster raster, rt_envelope *env)
Get raster&#39;s envelope.
Definition: rt_raster.c:873
rtdealloc
void rtdealloc(void *mem)
Definition: rt_context.c:186
_rti_rasterize_arg_t
Definition: rt_raster.c:2400
rt_raster_calc_phys_params
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)
Calculates the physically significant descriptors embodied in an arbitrary geotransform.
Definition: rt_raster.c:250
rt_raster_set_skews
void rt_raster_set_skews(rt_raster raster, double skewX, double skewY)
Set skews about the X and Y axis.
Definition: rt_raster.c:168
ROUND
#define ROUND(x, y)
Definition: librtcore.h:2192
FALSE
#define FALSE
Definition: dbfopen.c:168
rt_envelope::MinX
double MinX
Definition: librtcore.h:165
RASTER_DEBUG
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:295
rt_band_duplicate
rt_band rt_band_duplicate(rt_band band)
Create a new band duplicated from source band.
Definition: rt_band.c:186
librtcore.h
This library is the generic raster handling section of PostGIS.
rt_util_gdal_convert_sr
char * rt_util_gdal_convert_sr(const char *srs, int proj4)
Definition: rt_util.c:212
rt_raster_t::bands
rt_band * bands
Definition: librtcore.h:2254
rt_gdaldriver_t::create_options
char * create_options
Definition: librtcore.h:2426
rt_band_get_pixtype
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition: rt_band.c:498
rt_band_set_pixel_line
rt_errorstate rt_band_set_pixel_line(rt_band band, int x, int y, void *vals, uint32_t len)
Set values of multiple pixels.
Definition: rt_band.c:720
rt_raster_has_band
int rt_raster_has_band(rt_raster raster, int nband)
Return TRUE if the raster has a band of this number.
Definition: rt_raster.c:1351
rt_raster_get_y_skew
double rt_raster_get_y_skew(rt_raster raster)
Get skew about the Y axis.
Definition: rt_raster.c:190
rt_gdaldriver_t::long_name
char * long_name
Definition: librtcore.h:2425
rt_raster_generate_new_band
int rt_raster_generate_new_band(rt_raster raster, rt_pixtype pixtype, double initialvalue, uint32_t hasnodata, double nodatavalue, int index)
Generate a new inline band and add it to a raster.
Definition: rt_raster.c:485
rt_util_gdal_sr_auth_info
rt_errorstate rt_util_gdal_sr_auth_info(GDALDatasetH hds, char **authname, char **authcode)
Get auth name and code.
Definition: rt_util.c:270
rt_util_from_ogr_envelope
void rt_util_from_ogr_envelope(OGREnvelope env, rt_envelope *ext)
Definition: rt_util.c:406
ovdump.src_ds
src_ds
MAIN.
Definition: ovdump.py:53
rt_raster_to_gdal
uint8_t * rt_raster_to_gdal(rt_raster raster, const char *srs, char *format, char **options, uint64_t *gdalsize)
Return formatted GDAL raster from raster.
Definition: rt_raster.c:1602
rt_raster_t::ipX
double ipX
Definition: librtcore.h:2246
rt_util_clamp_to_2BUI
uint8_t rt_util_clamp_to_2BUI(double value)
Definition: rt_util.c:39
rt_raster_t
Definition: librtcore.h:2235
rt_util_clamp_to_32F
float rt_util_clamp_to_32F(double value)
Definition: rt_util.c:79
rt_band_set_pixel
rt_errorstate rt_band_set_pixel(rt_band band, int x, int y, double val, int *converted)
Set single pixel&#39;s value.
Definition: rt_band.c:841
genraster.value
int value
Definition: genraster.py:61
rt_raster_from_band
rt_raster rt_raster_from_band(rt_raster raster, uint32_t *bandNums, int count)
Construct a new rt_raster from an existing rt_raster and an array of band numbers.
Definition: rt_raster.c:1439
rt_raster_calc_gt_coeff
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)
Calculates the coefficients of a geotransform given the physically significant parameters describing ...
Definition: rt_raster.c:314
rt_raster_t::scaleX
double scaleX
Definition: librtcore.h:2244
lwgeom_from_wkb
LWGEOM * lwgeom_from_wkb(const uint8_t *wkb, const size_t wkb_size, const char check)
WKB inputs must have a declared size, to prevent malformed WKB from reading off the end of the memory...
Definition: lwin_wkb.c:772
uint8_t
unsigned char uint8_t
Definition: uthash.h:79
rt_envelope::UpperLeftY
double UpperLeftY
Definition: librtcore.h:171
rt_raster_get_height
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_raster.c:129
rt_raster_get_width
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_raster.c:121
TRUE
#define TRUE
Definition: dbfopen.c:169
return
return(const char *)
Definition: dbfopen.c:1613
rt_raster_t::scaleY
double scaleY
Definition: librtcore.h:2245
rt_raster_destroy
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:82
rt_raster_get_band
rt_band rt_raster_get_band(rt_raster raster, int n)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:381
_rti_rasterize_arg_t::hasnodata
uint8_t * hasnodata
Definition: rt_raster.c:2410
rt_util_clamp_to_16BSI
int16_t rt_util_clamp_to_16BSI(double value)
Definition: rt_util.c:59