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