PostGIS  2.5.7dev-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 
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
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 
92 }
93 
94 static void
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++) {
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
122 
123  assert(NULL != raster);
124 
125  return raster->width;
126 }
127 
128 uint16_t
130 
131  assert(NULL != raster);
132 
133  return raster->height;
134 }
135 
136 void
139  double scaleX, double scaleY
140 ) {
141  assert(NULL != raster);
142 
143  raster->scaleX = scaleX;
144  raster->scaleY = scaleY;
145 
147 }
148 
149 double
151 
152 
153  assert(NULL != raster);
154 
155  return raster->scaleX;
156 }
157 
158 double
160 
161 
162  assert(NULL != raster);
163 
164  return raster->scaleY;
165 }
166 
167 void
170  double skewX, double skewY
171 ) {
172  assert(NULL != raster);
173 
174  raster->skewX = skewX;
175  raster->skewY = skewY;
176 
178 }
179 
180 double
182 
183 
184  assert(NULL != raster);
185 
186  return raster->skewX;
187 }
188 
189 double
191 
192 
193  assert(NULL != raster);
194 
195  return raster->skewY;
196 }
197 
198 void
201  double x, double y
202 ) {
203 
204  assert(NULL != raster);
205 
206  raster->ipX = x;
207  raster->ipY = y;
208 
210 }
211 
212 double
214 
215 
216  assert(NULL != raster);
217 
218  return raster->ipX;
219 }
220 
221 double
223 
224 
225  assert(NULL != raster);
226 
227  return raster->ipY;
228 }
229 
230 void
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
357  assert(NULL != raster);
358 
359  return clamp_srid(raster->srid);
360 }
361 
362 void
364  assert(NULL != raster);
365 
366  raster->srid = clamp_srid(srid);
367 
369 }
370 
371 uint16_t
373 
374 
375  assert(NULL != raster);
376 
377  return raster->numBands;
378 }
379 
380 rt_band
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
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
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 */
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");
658  }
659 
660  /* set isnodata if hasnodata = TRUE and initial value = nodatavalue */
661  if (hasnodata && FLT_EQ(initialvalue, nodatavalue))
663 
664  return index;
665 }
666 
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)
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
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
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 
740 }
741 
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  ) {
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 
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  ) {
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 
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 
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 
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
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;
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  }
1070 
1071  /* get inverse geotransform matrix */
1072  if (!GDALInvGeoTransform(_gt, _igt)) {
1073  rterror("rt_raster_compute_skewed_raster: Could not compute inverse geotransform matrix");
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");
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 
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");
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 
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");
1178  return NULL;
1179  }
1180 
1181  /* adjust ul */
1182  if (i < 1)
1183  _gt[0] = _w[i];
1184  else
1185  _gt[3] = _w[i];
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");
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 */
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");
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");
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);
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);
1262  return NULL;
1263  }
1264 
1265  if (!covers)
1266  {
1267  raster->width++;
1268  raster->height++;
1269  }
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];
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);
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);
1310  return NULL;
1311  }
1312  } while (covers);
1313 
1314  if (i < 1)
1315  raster->width++;
1316  else
1317  raster->height++;
1318 
1319  }
1320 
1321  GEOSGeom_destroy(ngeom);
1322 
1323  return raster;
1324 }
1325 
1333 int
1335  return (!raster || raster->height == 0 || raster->width == 0);
1336 }
1337 
1346 int
1348  return !(NULL == raster || nband >= raster->numBands || nband < 0);
1349 }
1350 
1351 /******************************************************************************
1352 * rt_raster_copy_band()
1353 ******************************************************************************/
1354 
1369 int
1371  rt_raster torast, rt_raster fromrast,
1372  int fromindex, int toindex
1373 ) {
1374  rt_band srcband = NULL;
1375  rt_band dstband = NULL;
1376 
1377  assert(NULL != torast);
1378  assert(NULL != fromrast);
1379 
1380  /* Check raster dimensions */
1381  if (torast->height != fromrast->height || torast->width != fromrast->width) {
1382  rtwarn("rt_raster_copy_band: Attempting to add a band with different width or height");
1383  return -1;
1384  }
1385 
1386  /* Check bands limits */
1387  if (fromrast->numBands < 1) {
1388  rtwarn("rt_raster_copy_band: Second raster has no band");
1389  return -1;
1390  }
1391  else if (fromindex < 0) {
1392  rtwarn("rt_raster_copy_band: Band index for second raster < 0. Defaulted to 0");
1393  fromindex = 0;
1394  }
1395  else if (fromindex >= fromrast->numBands) {
1396  rtwarn("rt_raster_copy_band: Band index for second raster > number of bands, truncated from %u to %u", fromindex, fromrast->numBands - 1);
1397  fromindex = fromrast->numBands - 1;
1398  }
1399 
1400  if (toindex < 0) {
1401  rtwarn("rt_raster_copy_band: Band index for first raster < 0. Defaulted to 0");
1402  toindex = 0;
1403  }
1404  else if (toindex > torast->numBands) {
1405  rtwarn("rt_raster_copy_band: Band index for first raster > number of bands, truncated from %u to %u", toindex, torast->numBands);
1406  toindex = torast->numBands;
1407  }
1408 
1409  /* Get band from source raster */
1410  srcband = rt_raster_get_band(fromrast, fromindex);
1411 
1412  /* duplicate band */
1413  dstband = rt_band_duplicate(srcband);
1414 
1415  /* Add band to the second raster */
1416  return rt_raster_add_band(torast, dstband, toindex);
1417 }
1418 
1419 /******************************************************************************
1420 * rt_raster_from_band()
1421 ******************************************************************************/
1422 
1434 rt_raster
1436  rt_raster rast = NULL;
1437  int i = 0;
1438  int j = 0;
1439  int idx;
1440  int32_t flag;
1441  double gt[6] = {0.};
1442 
1443  assert(NULL != raster);
1444  assert(NULL != bandNums);
1445 
1446  RASTER_DEBUGF(3, "rt_raster_from_band: source raster has %d bands",
1448 
1449  /* create new raster */
1450  rast = rt_raster_new(raster->width, raster->height);
1451  if (NULL == rast) {
1452  rterror("rt_raster_from_band: Out of memory allocating new raster");
1453  return NULL;
1454  }
1455 
1456  /* copy raster attributes */
1459 
1460  /* srid */
1461  rt_raster_set_srid(rast, raster->srid);
1462 
1463  /* copy bands */
1464  for (i = 0; i < count; i++) {
1465  idx = bandNums[i];
1466  flag = rt_raster_copy_band(rast, raster, idx, i);
1467 
1468  if (flag < 0) {
1469  rterror("rt_raster_from_band: Could not copy band");
1470  for (j = 0; j < i; j++) rt_band_destroy(rast->bands[j]);
1472  return NULL;
1473  }
1474 
1475  RASTER_DEBUGF(3, "rt_raster_from_band: band created at index %d",
1476  flag);
1477  }
1478 
1479  RASTER_DEBUGF(3, "rt_raster_from_band: new raster has %d bands",
1481  return rast;
1482 }
1483 
1484 /******************************************************************************
1485 * rt_raster_replace_band()
1486 ******************************************************************************/
1487 
1497 rt_band
1499  rt_band oldband = NULL;
1500  assert(NULL != raster);
1501  assert(NULL != band);
1502 
1503  if (band->width != raster->width || band->height != raster->height) {
1504  rterror("rt_raster_replace_band: Band does not match raster's dimensions: %dx%d band to %dx%d raster",
1505  band->width, band->height, raster->width, raster->height);
1506  return 0;
1507  }
1508 
1509  if (index >= raster->numBands || index < 0) {
1510  rterror("rt_raster_replace_band: Band index is not valid");
1511  return 0;
1512  }
1513 
1514  oldband = rt_raster_get_band(raster, index);
1515  RASTER_DEBUGF(3, "rt_raster_replace_band: old band at %p", oldband);
1516  RASTER_DEBUGF(3, "rt_raster_replace_band: new band at %p", band);
1517 
1518  raster->bands[index] = band;
1519  RASTER_DEBUGF(3, "rt_raster_replace_band: new band at %p", raster->bands[index]);
1520 
1521  band->raster = raster;
1522  oldband->raster = NULL;
1523 
1524  return oldband;
1525 }
1526 
1527 /******************************************************************************
1528 * rt_raster_clone()
1529 ******************************************************************************/
1530 
1539 rt_raster
1541  rt_raster rtn = NULL;
1542  double gt[6] = {0};
1543 
1544  assert(NULL != raster);
1545 
1546  if (deep) {
1547  int numband = rt_raster_get_num_bands(raster);
1548  uint32_t *nband = NULL;
1549  int i = 0;
1550 
1551  nband = rtalloc(sizeof(uint32_t) * numband);
1552  if (nband == NULL) {
1553  rterror("rt_raster_clone: Could not allocate memory for deep clone");
1554  return NULL;
1555  }
1556  for (i = 0; i < numband; i++)
1557  nband[i] = i;
1558 
1559  rtn = rt_raster_from_band(raster, nband, numband);
1560  rtdealloc(nband);
1561 
1562  return rtn;
1563  }
1564 
1565  rtn = rt_raster_new(
1568  );
1569  if (rtn == NULL) {
1570  rterror("rt_raster_clone: Could not create cloned raster");
1571  return NULL;
1572  }
1573 
1577 
1578  return rtn;
1579 }
1580 
1581 /******************************************************************************
1582 * rt_raster_to_gdal()
1583 ******************************************************************************/
1584 
1597 uint8_t*
1599  rt_raster raster, const char *srs,
1600  char *format, char **options, uint64_t *gdalsize
1601 ) {
1602  const char *cc;
1603  const char *vio;
1604 
1605  GDALDriverH src_drv = NULL;
1606  int destroy_src_drv = 0;
1607  GDALDatasetH src_ds = NULL;
1608 
1609  vsi_l_offset rtn_lenvsi;
1610  uint64_t rtn_len = 0;
1611 
1612  GDALDriverH rtn_drv = NULL;
1613  GDALDatasetH rtn_ds = NULL;
1614  uint8_t *rtn = NULL;
1615 
1616  assert(NULL != raster);
1617  assert(NULL != gdalsize);
1618 
1619  /* any supported format is possible */
1621  RASTER_DEBUG(3, "loaded all supported GDAL formats");
1622 
1623  /* output format not specified */
1624  if (format == NULL || !strlen(format))
1625  format = "GTiff";
1626  RASTER_DEBUGF(3, "output format is %s", format);
1627 
1628  /* load raster into a GDAL MEM raster */
1629  src_ds = rt_raster_to_gdal_mem(raster, srs, NULL, NULL, 0, &src_drv, &destroy_src_drv);
1630  if (NULL == src_ds) {
1631  rterror("rt_raster_to_gdal: Could not convert raster to GDAL MEM format");
1632  return 0;
1633  }
1634 
1635  /* load driver */
1636  rtn_drv = GDALGetDriverByName(format);
1637  if (NULL == rtn_drv) {
1638  rterror("rt_raster_to_gdal: Could not load the output GDAL driver");
1639  GDALClose(src_ds);
1640  if (destroy_src_drv) GDALDestroyDriver(src_drv);
1641  return 0;
1642  }
1643  RASTER_DEBUG(3, "Output driver loaded");
1644 
1645  /* CreateCopy support */
1646  cc = GDALGetMetadataItem(rtn_drv, GDAL_DCAP_CREATECOPY, NULL);
1647  /* VirtualIO support */
1648  vio = GDALGetMetadataItem(rtn_drv, GDAL_DCAP_VIRTUALIO, NULL);
1649 
1650  if (cc == NULL || vio == NULL) {
1651  rterror("rt_raster_to_gdal: Output GDAL driver does not support CreateCopy and/or VirtualIO");
1652  GDALClose(src_ds);
1653  if (destroy_src_drv) GDALDestroyDriver(src_drv);
1654  return 0;
1655  }
1656 
1657  /* convert GDAL MEM raster to output format */
1658  RASTER_DEBUG(3, "Copying GDAL MEM raster to memory file in output format");
1659  rtn_ds = GDALCreateCopy(
1660  rtn_drv,
1661  "/vsimem/out.dat", /* should be fine assuming this is in a process */
1662  src_ds,
1663  FALSE, /* should copy be strictly equivelent? */
1664  options, /* format options */
1665  NULL, /* progress function */
1666  NULL /* progress data */
1667  );
1668 
1669  /* close source dataset */
1670  GDALClose(src_ds);
1671  if (destroy_src_drv) GDALDestroyDriver(src_drv);
1672  RASTER_DEBUG(3, "Closed GDAL MEM raster");
1673 
1674  if (NULL == rtn_ds) {
1675  rterror("rt_raster_to_gdal: Could not create the output GDAL dataset");
1676  return 0;
1677  }
1678 
1679  RASTER_DEBUGF(4, "dataset SRS: %s", GDALGetProjectionRef(rtn_ds));
1680 
1681  /* close dataset, this also flushes any pending writes */
1682  GDALClose(rtn_ds);
1683  RASTER_DEBUG(3, "Closed GDAL output raster");
1684 
1685  RASTER_DEBUG(3, "Done copying GDAL MEM raster to memory file in output format");
1686 
1687  /* from memory file to buffer */
1688  RASTER_DEBUG(3, "Copying GDAL memory file to buffer");
1689  rtn = VSIGetMemFileBuffer("/vsimem/out.dat", &rtn_lenvsi, TRUE);
1690  RASTER_DEBUG(3, "Done copying GDAL memory file to buffer");
1691  if (NULL == rtn) {
1692  rterror("rt_raster_to_gdal: Could not create the output GDAL raster");
1693  return 0;
1694  }
1695 
1696  rtn_len = (uint64_t) rtn_lenvsi;
1697  *gdalsize = rtn_len;
1698 
1699  return rtn;
1700 }
1701 
1702 /******************************************************************************
1703 * rt_raster_gdal_drivers()
1704 ******************************************************************************/
1705 
1716 rt_raster_gdal_drivers(uint32_t *drv_count, uint8_t can_write) {
1717  const char *cc;
1718  const char *vio;
1719  const char *txt;
1720  int txt_len;
1721  GDALDriverH *drv = NULL;
1722  rt_gdaldriver rtn = NULL;
1723  int count;
1724  int i;
1725  uint32_t j;
1726 
1727  assert(drv_count != NULL);
1728 
1730  count = GDALGetDriverCount();
1731  RASTER_DEBUGF(3, "%d drivers found", count);
1732 
1733  rtn = (rt_gdaldriver) rtalloc(count * sizeof(struct rt_gdaldriver_t));
1734  if (NULL == rtn) {
1735  rterror("rt_raster_gdal_drivers: Could not allocate memory for gdaldriver structure");
1736  return 0;
1737  }
1738 
1739  for (i = 0, j = 0; i < count; i++) {
1740  drv = GDALGetDriver(i);
1741 
1742 #ifdef GDAL_DCAP_RASTER
1743  /* Starting with GDAL 2.0, vector drivers can also be returned */
1744  /* Only keep raster drivers */
1745  const char *is_raster;
1746  is_raster = GDALGetMetadataItem(drv, GDAL_DCAP_RASTER, NULL);
1747  if (is_raster == NULL || !EQUAL(is_raster, "YES"))
1748  continue;
1749 #endif
1750 
1751  /* CreateCopy support */
1752  cc = GDALGetMetadataItem(drv, GDAL_DCAP_CREATECOPY, NULL);
1753 
1754  /* VirtualIO support */
1755  vio = GDALGetMetadataItem(drv, GDAL_DCAP_VIRTUALIO, NULL);
1756 
1757  if (can_write && (cc == NULL || vio == NULL))
1758  continue;
1759 
1760  /* we can always read what GDAL can load */
1761  rtn[j].can_read = 1;
1762  /* we require CreateCopy and VirtualIO support to write to GDAL */
1763  rtn[j].can_write = (cc != NULL && vio != NULL);
1764 
1765  if (rtn[j].can_write) {
1766  RASTER_DEBUGF(3, "driver %s (%d) supports CreateCopy() and VirtualIO()", txt, i);
1767  }
1768 
1769  /* index of driver */
1770  rtn[j].idx = i;
1771 
1772  /* short name */
1773  txt = GDALGetDriverShortName(drv);
1774  txt_len = strlen(txt);
1775 
1776  txt_len = (txt_len + 1) * sizeof(char);
1777  rtn[j].short_name = (char *) rtalloc(txt_len);
1778  memcpy(rtn[j].short_name, txt, txt_len);
1779 
1780  /* long name */
1781  txt = GDALGetDriverLongName(drv);
1782  txt_len = strlen(txt);
1783 
1784  txt_len = (txt_len + 1) * sizeof(char);
1785  rtn[j].long_name = (char *) rtalloc(txt_len);
1786  memcpy(rtn[j].long_name, txt, txt_len);
1787 
1788  /* creation options */
1789  txt = GDALGetDriverCreationOptionList(drv);
1790  txt_len = strlen(txt);
1791 
1792  txt_len = (txt_len + 1) * sizeof(char);
1793  rtn[j].create_options = (char *) rtalloc(txt_len);
1794  memcpy(rtn[j].create_options, txt, txt_len);
1795 
1796  j++;
1797  }
1798 
1799  /* free unused memory */
1800  rtn = rtrealloc(rtn, j * sizeof(struct rt_gdaldriver_t));
1801  *drv_count = j;
1802 
1803  return rtn;
1804 }
1805 
1806 /******************************************************************************
1807 * rt_raster_to_gdal_mem()
1808 ******************************************************************************/
1809 
1825 GDALDatasetH
1827  rt_raster raster,
1828  const char *srs,
1829  uint32_t *bandNums,
1830  int *excludeNodataValues,
1831  int count,
1832  GDALDriverH *rtn_drv, int *destroy_rtn_drv
1833 ) {
1834  GDALDriverH drv = NULL;
1835  GDALDatasetH ds = NULL;
1836  double gt[6] = {0.0};
1837  CPLErr cplerr;
1838  GDALDataType gdal_pt = GDT_Unknown;
1839  GDALRasterBandH band;
1840  void *pVoid;
1841  char *pszDataPointer;
1842  char szGDALOption[50];
1843  char *apszOptions[4];
1844  double nodata = 0.0;
1845  int allocBandNums = 0;
1846  int allocNodataValues = 0;
1847 
1848  int i;
1849  uint32_t numBands;
1850  uint32_t width = 0;
1851  uint32_t height = 0;
1852  rt_band rtband = NULL;
1853  rt_pixtype pt = PT_END;
1854 
1855  assert(NULL != raster);
1856  assert(NULL != rtn_drv);
1857  assert(NULL != destroy_rtn_drv);
1858 
1859  *destroy_rtn_drv = 0;
1860 
1861  /* store raster in GDAL MEM raster */
1862  if (!rt_util_gdal_driver_registered("MEM")) {
1863  RASTER_DEBUG(4, "Registering MEM driver");
1864  GDALRegister_MEM();
1865  *destroy_rtn_drv = 1;
1866  }
1867  drv = GDALGetDriverByName("MEM");
1868  if (NULL == drv) {
1869  rterror("rt_raster_to_gdal_mem: Could not load the MEM GDAL driver");
1870  return 0;
1871  }
1872  *rtn_drv = drv;
1873 
1874  /* unload driver from GDAL driver manager */
1875  if (*destroy_rtn_drv) {
1876  RASTER_DEBUG(4, "Deregistering MEM driver");
1877  GDALDeregisterDriver(drv);
1878  }
1879 
1880  width = rt_raster_get_width(raster);
1881  height = rt_raster_get_height(raster);
1882  ds = GDALCreate(
1883  drv, "",
1884  width, height,
1885  0, GDT_Byte, NULL
1886  );
1887  if (NULL == ds) {
1888  rterror("rt_raster_to_gdal_mem: Could not create a GDALDataset to convert into");
1889  return 0;
1890  }
1891 
1892  /* add geotransform */
1894  cplerr = GDALSetGeoTransform(ds, gt);
1895  if (cplerr != CE_None) {
1896  rterror("rt_raster_to_gdal_mem: Could not set geotransformation");
1897  GDALClose(ds);
1898  return 0;
1899  }
1900 
1901  /* set spatial reference */
1902  if (NULL != srs && strlen(srs)) {
1903  char *_srs = rt_util_gdal_convert_sr(srs, 0);
1904  if (_srs == NULL) {
1905  rterror("rt_raster_to_gdal_mem: Could not convert srs to GDAL accepted format");
1906  GDALClose(ds);
1907  return 0;
1908  }
1909 
1910  cplerr = GDALSetProjection(ds, _srs);
1911  CPLFree(_srs);
1912  if (cplerr != CE_None) {
1913  rterror("rt_raster_to_gdal_mem: Could not set projection");
1914  GDALClose(ds);
1915  return 0;
1916  }
1917  RASTER_DEBUGF(3, "Projection set to: %s", GDALGetProjectionRef(ds));
1918  }
1919 
1920  /* process bandNums */
1921  numBands = rt_raster_get_num_bands(raster);
1922  if (NULL != bandNums && count > 0) {
1923  for (i = 0; i < count; i++) {
1924  if (bandNums[i] >= numBands) {
1925  rterror("rt_raster_to_gdal_mem: The band index %d is invalid", bandNums[i]);
1926  GDALClose(ds);
1927  return 0;
1928  }
1929  }
1930  }
1931  else {
1932  count = numBands;
1933  bandNums = (uint32_t *) rtalloc(sizeof(uint32_t) * count);
1934  if (NULL == bandNums) {
1935  rterror("rt_raster_to_gdal_mem: Could not allocate memory for band indices");
1936  GDALClose(ds);
1937  return 0;
1938  }
1939  allocBandNums = 1;
1940  for (i = 0; i < count; i++) bandNums[i] = i;
1941  }
1942 
1943  /* process exclude_nodata_values */
1944  if (NULL == excludeNodataValues) {
1945  excludeNodataValues = (int *) rtalloc(sizeof(int) * count);
1946  if (NULL == excludeNodataValues) {
1947  rterror("rt_raster_to_gdal_mem: Could not allocate memory for NODATA flags");
1948  GDALClose(ds);
1949  return 0;
1950  }
1951  allocNodataValues = 1;
1952  for (i = 0; i < count; i++) excludeNodataValues[i] = 1;
1953  }
1954 
1955  /* add band(s) */
1956  for (i = 0; i < count; i++) {
1957  rtband = rt_raster_get_band(raster, bandNums[i]);
1958  if (NULL == rtband) {
1959  rterror("rt_raster_to_gdal_mem: Could not get requested band index %d", bandNums[i]);
1960  if (allocBandNums) rtdealloc(bandNums);
1961  if (allocNodataValues) rtdealloc(excludeNodataValues);
1962  GDALClose(ds);
1963  return 0;
1964  }
1965 
1966  pt = rt_band_get_pixtype(rtband);
1967  gdal_pt = rt_util_pixtype_to_gdal_datatype(pt);
1968  if (gdal_pt == GDT_Unknown)
1969  rtwarn("rt_raster_to_gdal_mem: Unknown pixel type for band");
1970 
1971  /*
1972  For all pixel types other than PT_8BSI, set pointer to start of data
1973  */
1974  if (pt != PT_8BSI) {
1975  pVoid = rt_band_get_data(rtband);
1976  RASTER_DEBUGF(4, "Band data is at pos %p", pVoid);
1977 
1978  pszDataPointer = (char *) rtalloc(20 * sizeof (char));
1979  sprintf(pszDataPointer, "%p", pVoid);
1980  RASTER_DEBUGF(4, "rt_raster_to_gdal_mem: szDatapointer is %p",
1981  pszDataPointer);
1982 
1983  if (strncasecmp(pszDataPointer, "0x", 2) == 0)
1984  sprintf(szGDALOption, "DATAPOINTER=%s", pszDataPointer);
1985  else
1986  sprintf(szGDALOption, "DATAPOINTER=0x%s", pszDataPointer);
1987 
1988  RASTER_DEBUG(3, "Storing info for GDAL MEM raster band");
1989 
1990  apszOptions[0] = szGDALOption;
1991  apszOptions[1] = NULL; /* pixel offset, not needed */
1992  apszOptions[2] = NULL; /* line offset, not needed */
1993  apszOptions[3] = NULL;
1994 
1995  /* free */
1996  rtdealloc(pszDataPointer);
1997 
1998  /* add band */
1999  if (GDALAddBand(ds, gdal_pt, apszOptions) == CE_Failure) {
2000  rterror("rt_raster_to_gdal_mem: Could not add GDAL raster band");
2001  if (allocBandNums) rtdealloc(bandNums);
2002  GDALClose(ds);
2003  return 0;
2004  }
2005  }
2006  /*
2007  PT_8BSI is special as GDAL has no equivalent pixel type.
2008  Must convert 8BSI to 16BSI so create basic band
2009  */
2010  else {
2011  /* add band */
2012  if (GDALAddBand(ds, gdal_pt, NULL) == CE_Failure) {
2013  rterror("rt_raster_to_gdal_mem: Could not add GDAL raster band");
2014  if (allocBandNums) rtdealloc(bandNums);
2015  if (allocNodataValues) rtdealloc(excludeNodataValues);
2016  GDALClose(ds);
2017  return 0;
2018  }
2019  }
2020 
2021  /* check band count */
2022  if (GDALGetRasterCount(ds) != i + 1) {
2023  rterror("rt_raster_to_gdal_mem: Error creating GDAL MEM raster band");
2024  if (allocBandNums) rtdealloc(bandNums);
2025  if (allocNodataValues) rtdealloc(excludeNodataValues);
2026  GDALClose(ds);
2027  return 0;
2028  }
2029 
2030  /* get new band */
2031  band = NULL;
2032  band = GDALGetRasterBand(ds, i + 1);
2033  if (NULL == band) {
2034  rterror("rt_raster_to_gdal_mem: Could not get GDAL band for additional processing");
2035  if (allocBandNums) rtdealloc(bandNums);
2036  if (allocNodataValues) rtdealloc(excludeNodataValues);
2037  GDALClose(ds);
2038  return 0;
2039  }
2040 
2041  /* PT_8BSI requires manual setting of pixels */
2042  if (pt == PT_8BSI) {
2043  uint32_t nXBlocks, nYBlocks;
2044  int nXBlockSize, nYBlockSize;
2045  uint32_t iXBlock, iYBlock;
2046  int nXValid, nYValid;
2047  int iX, iY;
2048  int iXMax, iYMax;
2049 
2050  int x, y, z;
2051  uint32_t valueslen = 0;
2052  int16_t *values = NULL;
2053  double value = 0.;
2054 
2055  /* this makes use of GDAL's "natural" blocks */
2056  GDALGetBlockSize(band, &nXBlockSize, &nYBlockSize);
2057  nXBlocks = (width + nXBlockSize - 1) / nXBlockSize;
2058  nYBlocks = (height + nYBlockSize - 1) / nYBlockSize;
2059  RASTER_DEBUGF(4, "(nXBlockSize, nYBlockSize) = (%d, %d)", nXBlockSize, nYBlockSize);
2060  RASTER_DEBUGF(4, "(nXBlocks, nYBlocks) = (%d, %d)", nXBlocks, nYBlocks);
2061 
2062  /* length is for the desired pixel type */
2063  valueslen = rt_pixtype_size(PT_16BSI) * nXBlockSize * nYBlockSize;
2064  values = rtalloc(valueslen);
2065  if (NULL == values) {
2066  rterror("rt_raster_to_gdal_mem: Could not allocate memory for GDAL band pixel values");
2067  if (allocBandNums) rtdealloc(bandNums);
2068  if (allocNodataValues) rtdealloc(excludeNodataValues);
2069  GDALClose(ds);
2070  return 0;
2071  }
2072 
2073  for (iYBlock = 0; iYBlock < nYBlocks; iYBlock++) {
2074  for (iXBlock = 0; iXBlock < nXBlocks; iXBlock++) {
2075  memset(values, 0, valueslen);
2076 
2077  x = iXBlock * nXBlockSize;
2078  y = iYBlock * nYBlockSize;
2079  RASTER_DEBUGF(4, "(iXBlock, iYBlock) = (%d, %d)", iXBlock, iYBlock);
2080  RASTER_DEBUGF(4, "(x, y) = (%d, %d)", x, y);
2081 
2082  /* valid block width */
2083  if ((iXBlock + 1) * nXBlockSize > width)
2084  nXValid = width - (iXBlock * nXBlockSize);
2085  else
2086  nXValid = nXBlockSize;
2087 
2088  /* valid block height */
2089  if ((iYBlock + 1) * nYBlockSize > height)
2090  nYValid = height - (iYBlock * nYBlockSize);
2091  else
2092  nYValid = nYBlockSize;
2093 
2094  RASTER_DEBUGF(4, "(nXValid, nYValid) = (%d, %d)", nXValid, nYValid);
2095 
2096  /* convert 8BSI values to 16BSI */
2097  z = 0;
2098  iYMax = y + nYValid;
2099  iXMax = x + nXValid;
2100  for (iY = y; iY < iYMax; iY++) {
2101  for (iX = x; iX < iXMax; iX++) {
2102  if (rt_band_get_pixel(rtband, iX, iY, &value, NULL) != ES_NONE) {
2103  rterror("rt_raster_to_gdal_mem: Could not get pixel value to convert from 8BSI to 16BSI");
2104  rtdealloc(values);
2105  if (allocBandNums) rtdealloc(bandNums);
2106  if (allocNodataValues) rtdealloc(excludeNodataValues);
2107  GDALClose(ds);
2108  return 0;
2109  }
2110 
2111  values[z++] = rt_util_clamp_to_16BSI(value);
2112  }
2113  }
2114 
2115  /* burn values */
2116  if (GDALRasterIO(
2117  band, GF_Write,
2118  x, y,
2119  nXValid, nYValid,
2120  values, nXValid, nYValid,
2121  gdal_pt,
2122  0, 0
2123  ) != CE_None) {
2124  rterror("rt_raster_to_gdal_mem: Could not write converted 8BSI to 16BSI values to GDAL band");
2125  rtdealloc(values);
2126  if (allocBandNums) rtdealloc(bandNums);
2127  if (allocNodataValues) rtdealloc(excludeNodataValues);
2128  GDALClose(ds);
2129  return 0;
2130  }
2131  }
2132  }
2133 
2134  rtdealloc(values);
2135  }
2136 
2137  /* Add nodata value for band */
2138  if (rt_band_get_hasnodata_flag(rtband) != FALSE && excludeNodataValues[i]) {
2139  rt_band_get_nodata(rtband, &nodata);
2140  if (GDALSetRasterNoDataValue(band, nodata) != CE_None)
2141  rtwarn("rt_raster_to_gdal_mem: Could not set nodata value for band");
2142  RASTER_DEBUGF(3, "nodata value set to %f", GDALGetRasterNoDataValue(band, NULL));
2143  }
2144 
2145 #if POSTGIS_DEBUG_LEVEL > 3
2146  {
2147  GDALRasterBandH _grb = NULL;
2148  double _min;
2149  double _max;
2150  double _mean;
2151  double _stddev;
2152 
2153  _grb = GDALGetRasterBand(ds, i + 1);
2154  GDALComputeRasterStatistics(_grb, FALSE, &_min, &_max, &_mean, &_stddev, NULL, NULL);
2155  RASTER_DEBUGF(4, "GDAL Band %d stats: %f, %f, %f, %f", i + 1, _min, _max, _mean, _stddev);
2156  }
2157 #endif
2158 
2159  }
2160 
2161  /* necessary??? */
2162  GDALFlushCache(ds);
2163 
2164  if (allocBandNums) rtdealloc(bandNums);
2165  if (allocNodataValues) rtdealloc(excludeNodataValues);
2166 
2167  return ds;
2168 }
2169 
2170 /******************************************************************************
2171 * rt_raster_from_gdal_dataset()
2172 ******************************************************************************/
2173 
2181 rt_raster
2183  rt_raster rast = NULL;
2184  double gt[6] = {0};
2185  CPLErr cplerr;
2186  uint32_t width = 0;
2187  uint32_t height = 0;
2188  uint32_t numBands = 0;
2189  uint32_t i = 0;
2190  char *authname = NULL;
2191  char *authcode = NULL;
2192 
2193  GDALRasterBandH gdband = NULL;
2194  GDALDataType gdpixtype = GDT_Unknown;
2195  rt_band band;
2196  int32_t idx;
2197  rt_pixtype pt = PT_END;
2198  uint32_t ptlen = 0;
2199  int hasnodata = 0;
2200  double nodataval;
2201 
2202  int x;
2203  int y;
2204 
2205  uint32_t nXBlocks, nYBlocks;
2206  int nXBlockSize, nYBlockSize;
2207  uint32_t iXBlock, iYBlock;
2208  uint32_t nXValid, nYValid;
2209  uint32_t iY;
2210 
2211  uint8_t *values = NULL;
2212  uint32_t valueslen = 0;
2213  uint8_t *ptr = NULL;
2214 
2215  assert(NULL != ds);
2216 
2217  /* raster size */
2218  width = GDALGetRasterXSize(ds);
2219  height = GDALGetRasterYSize(ds);
2220  RASTER_DEBUGF(3, "Raster dimensions (width x height): %d x %d", width, height);
2221 
2222  /* create new raster */
2223  RASTER_DEBUG(3, "Creating new raster");
2224  rast = rt_raster_new(width, height);
2225  if (NULL == rast) {
2226  rterror("rt_raster_from_gdal_dataset: Out of memory allocating new raster");
2227  return NULL;
2228  }
2229  RASTER_DEBUGF(3, "Created raster dimensions (width x height): %d x %d", rast->width, rast->height);
2230 
2231  /* get raster attributes */
2232  cplerr = GDALGetGeoTransform(ds, gt);
2233  if (GDALGetGeoTransform(ds, gt) != CE_None) {
2234  RASTER_DEBUG(4, "Using default geotransform matrix (0, 1, 0, 0, 0, -1)");
2235  gt[0] = 0;
2236  gt[1] = 1;
2237  gt[2] = 0;
2238  gt[3] = 0;
2239  gt[4] = 0;
2240  gt[5] = -1;
2241  }
2242 
2243  /* apply raster attributes */
2245 
2246  RASTER_DEBUGF(3, "Raster geotransform (%f, %f, %f, %f, %f, %f)",
2247  gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]);
2248 
2249  /* srid */
2250  if (rt_util_gdal_sr_auth_info(ds, &authname, &authcode) == ES_NONE) {
2251  if (
2252  authname != NULL &&
2253  strcmp(authname, "EPSG") == 0 &&
2254  authcode != NULL
2255  ) {
2256  rt_raster_set_srid(rast, atoi(authcode));
2257  RASTER_DEBUGF(3, "New raster's SRID = %d", rast->srid);
2258  }
2259 
2260  if (authname != NULL)
2261  rtdealloc(authname);
2262  if (authcode != NULL)
2263  rtdealloc(authcode);
2264  }
2265 
2266  numBands = GDALGetRasterCount(ds);
2267 
2268 #if POSTGIS_DEBUG_LEVEL > 3
2269  for (i = 1; i <= numBands; i++) {
2270  GDALRasterBandH _grb = NULL;
2271  double _min;
2272  double _max;
2273  double _mean;
2274  double _stddev;
2275 
2276  _grb = GDALGetRasterBand(ds, i);
2277  GDALComputeRasterStatistics(_grb, FALSE, &_min, &_max, &_mean, &_stddev, NULL, NULL);
2278  RASTER_DEBUGF(4, "GDAL Band %d stats: %f, %f, %f, %f", i, _min, _max, _mean, _stddev);
2279  }
2280 #endif
2281 
2282  /* copy bands */
2283  for (i = 1; i <= numBands; i++) {
2284  RASTER_DEBUGF(3, "Processing band %d of %d", i, numBands);
2285  gdband = NULL;
2286  gdband = GDALGetRasterBand(ds, i);
2287  if (NULL == gdband) {
2288  rterror("rt_raster_from_gdal_dataset: Could not get GDAL band");
2290  return NULL;
2291  }
2292  RASTER_DEBUGF(4, "gdband @ %p", gdband);
2293 
2294  /* pixtype */
2295  gdpixtype = GDALGetRasterDataType(gdband);
2296  RASTER_DEBUGF(4, "gdpixtype, size = %s, %d", GDALGetDataTypeName(gdpixtype), GDALGetDataTypeSize(gdpixtype) / 8);
2297  pt = rt_util_gdal_datatype_to_pixtype(gdpixtype);
2298  if (pt == PT_END) {
2299  rterror("rt_raster_from_gdal_dataset: Unknown pixel type for GDAL band");
2301  return NULL;
2302  }
2303  ptlen = rt_pixtype_size(pt);
2304 
2305  /* size: width and height */
2306  width = GDALGetRasterBandXSize(gdband);
2307  height = GDALGetRasterBandYSize(gdband);
2308  RASTER_DEBUGF(3, "GDAL band dimensions (width x height): %d x %d", width, height);
2309 
2310  /* nodata */
2311  nodataval = GDALGetRasterNoDataValue(gdband, &hasnodata);
2312  RASTER_DEBUGF(3, "(hasnodata, nodataval) = (%d, %f)", hasnodata, nodataval);
2313 
2314  /* create band object */
2316  rast, pt,
2317  (hasnodata ? nodataval : 0),
2318  hasnodata, nodataval, rt_raster_get_num_bands(rast)
2319  );
2320  if (idx < 0) {
2321  rterror("rt_raster_from_gdal_dataset: Could not allocate memory for raster band");
2323  return NULL;
2324  }
2325  band = rt_raster_get_band(rast, idx);
2326  RASTER_DEBUGF(3, "Created band of dimension (width x height): %d x %d", band->width, band->height);
2327 
2328  /* this makes use of GDAL's "natural" blocks */
2329  GDALGetBlockSize(gdband, &nXBlockSize, &nYBlockSize);
2330  nXBlocks = (width + nXBlockSize - 1) / nXBlockSize;
2331  nYBlocks = (height + nYBlockSize - 1) / nYBlockSize;
2332  RASTER_DEBUGF(4, "(nXBlockSize, nYBlockSize) = (%d, %d)", nXBlockSize, nYBlockSize);
2333  RASTER_DEBUGF(4, "(nXBlocks, nYBlocks) = (%d, %d)", nXBlocks, nYBlocks);
2334 
2335  /* allocate memory for values */
2336  valueslen = ptlen * nXBlockSize * nYBlockSize;
2337  values = rtalloc(valueslen);
2338  if (values == NULL) {
2339  rterror("rt_raster_from_gdal_dataset: Could not allocate memory for GDAL band pixel values");
2341  return NULL;
2342  }
2343  RASTER_DEBUGF(3, "values @ %p of length = %d", values, valueslen);
2344 
2345  for (iYBlock = 0; iYBlock < nYBlocks; iYBlock++) {
2346  for (iXBlock = 0; iXBlock < nXBlocks; iXBlock++) {
2347  x = iXBlock * nXBlockSize;
2348  y = iYBlock * nYBlockSize;
2349  RASTER_DEBUGF(4, "(iXBlock, iYBlock) = (%d, %d)", iXBlock, iYBlock);
2350  RASTER_DEBUGF(4, "(x, y) = (%d, %d)", x, y);
2351 
2352  memset(values, 0, valueslen);
2353 
2354  /* valid block width */
2355  if ((iXBlock + 1) * nXBlockSize > width)
2356  nXValid = width - (iXBlock * nXBlockSize);
2357  else
2358  nXValid = nXBlockSize;
2359 
2360  /* valid block height */
2361  if ((iYBlock + 1) * nYBlockSize > height)
2362  nYValid = height - (iYBlock * nYBlockSize);
2363  else
2364  nYValid = nYBlockSize;
2365 
2366  RASTER_DEBUGF(4, "(nXValid, nYValid) = (%d, %d)", nXValid, nYValid);
2367 
2368  cplerr = GDALRasterIO(
2369  gdband, GF_Read,
2370  x, y,
2371  nXValid, nYValid,
2372  values, nXValid, nYValid,
2373  gdpixtype,
2374  0, 0
2375  );
2376  if (cplerr != CE_None) {
2377  rterror("rt_raster_from_gdal_dataset: Could not get data from GDAL raster");
2378  rtdealloc(values);
2380  return NULL;
2381  }
2382 
2383  /* if block width is same as raster width, shortcut */
2384  if (nXBlocks == 1 && nYBlockSize > 1 && nXValid == width) {
2385  x = 0;
2386  y = nYBlockSize * iYBlock;
2387 
2388  RASTER_DEBUGF(4, "Setting set of pixel lines at (%d, %d) for %d pixels", x, y, nXValid * nYValid);
2389  rt_band_set_pixel_line(band, x, y, values, nXValid * nYValid);
2390  }
2391  else {
2392  ptr = values;
2393  x = nXBlockSize * iXBlock;
2394  for (iY = 0; iY < nYValid; iY++) {
2395  y = iY + (nYBlockSize * iYBlock);
2396 
2397  RASTER_DEBUGF(4, "Setting pixel line at (%d, %d) for %d pixels", x, y, nXValid);
2398  rt_band_set_pixel_line(band, x, y, ptr, nXValid);
2399  ptr += (nXValid * ptlen);
2400  }
2401  }
2402  }
2403  }
2404 
2405  /* free memory */
2406  rtdealloc(values);
2407  }
2408 
2409  return rast;
2410 }
2411 
2412 /******************************************************************************
2413 * rt_raster_gdal_rasterize()
2414 ******************************************************************************/
2415 
2419 
2421 
2422  OGRSpatialReferenceH src_sr;
2423 
2425  double *init;
2426  double *nodata;
2428  double *value;
2429  int *bandlist;
2430 };
2431 
2432 static _rti_rasterize_arg
2434  _rti_rasterize_arg arg = NULL;
2435 
2436  arg = rtalloc(sizeof(struct _rti_rasterize_arg_t));
2437  if (arg == NULL) {
2438  rterror("_rti_rasterize_arg_init: Could not allocate memory for _rti_rasterize_arg");
2439  return NULL;
2440  }
2441 
2442  arg->noband = 0;
2443 
2444  arg->numbands = 0;
2445 
2446  arg->src_sr = NULL;
2447 
2448  arg->pixtype = NULL;
2449  arg->init = NULL;
2450  arg->nodata = NULL;
2451  arg->hasnodata = NULL;
2452  arg->value = NULL;
2453  arg->bandlist = NULL;
2454 
2455  return arg;
2456 }
2457 
2458 static void
2460  if (arg->noband) {
2461  if (arg->pixtype != NULL)
2462  rtdealloc(arg->pixtype);
2463  if (arg->init != NULL)
2464  rtdealloc(arg->init);
2465  if (arg->nodata != NULL)
2466  rtdealloc(arg->nodata);
2467  if (arg->hasnodata != NULL)
2468  rtdealloc(arg->hasnodata);
2469  if (arg->value != NULL)
2470  rtdealloc(arg->value);
2471  }
2472 
2473  if (arg->bandlist != NULL)
2474  rtdealloc(arg->bandlist);
2475 
2476  if (arg->src_sr != NULL)
2477  OSRDestroySpatialReference(arg->src_sr);
2478 
2479  rtdealloc(arg);
2480 }
2481 
2508 rt_raster
2510  const unsigned char *wkb, uint32_t wkb_len,
2511  const char *srs,
2512  uint32_t num_bands, rt_pixtype *pixtype,
2513  double *init, double *value,
2514  double *nodata, uint8_t *hasnodata,
2515  int *width, int *height,
2516  double *scale_x, double *scale_y,
2517  double *ul_xw, double *ul_yw,
2518  double *grid_xw, double *grid_yw,
2519  double *skew_x, double *skew_y,
2520  char **options
2521 ) {
2522  rt_raster rast = NULL;
2523  uint32_t i = 0;
2524  int err = 0;
2525 
2526  _rti_rasterize_arg arg = NULL;
2527 
2528  int _dim[2] = {0};
2529  double _scale[2] = {0};
2530  double _skew[2] = {0};
2531 
2532  OGRErr ogrerr;
2533  OGRGeometryH src_geom;
2534  OGREnvelope src_env;
2535  rt_envelope extent;
2536  OGRwkbGeometryType wkbtype = wkbUnknown;
2537 
2538  int ul_user = 0;
2539 
2540  CPLErr cplerr;
2541  double _gt[6] = {0};
2542  GDALDriverH _drv = NULL;
2543  int unload_drv = 0;
2544  GDALDatasetH _ds = NULL;
2545  GDALRasterBandH _band = NULL;
2546 
2547  uint16_t _width = 0;
2548  uint16_t _height = 0;
2549 
2550  RASTER_DEBUG(3, "starting");
2551 
2552  assert(NULL != wkb);
2553  assert(0 != wkb_len);
2554 
2555  /* internal variables */
2556  arg = _rti_rasterize_arg_init();
2557  if (arg == NULL) {
2558  rterror("rt_raster_gdal_rasterize: Could not initialize internal variables");
2559  return NULL;
2560  }
2561 
2562  /* no bands, raster is a mask */
2563  if (num_bands < 1) {
2564  arg->noband = 1;
2565  arg->numbands = 1;
2566 
2567  arg->pixtype = (rt_pixtype *) rtalloc(sizeof(rt_pixtype));
2568  arg->pixtype[0] = PT_8BUI;
2569 
2570  arg->init = (double *) rtalloc(sizeof(double));
2571  arg->init[0] = 0;
2572 
2573  arg->nodata = (double *) rtalloc(sizeof(double));
2574  arg->nodata[0] = 0;
2575 
2576  arg->hasnodata = (uint8_t *) rtalloc(sizeof(uint8_t));
2577  arg->hasnodata[0] = 1;
2578 
2579  arg->value = (double *) rtalloc(sizeof(double));
2580  arg->value[0] = 1;
2581  }
2582  else {
2583  arg->noband = 0;
2584  arg->numbands = num_bands;
2585 
2586  arg->pixtype = pixtype;
2587  arg->init = init;
2588  arg->nodata = nodata;
2589  arg->hasnodata = hasnodata;
2590  arg->value = value;
2591  }
2592 
2593  /* OGR spatial reference */
2594  if (NULL != srs && strlen(srs)) {
2595  arg->src_sr = OSRNewSpatialReference(NULL);
2596  if (OSRSetFromUserInput(arg->src_sr, srs) != OGRERR_NONE) {
2597  rterror("rt_raster_gdal_rasterize: Could not create OSR spatial reference using the provided srs: %s", srs);
2599  return NULL;
2600  }
2601  }
2602 
2603  /* convert WKB to OGR Geometry */
2604  ogrerr = OGR_G_CreateFromWkb((unsigned char *) wkb, arg->src_sr, &src_geom, wkb_len);
2605  if (ogrerr != OGRERR_NONE) {
2606  rterror("rt_raster_gdal_rasterize: Could not create OGR Geometry from WKB");
2607 
2609  /* OGRCleanupAll(); */
2610 
2611  return NULL;
2612  }
2613 
2614  /* OGR Geometry is empty */
2615  if (OGR_G_IsEmpty(src_geom)) {
2616  rtinfo("Geometry provided is empty. Returning empty raster");
2617 
2618  OGR_G_DestroyGeometry(src_geom);
2620  /* OGRCleanupAll(); */
2621 
2622  return rt_raster_new(0, 0);
2623  }
2624 
2625  /* get envelope */
2626  OGR_G_GetEnvelope(src_geom, &src_env);
2627  rt_util_from_ogr_envelope(src_env, &extent);
2628 
2629  RASTER_DEBUGF(3, "Suggested raster envelope: %f, %f, %f, %f",
2630  extent.MinX, extent.MinY, extent.MaxX, extent.MaxY);
2631 
2632  /* user-defined scale */
2633  if (
2634  (NULL != scale_x) &&
2635  (NULL != scale_y) &&
2636  (FLT_NEQ(*scale_x, 0.0)) &&
2637  (FLT_NEQ(*scale_y, 0.0))
2638  ) {
2639  /* for now, force scale to be in left-right, top-down orientation */
2640  _scale[0] = fabs(*scale_x);
2641  _scale[1] = fabs(*scale_y);
2642  }
2643  /* user-defined width/height */
2644  else if (
2645  (NULL != width) &&
2646  (NULL != height) &&
2647  (FLT_NEQ(*width, 0.0)) &&
2648  (FLT_NEQ(*height, 0.0))
2649  ) {
2650  _dim[0] = abs(*width);
2651  _dim[1] = abs(*height);
2652 
2653  if (FLT_NEQ(extent.MaxX, extent.MinX))
2654  _scale[0] = fabs((extent.MaxX - extent.MinX) / _dim[0]);
2655  else
2656  _scale[0] = 1.;
2657 
2658  if (FLT_NEQ(extent.MaxY, extent.MinY))
2659  _scale[1] = fabs((extent.MaxY - extent.MinY) / _dim[1]);
2660  else
2661  _scale[1] = 1.;
2662  }
2663  else {
2664  rterror("rt_raster_gdal_rasterize: Values must be provided for width and height or X and Y of scale");
2665 
2666  OGR_G_DestroyGeometry(src_geom);
2668  /* OGRCleanupAll(); */
2669 
2670  return NULL;
2671  }
2672  RASTER_DEBUGF(3, "scale (x, y) = %f, %f", _scale[0], -1 * _scale[1]);
2673  RASTER_DEBUGF(3, "dim (x, y) = %d, %d", _dim[0], _dim[1]);
2674 
2675  /* user-defined skew */
2676  if (NULL != skew_x) {
2677  _skew[0] = *skew_x;
2678 
2679  /*
2680  negative scale-x affects skew
2681  for now, force skew to be in left-right, top-down orientation
2682  */
2683  if (
2684  NULL != scale_x &&
2685  *scale_x < 0.
2686  ) {
2687  _skew[0] *= -1;
2688  }
2689  }
2690  if (NULL != skew_y) {
2691  _skew[1] = *skew_y;
2692 
2693  /*
2694  positive scale-y affects skew
2695  for now, force skew to be in left-right, top-down orientation
2696  */
2697  if (
2698  NULL != scale_y &&
2699  *scale_y > 0.
2700  ) {
2701  _skew[1] *= -1;
2702  }
2703  }
2704 
2705  /*
2706  if geometry is a point, a linestring or set of either and bounds not set,
2707  increase extent by a pixel to avoid missing points on border
2708 
2709  a whole pixel is used instead of half-pixel due to backward
2710  compatibility with GDAL 1.6, 1.7 and 1.8. 1.9+ works fine with half-pixel.
2711  */
2712  wkbtype = wkbFlatten(OGR_G_GetGeometryType(src_geom));
2713  if ((
2714  (wkbtype == wkbPoint) ||
2715  (wkbtype == wkbMultiPoint) ||
2716  (wkbtype == wkbLineString) ||
2717  (wkbtype == wkbMultiLineString)
2718  ) &&
2719  _dim[0] == 0 &&
2720  _dim[1] == 0
2721  ) {
2722 
2723 #if POSTGIS_GDAL_VERSION > 18
2724 
2725  RASTER_DEBUG(3, "Adjusting extent for GDAL > 1.8 by half the scale on X-axis");
2726  extent.MinX -= (_scale[0] / 2.);
2727  extent.MaxX += (_scale[0] / 2.);
2728 
2729  RASTER_DEBUG(3, "Adjusting extent for GDAL > 1.8 by half the scale on Y-axis");
2730  extent.MinY -= (_scale[1] / 2.);
2731  extent.MaxY += (_scale[1] / 2.);
2732 
2733 #else
2734 
2735  RASTER_DEBUG(3, "Adjusting extent for GDAL <= 1.8 by the scale on X-axis");
2736  extent.MinX -= _scale[0];
2737  extent.MaxX += _scale[0];
2738 
2739  RASTER_DEBUG(3, "Adjusting extent for GDAL <= 1.8 by the scale on Y-axis");
2740  extent.MinY -= _scale[1];
2741  extent.MaxY += _scale[1];
2742 
2743 #endif
2744 
2745 
2746  RASTER_DEBUGF(3, "Adjusted extent: %f, %f, %f, %f",
2747  extent.MinX, extent.MinY, extent.MaxX, extent.MaxY);
2748 
2749  extent.UpperLeftX = extent.MinX;
2750  extent.UpperLeftY = extent.MaxY;
2751  }
2752 
2753  /* reprocess extent if skewed */
2754  if (
2755  FLT_NEQ(_skew[0], 0) ||
2756  FLT_NEQ(_skew[1], 0)
2757  ) {
2758  rt_raster skewedrast;
2759 
2760  RASTER_DEBUG(3, "Computing skewed extent's envelope");
2761 
2762  skewedrast = rt_raster_compute_skewed_raster(
2763  extent,
2764  _skew,
2765  _scale,
2766  0.01
2767  );
2768  if (skewedrast == NULL) {
2769  rterror("rt_raster_gdal_rasterize: Could not compute skewed raster");
2770 
2771  OGR_G_DestroyGeometry(src_geom);
2773  /* OGRCleanupAll(); */
2774 
2775  return NULL;
2776  }
2777 
2778  _dim[0] = skewedrast->width;
2779  _dim[1] = skewedrast->height;
2780 
2781  extent.UpperLeftX = skewedrast->ipX;
2782  extent.UpperLeftY = skewedrast->ipY;
2783 
2784  rt_raster_destroy(skewedrast);
2785  }
2786 
2787  /* raster dimensions */
2788  if (!_dim[0])
2789  _dim[0] = (int) fmax((fabs(extent.MaxX - extent.MinX) + (_scale[0] / 2.)) / _scale[0], 1);
2790  if (!_dim[1])
2791  _dim[1] = (int) fmax((fabs(extent.MaxY - extent.MinY) + (_scale[1] / 2.)) / _scale[1], 1);
2792 
2793  /* temporary raster */
2794  rast = rt_raster_new(_dim[0], _dim[1]);
2795  if (rast == NULL) {
2796  rterror("rt_raster_gdal_rasterize: Out of memory allocating temporary raster");
2797 
2798  OGR_G_DestroyGeometry(src_geom);
2800  /* OGRCleanupAll(); */
2801 
2802  return NULL;
2803  }
2804 
2805  /* set raster's spatial attributes */
2807  rt_raster_set_scale(rast, _scale[0], -1 * _scale[1]);
2808  rt_raster_set_skews(rast, _skew[0], _skew[1]);
2809 
2811  RASTER_DEBUGF(3, "Temp raster's geotransform: %f, %f, %f, %f, %f, %f",
2812  _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
2813  RASTER_DEBUGF(3, "Temp raster's dimensions (width x height): %d x %d",
2814  _dim[0], _dim[1]);
2815 
2816  /* user-specified upper-left corner */
2817  if (
2818  NULL != ul_xw &&
2819  NULL != ul_yw
2820  ) {
2821  ul_user = 1;
2822 
2823  RASTER_DEBUGF(4, "Using user-specified upper-left corner: %f, %f", *ul_xw, *ul_yw);
2824 
2825  /* set upper-left corner */
2826  rt_raster_set_offsets(rast, *ul_xw, *ul_yw);
2827  extent.UpperLeftX = *ul_xw;
2828  extent.UpperLeftY = *ul_yw;
2829  }
2830  else if (
2831  ((NULL != ul_xw) && (NULL == ul_yw)) ||
2832  ((NULL == ul_xw) && (NULL != ul_yw))
2833  ) {
2834  rterror("rt_raster_gdal_rasterize: Both X and Y upper-left corner values must be provided");
2835 
2837  OGR_G_DestroyGeometry(src_geom);
2839  /* OGRCleanupAll(); */
2840 
2841  return NULL;
2842  }
2843 
2844  /* alignment only considered if upper-left corner not provided */
2845  if (
2846  !ul_user && (
2847  (NULL != grid_xw) || (NULL != grid_yw)
2848  )
2849  ) {
2850 
2851  if (
2852  ((NULL != grid_xw) && (NULL == grid_yw)) ||
2853  ((NULL == grid_xw) && (NULL != grid_yw))
2854  ) {
2855  rterror("rt_raster_gdal_rasterize: Both X and Y alignment values must be provided");
2856 
2858  OGR_G_DestroyGeometry(src_geom);
2860  /* OGRCleanupAll(); */
2861 
2862  return NULL;
2863  }
2864 
2865  RASTER_DEBUGF(4, "Aligning extent to user-specified grid: %f, %f", *grid_xw, *grid_yw);
2866 
2867  do {
2868  double _r[2] = {0};
2869  double _w[2] = {0};
2870 
2871  /* raster is already aligned */
2872  if (FLT_EQ(*grid_xw, extent.UpperLeftX) && FLT_EQ(*grid_yw, extent.UpperLeftY)) {
2873  RASTER_DEBUG(3, "Skipping raster alignment as it is already aligned to grid");
2874  break;
2875  }
2876 
2877  extent.UpperLeftX = rast->ipX;
2878  extent.UpperLeftY = rast->ipY;
2879  rt_raster_set_offsets(rast, *grid_xw, *grid_yw);
2880 
2881  /* process upper-left corner */
2883  rast,
2884  extent.UpperLeftX, extent.UpperLeftY,
2885  &(_r[0]), &(_r[1]),
2886  NULL
2887  ) != ES_NONE) {
2888  rterror("rt_raster_gdal_rasterize: Could not compute raster pixel for spatial coordinates");
2889 
2891  OGR_G_DestroyGeometry(src_geom);
2893  /* OGRCleanupAll(); */
2894 
2895  return NULL;
2896  }
2897 
2899  rast,
2900  _r[0], _r[1],
2901  &(_w[0]), &(_w[1]),
2902  NULL
2903  ) != ES_NONE) {
2904  rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
2905 
2907  OGR_G_DestroyGeometry(src_geom);
2909  /* OGRCleanupAll(); */
2910 
2911  return NULL;
2912  }
2913 
2914  /* shift occurred */
2915  if (FLT_NEQ(_w[0], extent.UpperLeftX)) {
2916  if (NULL == width)
2917  rast->width++;
2918  else if (NULL == scale_x) {
2919  double _c[2] = {0};
2920 
2922 
2923  /* get upper-right corner */
2925  rast,
2926  rast->width, 0,
2927  &(_c[0]), &(_c[1]),
2928  NULL
2929  ) != ES_NONE) {
2930  rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
2931 
2933  OGR_G_DestroyGeometry(src_geom);
2935  /* OGRCleanupAll(); */
2936 
2937  return NULL;
2938  }
2939 
2940  rast->scaleX = fabs((_c[0] - _w[0]) / ((double) rast->width));
2941  }
2942  }
2943  if (FLT_NEQ(_w[1], extent.UpperLeftY)) {
2944  if (NULL == height)
2945  rast->height++;
2946  else if (NULL == scale_y) {
2947  double _c[2] = {0};
2948 
2950 
2951  /* get upper-right corner */
2953  rast,
2954  0, rast->height,
2955  &(_c[0]), &(_c[1]),
2956  NULL
2957  ) != ES_NONE) {
2958  rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
2959 
2961  OGR_G_DestroyGeometry(src_geom);
2963  /* OGRCleanupAll(); */
2964 
2965  return NULL;
2966  }
2967 
2968  rast->scaleY = -1 * fabs((_c[1] - _w[1]) / ((double) rast->height));
2969  }
2970  }
2971 
2972  rt_raster_set_offsets(rast, _w[0], _w[1]);
2973  }
2974  while (0);
2975  }
2976 
2977  /*
2978  after this point, rt_envelope extent is no longer used
2979  */
2980 
2981  /* get key attributes from rast */
2982  _dim[0] = rast->width;
2983  _dim[1] = rast->height;
2985 
2986  /* scale-x is negative or scale-y is positive */
2987  if ((
2988  (NULL != scale_x) && (*scale_x < 0.)
2989  ) || (
2990  (NULL != scale_y) && (*scale_y > 0)
2991  )) {
2992  double _w[2] = {0};
2993 
2994  /* negative scale-x */
2995  if (
2996  (NULL != scale_x) &&
2997  (*scale_x < 0.)
2998  ) {
2999  RASTER_DEBUG(3, "Processing negative scale-x");
3000 
3002  rast,
3003  _dim[0], 0,
3004  &(_w[0]), &(_w[1]),
3005  NULL
3006  ) != ES_NONE) {
3007  rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
3008 
3010  OGR_G_DestroyGeometry(src_geom);
3012  /* OGRCleanupAll(); */
3013 
3014  return NULL;
3015  }
3016 
3017  _gt[0] = _w[0];
3018  _gt[1] = *scale_x;
3019 
3020  /* check for skew */
3021  if (NULL != skew_x && FLT_NEQ(*skew_x, 0))
3022  _gt[2] = *skew_x;
3023  }
3024  /* positive scale-y */
3025  if (
3026  (NULL != scale_y) &&
3027  (*scale_y > 0)
3028  ) {
3029  RASTER_DEBUG(3, "Processing positive scale-y");
3030 
3032  rast,
3033  0, _dim[1],
3034  &(_w[0]), &(_w[1]),
3035  NULL
3036  ) != ES_NONE) {
3037  rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
3038 
3040  OGR_G_DestroyGeometry(src_geom);
3042  /* OGRCleanupAll(); */
3043 
3044  return NULL;
3045  }
3046 
3047  _gt[3] = _w[1];
3048  _gt[5] = *scale_y;
3049 
3050  /* check for skew */
3051  if (NULL != skew_y && FLT_NEQ(*skew_y, 0))
3052  _gt[4] = *skew_y;
3053  }
3054  }
3055 
3057  rast = NULL;
3058 
3059  RASTER_DEBUGF(3, "Applied geotransform: %f, %f, %f, %f, %f, %f",
3060  _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
3061  RASTER_DEBUGF(3, "Raster dimensions (width x height): %d x %d",
3062  _dim[0], _dim[1]);
3063 
3064  /* load GDAL mem */
3065  if (!rt_util_gdal_driver_registered("MEM")) {
3066  RASTER_DEBUG(4, "Registering MEM driver");
3067  GDALRegister_MEM();
3068  unload_drv = 1;
3069  }
3070  _drv = GDALGetDriverByName("MEM");
3071  if (NULL == _drv) {
3072  rterror("rt_raster_gdal_rasterize: Could not load the MEM GDAL driver");
3073 
3074  OGR_G_DestroyGeometry(src_geom);
3076  /* OGRCleanupAll(); */
3077 
3078  return NULL;
3079  }
3080 
3081  /* unload driver from GDAL driver manager */
3082  if (unload_drv) {
3083  RASTER_DEBUG(4, "Deregistering MEM driver");
3084  GDALDeregisterDriver(_drv);
3085  }
3086 
3087  _ds = GDALCreate(_drv, "", _dim[0], _dim[1], 0, GDT_Byte, NULL);
3088  if (NULL == _ds) {
3089  rterror("rt_raster_gdal_rasterize: Could not create a GDALDataset to rasterize the geometry into");
3090 
3091  OGR_G_DestroyGeometry(src_geom);
3093  /* OGRCleanupAll(); */
3094  if (unload_drv) GDALDestroyDriver(_drv);
3095 
3096  return NULL;
3097  }
3098 
3099  /* set geotransform */
3100  cplerr = GDALSetGeoTransform(_ds, _gt);
3101  if (cplerr != CE_None) {
3102  rterror("rt_raster_gdal_rasterize: Could not set geotransform on GDALDataset");
3103 
3104  OGR_G_DestroyGeometry(src_geom);
3106  /* OGRCleanupAll(); */
3107 
3108  GDALClose(_ds);
3109  if (unload_drv) GDALDestroyDriver(_drv);
3110 
3111  return NULL;
3112  }
3113 
3114  /* set SRS */
3115  if (NULL != arg->src_sr) {
3116  char *_srs = NULL;
3117  OSRExportToWkt(arg->src_sr, &_srs);
3118 
3119  cplerr = GDALSetProjection(_ds, _srs);
3120  CPLFree(_srs);
3121  if (cplerr != CE_None) {
3122  rterror("rt_raster_gdal_rasterize: Could not set projection on GDALDataset");
3123 
3124  OGR_G_DestroyGeometry(src_geom);
3126  /* OGRCleanupAll(); */
3127 
3128  GDALClose(_ds);
3129  if (unload_drv) GDALDestroyDriver(_drv);
3130 
3131  return NULL;
3132  }
3133  }
3134 
3135  /* set bands */
3136  for (i = 0; i < arg->numbands; i++) {
3137  err = 0;
3138 
3139  do {
3140  /* add band */
3141  cplerr = GDALAddBand(_ds, rt_util_pixtype_to_gdal_datatype(arg->pixtype[i]), NULL);
3142  if (cplerr != CE_None) {
3143  rterror("rt_raster_gdal_rasterize: Could not add band to GDALDataset");
3144  err = 1;
3145  break;
3146  }
3147 
3148  _band = GDALGetRasterBand(_ds, i + 1);
3149  if (NULL == _band) {
3150  rterror("rt_raster_gdal_rasterize: Could not get band %d from GDALDataset", i + 1);
3151  err = 1;
3152  break;
3153  }
3154 
3155  /* nodata value */
3156  if (arg->hasnodata[i]) {
3157  RASTER_DEBUGF(4, "Setting NODATA value of band %d to %f", i, arg->nodata[i]);
3158  cplerr = GDALSetRasterNoDataValue(_band, arg->nodata[i]);
3159  if (cplerr != CE_None) {
3160  rterror("rt_raster_gdal_rasterize: Could not set nodata value");
3161  err = 1;
3162  break;
3163  }
3164  RASTER_DEBUGF(4, "NODATA value set to %f", GDALGetRasterNoDataValue(_band, NULL));
3165  }
3166 
3167  /* initial value */
3168  cplerr = GDALFillRaster(_band, arg->init[i], 0);
3169  if (cplerr != CE_None) {
3170  rterror("rt_raster_gdal_rasterize: Could not set initial value");
3171  err = 1;
3172  break;
3173  }
3174  }
3175  while (0);
3176 
3177  if (err) {
3178 
3179  OGR_G_DestroyGeometry(src_geom);
3181 
3182  /* OGRCleanupAll(); */
3183 
3184  GDALClose(_ds);
3185  if (unload_drv) GDALDestroyDriver(_drv);
3186 
3187  return NULL;
3188  }
3189  }
3190 
3191  arg->bandlist = (int *) rtalloc(sizeof(int) * arg->numbands);
3192  for (i = 0; i < arg->numbands; i++) arg->bandlist[i] = i + 1;
3193 
3194  /* burn geometry */
3195  cplerr = GDALRasterizeGeometries(
3196  _ds,
3197  arg->numbands, arg->bandlist,
3198  1, &src_geom,
3199  NULL, NULL,
3200  arg->value,
3201  options,
3202  NULL, NULL
3203  );
3204  if (cplerr != CE_None) {
3205  rterror("rt_raster_gdal_rasterize: Could not rasterize geometry");
3206 
3207  OGR_G_DestroyGeometry(src_geom);
3209  /* OGRCleanupAll(); */
3210 
3211  GDALClose(_ds);
3212  if (unload_drv) GDALDestroyDriver(_drv);
3213 
3214  return NULL;
3215  }
3216 
3217  /* convert gdal dataset to raster */
3218  GDALFlushCache(_ds);
3219  RASTER_DEBUG(3, "Converting GDAL dataset to raster");
3221 
3222  OGR_G_DestroyGeometry(src_geom);
3223  /* OGRCleanupAll(); */
3224 
3225  GDALClose(_ds);
3226  if (unload_drv) GDALDestroyDriver(_drv);
3227 
3228  if (NULL == rast) {
3229  rterror("rt_raster_gdal_rasterize: Could not rasterize geometry");
3230  return NULL;
3231  }
3232 
3233  /* width, height */
3234  _width = rt_raster_get_width(rast);
3235  _height = rt_raster_get_height(rast);
3236 
3237  /* check each band for pixtype */
3238  for (i = 0; i < arg->numbands; i++) {
3239  uint8_t *data = NULL;
3240  rt_band band = NULL;
3241  rt_band oldband = NULL;
3242 
3243  double val = 0;
3244  int nodata = 0;
3245  int hasnodata = 0;
3246  double nodataval = 0;
3247  int x = 0;
3248  int y = 0;
3249 
3250  oldband = rt_raster_get_band(rast, i);
3251  if (oldband == NULL) {
3252  rterror("rt_raster_gdal_rasterize: Could not get band %d of output raster", i);
3255  return NULL;
3256  }
3257 
3258  /* band is of user-specified type */
3259  if (rt_band_get_pixtype(oldband) == arg->pixtype[i])
3260  continue;
3261 
3262  /* hasnodata, nodataval */
3263  hasnodata = rt_band_get_hasnodata_flag(oldband);
3264  if (hasnodata)
3265  rt_band_get_nodata(oldband, &nodataval);
3266 
3267  /* allocate data */
3268  data = rtalloc(rt_pixtype_size(arg->pixtype[i]) * _width * _height);
3269  if (data == NULL) {
3270  rterror("rt_raster_gdal_rasterize: Could not allocate memory for band data");
3273  return NULL;
3274  }
3275  memset(data, 0, rt_pixtype_size(arg->pixtype[i]) * _width * _height);
3276 
3277  /* create new band of correct type */
3279  _width, _height,
3280  arg->pixtype[i],
3281  hasnodata, nodataval,
3282  data
3283  );
3284  if (band == NULL) {
3285  rterror("rt_raster_gdal_rasterize: Could not create band");
3286  rtdealloc(data);
3289  return NULL;
3290  }
3291 
3292  /* give ownership of data to band */
3294 
3295  /* copy pixel by pixel */
3296  for (x = 0; x < _width; x++) {
3297  for (y = 0; y < _height; y++) {
3298  err = rt_band_get_pixel(oldband, x, y, &val, &nodata);
3299  if (err != ES_NONE) {
3300  rterror("rt_raster_gdal_rasterize: Could not get pixel value");
3304  return NULL;
3305  }
3306 
3307  if (nodata)
3308  val = nodataval;
3309 
3310  err = rt_band_set_pixel(band, x, y, val, NULL);
3311  if (err != ES_NONE) {
3312  rterror("rt_raster_gdal_rasterize: Could not set pixel value");
3316  return NULL;
3317  }
3318  }
3319  }
3320 
3321  /* replace band */
3322  oldband = rt_raster_replace_band(rast, band, i);
3323  if (oldband == NULL) {
3324  rterror("rt_raster_gdal_rasterize: Could not replace band %d of output raster", i);
3328  return NULL;
3329  }
3330 
3331  /* free oldband */
3332  rt_band_destroy(oldband);
3333  }
3334 
3336 
3337  RASTER_DEBUG(3, "done");
3338 
3339  return rast;
3340 }
3341 
3342 /******************************************************************************
3343 * rt_raster_from_two_rasters()
3344 ******************************************************************************/
3345 
3346 /*
3347  * Return raster of computed extent specified extenttype applied
3348  * on two input rasters. The raster returned should be freed by
3349  * the caller
3350  *
3351  * @param rast1 : the first raster
3352  * @param rast2 : the second raster
3353  * @param extenttype : type of extent for the output raster
3354  * @param *rtnraster : raster of computed extent
3355  * @param *offset : 4-element array indicating the X,Y offsets
3356  * for each raster. 0,1 for rast1 X,Y. 2,3 for rast2 X,Y.
3357  *
3358  * @return ES_NONE if success, ES_ERROR if error
3359  */
3362  rt_raster rast1, rt_raster rast2,
3363  rt_extenttype extenttype,
3364  rt_raster *rtnraster, double *offset
3365 ) {
3366  int i;
3367 
3368  rt_raster _rast[2] = {rast1, rast2};
3369  double _offset[2][4] = {{0.}};
3370  uint16_t _dim[2][2] = {{0}};
3371 
3372  rt_raster raster = NULL;
3373  int aligned = 0;
3374  int dim[2] = {0};
3375  double gt[6] = {0.};
3376 
3377  assert(NULL != rast1);
3378  assert(NULL != rast2);
3379  assert(NULL != rtnraster);
3380 
3381  /* set rtnraster to NULL */
3382  *rtnraster = NULL;
3383 
3384  /* rasters must be aligned */
3385  if (rt_raster_same_alignment(rast1, rast2, &aligned, NULL) != ES_NONE) {
3386  rterror("rt_raster_from_two_rasters: Could not test for alignment on the two rasters");
3387  return ES_ERROR;
3388  }
3389  if (!aligned) {
3390  rterror("rt_raster_from_two_rasters: The two rasters provided do not have the same alignment");
3391  return ES_ERROR;
3392  }
3393 
3394  /* dimensions */
3395  _dim[0][0] = rast1->width;
3396  _dim[0][1] = rast1->height;
3397  _dim[1][0] = rast2->width;
3398  _dim[1][1] = rast2->height;
3399 
3400  /* get raster offsets */
3402  _rast[1],
3403  _rast[0]->ipX, _rast[0]->ipY,
3404  &(_offset[1][0]), &(_offset[1][1]),
3405  NULL
3406  ) != ES_NONE) {
3407  rterror("rt_raster_from_two_rasters: Could not compute offsets of the second raster relative to the first raster");
3408  return ES_ERROR;
3409  }
3410  _offset[1][0] = -1 * _offset[1][0];
3411  _offset[1][1] = -1 * _offset[1][1];
3412  _offset[1][2] = _offset[1][0] + _dim[1][0] - 1;
3413  _offset[1][3] = _offset[1][1] + _dim[1][1] - 1;
3414 
3415  i = -1;
3416  switch (extenttype) {
3417  case ET_FIRST:
3418  i = 0;
3419  _offset[0][0] = 0.;
3420  _offset[0][1] = 0.;
3421  /* FALLTHROUGH */
3422  case ET_LAST:
3423  case ET_SECOND:
3424  if (i < 0) {
3425  i = 1;
3426  _offset[0][0] = -1 * _offset[1][0];
3427  _offset[0][1] = -1 * _offset[1][1];
3428  _offset[1][0] = 0.;
3429  _offset[1][1] = 0.;
3430  }
3431 
3432  dim[0] = _dim[i][0];
3433  dim[1] = _dim[i][1];
3435  dim[0],
3436  dim[1]
3437  );
3438  if (raster == NULL) {
3439  rterror("rt_raster_from_two_rasters: Could not create output raster");
3440  return ES_ERROR;
3441  }
3442  rt_raster_set_srid(raster, _rast[i]->srid);
3445  break;
3446  case ET_UNION: {
3447  double off[4] = {0};
3448 
3450  RASTER_DEBUGF(4, "gt = (%f, %f, %f, %f, %f, %f)",
3451  gt[0],
3452  gt[1],
3453  gt[2],
3454  gt[3],
3455  gt[4],
3456  gt[5]
3457  );
3458 
3459  /* new raster upper left offset */
3460  off[0] = 0;
3461  if (_offset[1][0] < 0)
3462  off[0] = _offset[1][0];
3463  off[1] = 0;
3464  if (_offset[1][1] < 0)
3465  off[1] = _offset[1][1];
3466 
3467  /* new raster lower right offset */
3468  off[2] = _dim[0][0] - 1;
3469  if ((int) _offset[1][2] >= _dim[0][0])
3470  off[2] = _offset[1][2];
3471  off[3] = _dim[0][1] - 1;
3472  if ((int) _offset[1][3] >= _dim[0][1])
3473  off[3] = _offset[1][3];
3474 
3475  /* upper left corner */
3477  _rast[0],
3478  off[0], off[1],
3479  &(gt[0]), &(gt[3]),
3480  NULL
3481  ) != ES_NONE) {
3482  rterror("rt_raster_from_two_rasters: Could not get spatial coordinates of upper-left pixel of output raster");
3483  return ES_ERROR;
3484  }
3485 
3486  dim[0] = off[2] - off[0] + 1;
3487  dim[1] = off[3] - off[1] + 1;
3488  RASTER_DEBUGF(4, "off = (%f, %f, %f, %f)",
3489  off[0],
3490  off[1],
3491  off[2],
3492  off[3]
3493  );
3494  RASTER_DEBUGF(4, "dim = (%d, %d)", dim[0], dim[1]);
3495 
3497  dim[0],
3498  dim[1]
3499  );
3500  if (raster == NULL) {
3501  rterror("rt_raster_from_two_rasters: Could not create output raster");
3502  return ES_ERROR;
3503  }
3504  rt_raster_set_srid(raster, _rast[0]->srid);
3506  RASTER_DEBUGF(4, "gt = (%f, %f, %f, %f, %f, %f)",
3507  gt[0],
3508  gt[1],
3509  gt[2],
3510  gt[3],
3511  gt[4],
3512  gt[5]
3513  );
3514 
3515  /* get offsets */
3517  _rast[0],
3518  gt[0], gt[3],
3519  &(_offset[0][0]), &(_offset[0][1]),
3520  NULL
3521  ) != ES_NONE) {
3522  rterror("rt_raster_from_two_rasters: Could not get offsets of the FIRST raster relative to the output raster");
3524  return ES_ERROR;
3525  }
3526  _offset[0][0] *= -1;
3527  _offset[0][1] *= -1;
3528 
3530  _rast[1],
3531  gt[0], gt[3],
3532  &(_offset[1][0]), &(_offset[1][1]),
3533  NULL
3534  ) != ES_NONE) {
3535  rterror("rt_raster_from_two_rasters: Could not get offsets of the SECOND raster relative to the output raster");
3537  return ES_ERROR;
3538  }
3539  _offset[1][0] *= -1;
3540  _offset[1][1] *= -1;
3541  break;
3542  }
3543  case ET_INTERSECTION: {
3544  double off[4] = {0};
3545 
3546  /* no intersection */
3547  if (
3548  (_offset[1][2] < 0 || _offset[1][0] > (_dim[0][0] - 1)) ||
3549  (_offset[1][3] < 0 || _offset[1][1] > (_dim[0][1] - 1))
3550  ) {
3551  RASTER_DEBUG(3, "The two rasters provided have no intersection. Returning no band raster");
3552 
3553  raster = rt_raster_new(0, 0);
3554  if (raster == NULL) {
3555  rterror("rt_raster_from_two_rasters: Could not create output raster");
3556  return ES_ERROR;
3557  }
3558  rt_raster_set_srid(raster, _rast[0]->srid);
3559  rt_raster_set_scale(raster, 0, 0);
3560 
3561  /* set offsets if provided */
3562  if (NULL != offset) {
3563  for (i = 0; i < 4; i++)
3564  offset[i] = _offset[i / 2][i % 2];
3565  }
3566 
3567  *rtnraster = raster;
3568  return ES_NONE;
3569  }
3570 
3571  if (_offset[1][0] > 0)
3572  off[0] = _offset[1][0];
3573  if (_offset[1][1] > 0)
3574  off[1] = _offset[1][1];
3575 
3576  off[2] = _dim[0][0] - 1;
3577  if (_offset[1][2] < _dim[0][0])
3578  off[2] = _offset[1][2];
3579  off[3] = _dim[0][1] - 1;
3580  if (_offset[1][3] < _dim[0][1])
3581  off[3] = _offset[1][3];
3582 
3583  dim[0] = off[2] - off[0] + 1;
3584  dim[1] = off[3] - off[1] + 1;
3586  dim[0],
3587  dim[1]
3588  );
3589  if (raster == NULL) {
3590  rterror("rt_raster_from_two_rasters: Could not create output raster");
3591  return ES_ERROR;
3592  }
3593  rt_raster_set_srid(raster, _rast[0]->srid);
3594 
3595  /* get upper-left corner */
3598  _rast[0],
3599  off[0], off[1],
3600  &(gt[0]), &(gt[3]),
3601  gt
3602  ) != ES_NONE) {
3603  rterror("rt_raster_from_two_rasters: Could not get spatial coordinates of upper-left pixel of output raster");
3605  return ES_ERROR;
3606  }
3607 
3609 
3610  /* get offsets */
3612  _rast[0],
3613  gt[0], gt[3],
3614  &(_offset[0][0]), &(_offset[0][1]),
3615  NULL
3616  ) != ES_NONE) {
3617  rterror("rt_raster_from_two_rasters: Could not get pixel coordinates to compute the offsets of the FIRST raster relative to the output raster");
3619  return ES_ERROR;
3620  }
3621  _offset[0][0] *= -1;
3622  _offset[0][1] *= -1;
3623 
3625  _rast[1],
3626  gt[0], gt[3],
3627  &(_offset[1][0]), &(_offset[1][1]),
3628  NULL
3629  ) != ES_NONE) {
3630  rterror("rt_raster_from_two_rasters: Could not get pixel coordinates to compute the offsets of the SECOND raster relative to the output raster");
3632  return ES_ERROR;
3633  }
3634  _offset[1][0] *= -1;
3635  _offset[1][1] *= -1;
3636  break;
3637  }
3638  case ET_CUSTOM:
3639  rterror("rt_raster_from_two_rasters: Extent type ET_CUSTOM is not supported by this function");
3640  break;
3641  }
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 }
return(const char *)
Definition: dbfopen.c:1613
#define TRUE
Definition: dbfopen.c:169
#define FALSE
Definition: dbfopen.c:168
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, uint8_t autofix)
void lwgeom_geos_error(const char *fmt,...)
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1144
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:320
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:347
void lwpoly_free(LWPOLY *poly)
Definition: lwpoly.c:175
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:188
rt_band rt_band_new_inline(uint16_t width, uint16_t height, rt_pixtype pixtype, uint32_t hasnodata, double nodataval, uint8_t *data)
Create an in-db rt_band with no data.
Definition: rt_band.c:63
void rt_band_set_ownsdata_flag(rt_band band, int flag)
Definition: rt_band.c:667
void rterror(const char *fmt,...)
Wrappers used for reporting errors and info.
Definition: rt_context.c:199
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition: rt_context.c:171
rt_pixtype rt_util_gdal_datatype_to_pixtype(GDALDataType gdt)
Convert GDALDataType to rt_pixtype.
Definition: rt_util.c:153
rt_band rt_band_duplicate(rt_band band)
Create a new band duplicated from source band.
Definition: rt_band.c:287
rt_errorstate rt_band_set_isnodata_flag(rt_band band, int flag)
Set isnodata flag value.
Definition: rt_band.c:695
#define FLT_NEQ(x, y)
Definition: librtcore.h:2233
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:295
int rt_util_gdal_register_all(int force_register_all)
Definition: rt_util.c:336
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:299
int8_t rt_util_clamp_to_8BSI(double value)
Definition: rt_util.c:49
uint8_t rt_util_clamp_to_1BB(double value)
Definition: rt_util.c:34
void rtinfo(const char *fmt,...)
Definition: rt_context.c:211
int32_t rt_util_clamp_to_32BSI(double value)
Definition: rt_util.c:69
char * rt_util_gdal_convert_sr(const char *srs, int proj4)
Definition: rt_util.c:212
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:674
#define ROUND(x, y)
Definition: librtcore.h:2241
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_band.c:1221
rt_pixtype
Definition: librtcore.h:185
@ PT_32BUI
Definition: librtcore.h:194
@ PT_2BUI
Definition: librtcore.h:187
@ PT_32BSI
Definition: librtcore.h:193
@ PT_END
Definition: librtcore.h:197
@ PT_4BUI
Definition: librtcore.h:188
@ PT_32BF
Definition: librtcore.h:195
@ PT_1BB
Definition: librtcore.h:186
@ PT_16BUI
Definition: librtcore.h:192
@ PT_8BSI
Definition: librtcore.h:189
@ PT_16BSI
Definition: librtcore.h:191
@ PT_64BF
Definition: librtcore.h:196
@ PT_8BUI
Definition: librtcore.h:190
int rt_util_dbl_trunc_warning(double initialvalue, int32_t checkvalint, uint32_t checkvaluint, float checkvalfloat, double checkvaldouble, rt_pixtype pixtype)
Definition: rt_util.c:627
GDALDataType rt_util_pixtype_to_gdal_datatype(rt_pixtype pt)
Convert rt_pixtype to GDALDataType.
Definition: rt_util.c:118
#define FLT_EQ(x, y)
Definition: librtcore.h:2234
uint8_t rt_util_clamp_to_2BUI(double value)
Definition: rt_util.c:39
void rtwarn(const char *fmt,...)
Definition: rt_context.c:224
uint8_t rt_util_clamp_to_8BUI(double value)
Definition: rt_util.c:54
rt_errorstate rt_band_set_pixel(rt_band band, int x, int y, double val, int *converted)
Set single pixel's value.
Definition: rt_band.c:974
rt_errorstate
Enum definitions.
Definition: librtcore.h:179
@ ES_NONE
Definition: librtcore.h:180
@ ES_ERROR
Definition: librtcore.h:181
rt_errorstate rt_band_set_pixel_line(rt_band band, int x, int y, void *vals, uint32_t len)
Set values of multiple pixels.
Definition: rt_band.c:853
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition: rt_band.c:340
int16_t rt_util_clamp_to_16BSI(double value)
Definition: rt_util.c:59
void * rtrealloc(void *mem, size_t size)
Definition: rt_context.c:179
struct rt_gdaldriver_t * rt_gdaldriver
Definition: librtcore.h:154
rt_errorstate rt_raster_get_convex_hull(rt_raster raster, LWGEOM **hull)
Get raster's convex hull.
Definition: rt_geometry.c:803
int rt_util_gdal_driver_registered(const char *drv)
Definition: rt_util.c:355
rt_extenttype
Definition: librtcore.h:200
@ ET_CUSTOM
Definition: librtcore.h:206
@ ET_LAST
Definition: librtcore.h:205
@ ET_INTERSECTION
Definition: librtcore.h:201
@ ET_UNION
Definition: librtcore.h:202
@ ET_SECOND
Definition: librtcore.h:204
@ ET_FIRST
Definition: librtcore.h:203
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1730
uint8_t rt_util_clamp_to_4BUI(double value)
Definition: rt_util.c:44
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition: rt_band.c:631
void rtdealloc(void *mem)
Definition: rt_context.c:186
void * rt_band_get_data(rt_band band)
Get pointer to raster band data.
Definition: rt_band.c:400
rt_errorstate rt_util_gdal_sr_auth_info(GDALDatasetH hds, char **authname, char **authcode)
Get auth name and code.
Definition: rt_util.c:270
uint16_t rt_util_clamp_to_16BUI(double value)
Definition: rt_util.c:64
uint32_t rt_util_clamp_to_32BUI(double value)
Definition: rt_util.c:74
int rt_band_is_offline(rt_band band)
Return non-zero if the given band data is on the filesystem.
Definition: rt_band.c:329
LWPOLY * rt_util_envelope_to_lwpoly(rt_envelope ext)
Definition: rt_util.c:437
struct rt_raster_t * rt_raster
Types definitions.
Definition: librtcore.h:145
float rt_util_clamp_to_32F(double value)
Definition: rt_util.c:79
rt_errorstate rt_raster_same_alignment(rt_raster rast1, rt_raster rast2, int *aligned, char **reason)
void rt_util_from_ogr_envelope(OGREnvelope env, rt_envelope *ext)
Definition: rt_util.c:408
int rt_pixtype_size(rt_pixtype pixtype)
Return size in bytes of a value in the given pixtype.
Definition: rt_pixel.c:39
This library is the generic raster handling section of PostGIS.
int value
Definition: genraster.py:61
int count
Definition: genraster.py:56
band
Definition: ovdump.py:57
src_ds
MAIN.
Definition: ovdump.py:53
data
Definition: ovdump.py:103
ds
Definition: pixval.py:66
nband
Definition: pixval.py:52
raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
gt
Definition: window.py:77
Datum covers(PG_FUNCTION_ARGS)
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
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition: rt_raster.c:356
double rt_raster_get_x_skew(rt_raster raster)
Get skew about the X axis.
Definition: rt_raster.c:181
double rt_raster_get_x_offset(rt_raster raster)
Get raster x offset, in projection units.
Definition: rt_raster.c:213
int rt_raster_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
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
void rt_raster_set_geotransform_matrix(rt_raster raster, double *gt)
Set raster's geotransform using 6-element array.
Definition: rt_raster.c:727
void rt_raster_set_scale(rt_raster raster, double scaleX, double scaleY)
Set scale in projection units.
Definition: rt_raster.c:137
int rt_raster_add_band(rt_raster raster, rt_band band, int index)
Add band data to a raster.
Definition: rt_raster.c:405
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
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:1598
rt_errorstate rt_raster_geopoint_to_cell(rt_raster raster, double xw, double yw, double *xr, double *yr, double *igt)
Convert an xw,yw map point to a xr,yr raster point.
Definition: rt_raster.c:806
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:82
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
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_raster.c:48
int rt_raster_has_band(rt_raster raster, int nband)
Return TRUE if the raster has a band of this number.
Definition: rt_raster.c:1347
rt_raster rt_raster_compute_skewed_raster(rt_envelope extent, double *skew, double *scale, double tolerance)
Definition: rt_raster.c:960
rt_gdaldriver rt_raster_gdal_drivers(uint32_t *drv_count, uint8_t can_write)
Returns a set of available GDAL drivers.
Definition: rt_raster.c:1716
rt_raster rt_raster_from_gdal_dataset(GDALDatasetH ds)
Return a raster from a GDAL dataset.
Definition: rt_raster.c:2182
static void _rti_rasterize_arg_destroy(_rti_rasterize_arg arg)
Definition: rt_raster.c:2459
double rt_raster_get_x_scale(rt_raster raster)
Get scale X in projection units.
Definition: rt_raster.c:150
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:1498
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:2509
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
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:372
rt_raster rt_raster_clone(rt_raster raster, uint8_t deep)
Clone an existing raster.
Definition: rt_raster.c:1540
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_raster.c:129
rt_band rt_raster_get_band(rt_raster raster, int n)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:381
static void _rt_raster_geotransform_warn_offline_band(rt_raster raster)
Definition: rt_raster.c:95
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster's SRID.
Definition: rt_raster.c:363
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 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:1435
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_raster.c:121
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:1826
struct _rti_rasterize_arg_t * _rti_rasterize_arg
Definition: rt_raster.c:2416
rt_errorstate rt_raster_from_two_rasters(rt_raster rast1, rt_raster rast2, rt_extenttype extenttype, rt_raster *rtnraster, double *offset)
Definition: rt_raster.c:3361
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:1370
double rt_raster_get_y_scale(rt_raster raster)
Get scale Y in projection units.
Definition: rt_raster.c:159
rt_errorstate rt_raster_get_envelope(rt_raster raster, rt_envelope *env)
Get raster's envelope.
Definition: rt_raster.c:873
static _rti_rasterize_arg _rti_rasterize_arg_init()
Definition: rt_raster.c:2433
double rt_raster_get_y_skew(rt_raster raster)
Get skew about the Y axis.
Definition: rt_raster.c:190
void rt_raster_get_geotransform_matrix(rt_raster raster, double *gt)
Get 6-element array of raster geotransform matrix.
Definition: rt_raster.c:706
void rt_raster_set_offsets(rt_raster raster, double x, double y)
Set insertion points in projection units.
Definition: rt_raster.c:199
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
double rt_raster_get_y_offset(rt_raster raster)
Get raster y offset, in projection units.
Definition: rt_raster.c:222
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
Definition: rt_raster.c:1334
rt_pixtype * pixtype
Definition: rt_raster.c:2424
OGRSpatialReferenceH src_sr
Definition: rt_raster.c:2422
rt_raster raster
Definition: librtcore.h:2324
double MinX
Definition: librtcore.h:165
double UpperLeftY
Definition: librtcore.h:171
double UpperLeftX
Definition: librtcore.h:170
double MaxX
Definition: librtcore.h:166
double MinY
Definition: librtcore.h:167
double MaxY
Definition: librtcore.h:168
uint8_t can_write
Definition: librtcore.h:2477
char * long_name
Definition: librtcore.h:2474
char * create_options
Definition: librtcore.h:2475
char * short_name
Definition: librtcore.h:2473
uint8_t can_read
Definition: librtcore.h:2476
double scaleX
Definition: librtcore.h:2293
rt_band * bands
Definition: librtcore.h:2303
uint16_t width
Definition: librtcore.h:2301
double skewY
Definition: librtcore.h:2298
uint16_t numBands
Definition: librtcore.h:2290
int32_t srid
Definition: librtcore.h:2300
double ipX
Definition: librtcore.h:2295
uint16_t height
Definition: librtcore.h:2302
double scaleY
Definition: librtcore.h:2294
double ipY
Definition: librtcore.h:2296
double skewX
Definition: librtcore.h:2297
unsigned int uint32_t
Definition: uthash.h:78
unsigned char uint8_t
Definition: uthash.h:79