PostGIS  2.4.9dev-r@@SVN_REVISION@@
rtpg_create.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) 2011-2013 Regents of the University of California
7  * <bkpark@ucdavis.edu>
8  * Copyright (C) 2010-2011 Jorge Arevalo <jorge.arevalo@deimos-space.com>
9  * Copyright (C) 2010-2011 David Zwarg <dzwarg@azavea.com>
10  * Copyright (C) 2009-2011 Pierre Racine <pierre.racine@sbf.ulaval.ca>
11  * Copyright (C) 2009-2011 Mateusz Loskot <mateusz@loskot.net>
12  * Copyright (C) 2008-2009 Sandro Santilli <strk@kbt.io>
13  *
14  * This program is free software; you can redistribute it and/or
15  * modify it under the terms of the GNU General Public License
16  * as published by the Free Software Foundation; either version 2
17  * of the License, or (at your option) any later version.
18  *
19  * This program is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU General Public License for more details.
23  *
24  * You should have received a copy of the GNU General Public License
25  * along with this program; if not, write to the Free Software Foundation,
26  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
27  *
28  */
29 
30 #include <postgres.h>
31 #include <fmgr.h>
32 #include <funcapi.h>
33 #include <utils/builtins.h> /* for text_to_cstring() */
34 #include "utils/lsyscache.h" /* for get_typlenbyvalalign */
35 #include "utils/array.h" /* for ArrayType */
36 #include "catalog/pg_type.h" /* for INT2OID, INT4OID, FLOAT4OID, FLOAT8OID and TEXTOID */
37 
38 #include "rtpostgis.h"
39 
40 /* Raster and band creation */
41 Datum RASTER_makeEmpty(PG_FUNCTION_ARGS);
42 Datum RASTER_addBand(PG_FUNCTION_ARGS);
43 Datum RASTER_addBandRasterArray(PG_FUNCTION_ARGS);
44 Datum RASTER_addBandOutDB(PG_FUNCTION_ARGS);
45 Datum RASTER_copyBand(PG_FUNCTION_ARGS);
46 Datum RASTER_tile(PG_FUNCTION_ARGS);
47 
48 /* create new raster from existing raster's bands */
49 Datum RASTER_band(PG_FUNCTION_ARGS);
50 
54 PG_FUNCTION_INFO_V1(RASTER_makeEmpty);
55 Datum RASTER_makeEmpty(PG_FUNCTION_ARGS)
56 {
57  uint16 width = 0, height = 0;
58  double ipx = 0, ipy = 0, scalex = 0, scaley = 0, skewx = 0, skewy = 0;
59  int32_t srid = SRID_UNKNOWN;
60  rt_pgraster *pgraster = NULL;
61  rt_raster raster;
62 
63  if (PG_NARGS() < 9) {
64  elog(ERROR, "RASTER_makeEmpty: ST_MakeEmptyRaster requires 9 args");
65  PG_RETURN_NULL();
66  }
67 
68  if (!PG_ARGISNULL(0))
69  width = PG_GETARG_UINT16(0);
70 
71  if (!PG_ARGISNULL(1))
72  height = PG_GETARG_UINT16(1);
73 
74  if (!PG_ARGISNULL(2))
75  ipx = PG_GETARG_FLOAT8(2);
76 
77  if (!PG_ARGISNULL(3))
78  ipy = PG_GETARG_FLOAT8(3);
79 
80  if (!PG_ARGISNULL(4))
81  scalex = PG_GETARG_FLOAT8(4);
82 
83  if (!PG_ARGISNULL(5))
84  scaley = PG_GETARG_FLOAT8(5);
85 
86  if (!PG_ARGISNULL(6))
87  skewx = PG_GETARG_FLOAT8(6);
88 
89  if (!PG_ARGISNULL(7))
90  skewy = PG_GETARG_FLOAT8(7);
91 
92  if (!PG_ARGISNULL(8))
93  srid = PG_GETARG_INT32(8);
94 
95  POSTGIS_RT_DEBUGF(4, "%dx%d, ip:%g,%g, scale:%g,%g, skew:%g,%g srid:%d",
96  width, height, ipx, ipy, scalex, scaley,
97  skewx, skewy, srid);
98 
99  raster = rt_raster_new(width, height);
100  if (raster == NULL)
101  PG_RETURN_NULL(); /* error was supposedly printed already */
102 
103  rt_raster_set_scale(raster, scalex, scaley);
104  rt_raster_set_offsets(raster, ipx, ipy);
105  rt_raster_set_skews(raster, skewx, skewy);
106  rt_raster_set_srid(raster, srid);
107 
108  pgraster = rt_raster_serialize(raster);
109  rt_raster_destroy(raster);
110  if (!pgraster)
111  PG_RETURN_NULL();
112 
113  SET_VARSIZE(pgraster, pgraster->size);
114  PG_RETURN_POINTER(pgraster);
115 }
116 
120 PG_FUNCTION_INFO_V1(RASTER_addBand);
121 Datum RASTER_addBand(PG_FUNCTION_ARGS)
122 {
123  rt_pgraster *pgraster = NULL;
124  rt_pgraster *pgrtn = NULL;
125  rt_raster raster = NULL;
126  int bandindex = 0;
127  int maxbandindex = 0;
128  int numbands = 0;
129  int lastnumbands = 0;
130 
131  text *text_pixtype = NULL;
132  char *char_pixtype = NULL;
133 
134  struct addbandarg {
135  int index;
136  bool append;
137  rt_pixtype pixtype;
138  double initialvalue;
139  bool hasnodata;
140  double nodatavalue;
141  };
142  struct addbandarg *arg = NULL;
143 
144  ArrayType *array;
145  Oid etype;
146  Datum *e;
147  bool *nulls;
148  int16 typlen;
149  bool typbyval;
150  char typalign;
151  int n = 0;
152 
153  HeapTupleHeader tup;
154  bool isnull;
155  Datum tupv;
156 
157  int i = 0;
158 
159  /* pgraster is null, return null */
160  if (PG_ARGISNULL(0))
161  PG_RETURN_NULL();
162  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
163 
164  /* raster */
165  raster = rt_raster_deserialize(pgraster, FALSE);
166  if (!raster) {
167  PG_FREE_IF_COPY(pgraster, 0);
168  elog(ERROR, "RASTER_addBand: Could not deserialize raster");
169  PG_RETURN_NULL();
170  }
171 
172  /* process set of addbandarg */
173  POSTGIS_RT_DEBUG(3, "Processing Arg 1 (addbandargset)");
174  array = PG_GETARG_ARRAYTYPE_P(1);
175  etype = ARR_ELEMTYPE(array);
176  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
177 
178  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
179  &nulls, &n);
180 
181  if (!n) {
182  PG_FREE_IF_COPY(pgraster, 0);
183  elog(ERROR, "RASTER_addBand: Invalid argument for addbandargset");
184  PG_RETURN_NULL();
185  }
186 
187  /* allocate addbandarg */
188  arg = (struct addbandarg *) palloc(sizeof(struct addbandarg) * n);
189  if (arg == NULL) {
190  rt_raster_destroy(raster);
191  PG_FREE_IF_COPY(pgraster, 0);
192  elog(ERROR, "RASTER_addBand: Could not allocate memory for addbandarg");
193  PG_RETURN_NULL();
194  }
195 
196  /*
197  process each element of addbandargset
198  each element is the index of where to add the new band,
199  new band's pixeltype, the new band's initial value and
200  the new band's NODATA value if NOT NULL
201  */
202  for (i = 0; i < n; i++) {
203  if (nulls[i]) continue;
204 
205  POSTGIS_RT_DEBUGF(4, "Processing addbandarg at index %d", i);
206 
207  /* each element is a tuple */
208  tup = (HeapTupleHeader) DatumGetPointer(e[i]);
209  if (NULL == tup) {
210  pfree(arg);
211  rt_raster_destroy(raster);
212  PG_FREE_IF_COPY(pgraster, 0);
213  elog(ERROR, "RASTER_addBand: Invalid argument for addbandargset");
214  PG_RETURN_NULL();
215  }
216 
217  /* new band index, 1-based */
218  arg[i].index = 0;
219  arg[i].append = TRUE;
220  tupv = GetAttributeByName(tup, "index", &isnull);
221  if (!isnull) {
222  arg[i].index = DatumGetInt32(tupv);
223  arg[i].append = FALSE;
224  }
225 
226  /* for now, only check that band index is 1-based */
227  if (!arg[i].append && arg[i].index < 1) {
228  pfree(arg);
229  rt_raster_destroy(raster);
230  PG_FREE_IF_COPY(pgraster, 0);
231  elog(ERROR, "RASTER_addBand: Invalid argument for addbandargset. Invalid band index (must be 1-based) for addbandarg of index %d", i);
232  PG_RETURN_NULL();
233  }
234 
235  /* new band pixeltype */
236  arg[i].pixtype = PT_END;
237  tupv = GetAttributeByName(tup, "pixeltype", &isnull);
238  if (isnull) {
239  pfree(arg);
240  rt_raster_destroy(raster);
241  PG_FREE_IF_COPY(pgraster, 0);
242  elog(ERROR, "RASTER_addBand: Invalid argument for addbandargset. Pixel type cannot be NULL for addbandarg of index %d", i);
243  PG_RETURN_NULL();
244  }
245  text_pixtype = (text *) DatumGetPointer(tupv);
246  if (text_pixtype == NULL) {
247  pfree(arg);
248  rt_raster_destroy(raster);
249  PG_FREE_IF_COPY(pgraster, 0);
250  elog(ERROR, "RASTER_addBand: Invalid argument for addbandargset. Pixel type cannot be NULL for addbandarg of index %d", i);
251  PG_RETURN_NULL();
252  }
253  char_pixtype = text_to_cstring(text_pixtype);
254 
255  arg[i].pixtype = rt_pixtype_index_from_name(char_pixtype);
256  pfree(char_pixtype);
257  if (arg[i].pixtype == PT_END) {
258  pfree(arg);
259  rt_raster_destroy(raster);
260  PG_FREE_IF_COPY(pgraster, 0);
261  elog(ERROR, "RASTER_addBand: Invalid argument for addbandargset. Invalid pixel type for addbandarg of index %d", i);
262  PG_RETURN_NULL();
263  }
264 
265  /* new band initialvalue */
266  arg[i].initialvalue = 0;
267  tupv = GetAttributeByName(tup, "initialvalue", &isnull);
268  if (!isnull)
269  arg[i].initialvalue = DatumGetFloat8(tupv);
270 
271  /* new band NODATA value */
272  arg[i].hasnodata = FALSE;
273  arg[i].nodatavalue = 0;
274  tupv = GetAttributeByName(tup, "nodataval", &isnull);
275  if (!isnull) {
276  arg[i].hasnodata = TRUE;
277  arg[i].nodatavalue = DatumGetFloat8(tupv);
278  }
279  }
280 
281  /* add new bands to raster */
282  lastnumbands = rt_raster_get_num_bands(raster);
283  for (i = 0; i < n; i++) {
284  if (nulls[i]) continue;
285 
286  POSTGIS_RT_DEBUGF(3, "%d bands in old raster", lastnumbands);
287  maxbandindex = lastnumbands + 1;
288 
289  /* check that new band's index doesn't exceed maxbandindex */
290  if (!arg[i].append) {
291  if (arg[i].index > maxbandindex) {
292  elog(NOTICE, "Band index for addbandarg of index %d exceeds possible value. Adding band at index %d", i, maxbandindex);
293  arg[i].index = maxbandindex;
294  }
295  }
296  /* append, so use maxbandindex */
297  else
298  arg[i].index = maxbandindex;
299 
300  POSTGIS_RT_DEBUGF(4, "new band (index, pixtype, initialvalue, hasnodata, nodatavalue) = (%d, %s, %f, %s, %f)",
301  arg[i].index,
302  rt_pixtype_name(arg[i].pixtype),
303  arg[i].initialvalue,
304  arg[i].hasnodata ? "TRUE" : "FALSE",
305  arg[i].nodatavalue
306  );
307 
308  bandindex = rt_raster_generate_new_band(
309  raster,
310  arg[i].pixtype, arg[i].initialvalue,
311  arg[i].hasnodata, arg[i].nodatavalue,
312  arg[i].index - 1
313  );
314 
315  numbands = rt_raster_get_num_bands(raster);
316  if (numbands == lastnumbands || bandindex == -1) {
317  pfree(arg);
318  rt_raster_destroy(raster);
319  PG_FREE_IF_COPY(pgraster, 0);
320  elog(ERROR, "RASTER_addBand: Could not add band defined by addbandarg of index %d to raster", i);
321  PG_RETURN_NULL();
322  }
323 
324  lastnumbands = numbands;
325  POSTGIS_RT_DEBUGF(3, "%d bands in new raster", lastnumbands);
326  }
327 
328  pfree(arg);
329 
330  pgrtn = rt_raster_serialize(raster);
331  rt_raster_destroy(raster);
332  PG_FREE_IF_COPY(pgraster, 0);
333  if (!pgrtn)
334  PG_RETURN_NULL();
335 
336  SET_VARSIZE(pgrtn, pgrtn->size);
337  PG_RETURN_POINTER(pgrtn);
338 }
339 
343 PG_FUNCTION_INFO_V1(RASTER_addBandRasterArray);
344 Datum RASTER_addBandRasterArray(PG_FUNCTION_ARGS)
345 {
346  rt_pgraster *pgraster = NULL;
347  rt_pgraster *pgsrc = NULL;
348  rt_pgraster *pgrtn = NULL;
349 
350  rt_raster raster = NULL;
351  rt_raster src = NULL;
352 
353  int srcnband = 1;
354  bool appendband = FALSE;
355  int dstnband = 1;
356  int srcnumbands = 0;
357  int dstnumbands = 0;
358 
359  ArrayType *array;
360  Oid etype;
361  Datum *e;
362  bool *nulls;
363  int16 typlen;
364  bool typbyval;
365  char typalign;
366  int n = 0;
367 
368  int rtn = 0;
369  int i = 0;
370 
371  /* destination raster */
372  if (!PG_ARGISNULL(0)) {
373  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
374 
375  /* raster */
376  raster = rt_raster_deserialize(pgraster, FALSE);
377  if (!raster) {
378  PG_FREE_IF_COPY(pgraster, 0);
379  elog(ERROR, "RASTER_addBandRasterArray: Could not deserialize destination raster");
380  PG_RETURN_NULL();
381  }
382 
383  POSTGIS_RT_DEBUG(4, "destination raster isn't NULL");
384  }
385 
386  /* source rasters' band index, 1-based */
387  if (!PG_ARGISNULL(2))
388  srcnband = PG_GETARG_INT32(2);
389  if (srcnband < 1) {
390  elog(NOTICE, "Invalid band index for source rasters (must be 1-based). Returning original raster");
391  if (raster != NULL) {
392  rt_raster_destroy(raster);
393  PG_RETURN_POINTER(pgraster);
394  }
395  else
396  PG_RETURN_NULL();
397  }
398  POSTGIS_RT_DEBUGF(4, "srcnband = %d", srcnband);
399 
400  /* destination raster's band index, 1-based */
401  if (!PG_ARGISNULL(3)) {
402  dstnband = PG_GETARG_INT32(3);
403  appendband = FALSE;
404 
405  if (dstnband < 1) {
406  elog(NOTICE, "Invalid band index for destination raster (must be 1-based). Returning original raster");
407  if (raster != NULL) {
408  rt_raster_destroy(raster);
409  PG_RETURN_POINTER(pgraster);
410  }
411  else
412  PG_RETURN_NULL();
413  }
414  }
415  else
416  appendband = TRUE;
417 
418  /* additional processing of dstnband */
419  if (raster != NULL) {
420  dstnumbands = rt_raster_get_num_bands(raster);
421 
422  if (dstnumbands < 1) {
423  appendband = TRUE;
424  dstnband = 1;
425  }
426  else if (appendband)
427  dstnband = dstnumbands + 1;
428  else if (dstnband > dstnumbands) {
429  elog(NOTICE, "Band index provided for destination raster is greater than the number of bands in the raster. Bands will be appended");
430  appendband = TRUE;
431  dstnband = dstnumbands + 1;
432  }
433  }
434  POSTGIS_RT_DEBUGF(4, "appendband = %d", appendband);
435  POSTGIS_RT_DEBUGF(4, "dstnband = %d", dstnband);
436 
437  /* process set of source rasters */
438  POSTGIS_RT_DEBUG(3, "Processing array of source rasters");
439  array = PG_GETARG_ARRAYTYPE_P(1);
440  etype = ARR_ELEMTYPE(array);
441  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
442 
443  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
444  &nulls, &n);
445 
446  /* decrement srcnband and dstnband by 1, now 0-based */
447  srcnband--;
448  dstnband--;
449  POSTGIS_RT_DEBUGF(4, "0-based nband (src, dst) = (%d, %d)", srcnband, dstnband);
450 
451  /* time to copy bands */
452  for (i = 0; i < n; i++) {
453  if (nulls[i]) continue;
454  src = NULL;
455 
456  pgsrc = (rt_pgraster *) PG_DETOAST_DATUM(e[i]);
457  src = rt_raster_deserialize(pgsrc, FALSE);
458  if (src == NULL) {
459  pfree(nulls);
460  pfree(e);
461  if (raster != NULL)
462  rt_raster_destroy(raster);
463  if (pgraster != NULL)
464  PG_FREE_IF_COPY(pgraster, 0);
465  elog(ERROR, "RASTER_addBandRasterArray: Could not deserialize source raster at index %d", i + 1);
466  PG_RETURN_NULL();
467  }
468 
469  srcnumbands = rt_raster_get_num_bands(src);
470  POSTGIS_RT_DEBUGF(4, "source raster %d has %d bands", i + 1, srcnumbands);
471 
472  /* band index isn't valid */
473  if (srcnband > srcnumbands - 1) {
474  elog(NOTICE, "Invalid band index for source raster at index %d. Returning original raster", i + 1);
475  pfree(nulls);
476  pfree(e);
477  rt_raster_destroy(src);
478  if (raster != NULL) {
479  rt_raster_destroy(raster);
480  PG_RETURN_POINTER(pgraster);
481  }
482  else
483  PG_RETURN_NULL();
484  }
485 
486  /* destination raster is empty, new raster */
487  if (raster == NULL) {
488  uint32_t srcnbands[1] = {srcnband};
489 
490  POSTGIS_RT_DEBUG(4, "empty destination raster, using rt_raster_from_band");
491 
492  raster = rt_raster_from_band(src, srcnbands, 1);
493  rt_raster_destroy(src);
494  if (raster == NULL) {
495  pfree(nulls);
496  pfree(e);
497  if (pgraster != NULL)
498  PG_FREE_IF_COPY(pgraster, 0);
499  elog(ERROR, "RASTER_addBandRasterArray: Could not create raster from source raster at index %d", i + 1);
500  PG_RETURN_NULL();
501  }
502  }
503  /* copy band */
504  else {
505  rtn = rt_raster_copy_band(
506  raster, src,
507  srcnband, dstnband
508  );
509  rt_raster_destroy(src);
510 
511  if (rtn == -1 || rt_raster_get_num_bands(raster) == dstnumbands) {
512  elog(NOTICE, "Could not add band from source raster at index %d to destination raster. Returning original raster", i + 1);
513  rt_raster_destroy(raster);
514  pfree(nulls);
515  pfree(e);
516  if (pgraster != NULL)
517  PG_RETURN_POINTER(pgraster);
518  else
519  PG_RETURN_NULL();
520  }
521  }
522 
523  dstnband++;
524  dstnumbands++;
525  }
526 
527  if (raster != NULL) {
528  pgrtn = rt_raster_serialize(raster);
529  rt_raster_destroy(raster);
530  if (pgraster != NULL)
531  PG_FREE_IF_COPY(pgraster, 0);
532  if (!pgrtn)
533  PG_RETURN_NULL();
534 
535  SET_VARSIZE(pgrtn, pgrtn->size);
536  PG_RETURN_POINTER(pgrtn);
537  }
538 
539  PG_RETURN_NULL();
540 }
541 
545 PG_FUNCTION_INFO_V1(RASTER_addBandOutDB);
546 Datum RASTER_addBandOutDB(PG_FUNCTION_ARGS)
547 {
548  rt_pgraster *pgraster = NULL;
549  rt_pgraster *pgrtn = NULL;
550 
551  rt_raster raster = NULL;
552  rt_band band = NULL;
553  int numbands = 0;
554  int dstnband = 1; /* 1-based */
555  int appendband = FALSE;
556  char *outdbfile = NULL;
557  int *srcnband = NULL; /* 1-based */
558  int numsrcnband = 0;
559  int allbands = FALSE;
560  int hasnodata = FALSE;
561  double nodataval = 0.;
562  uint16_t width = 0;
563  uint16_t height = 0;
564  char *authname = NULL;
565  char *authcode = NULL;
566 
567  int i = 0;
568  int j = 0;
569 
570  GDALDatasetH hdsOut;
571  GDALRasterBandH hbandOut;
572  GDALDataType gdpixtype;
573 
574  rt_pixtype pt = PT_END;
575  double gt[6] = {0.};
576  double ogt[6] = {0.};
577  rt_raster _rast = NULL;
578  int aligned = 0;
579  int err = 0;
580 
581  /* destination raster */
582  if (!PG_ARGISNULL(0)) {
583  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
584 
585  /* raster */
586  raster = rt_raster_deserialize(pgraster, FALSE);
587  if (!raster) {
588  PG_FREE_IF_COPY(pgraster, 0);
589  elog(ERROR, "RASTER_addBandOutDB: Could not deserialize destination raster");
590  PG_RETURN_NULL();
591  }
592 
593  POSTGIS_RT_DEBUG(4, "destination raster isn't NULL");
594  }
595 
596  /* destination band index (1) */
597  if (!PG_ARGISNULL(1))
598  dstnband = PG_GETARG_INT32(1);
599  else
600  appendband = TRUE;
601 
602  /* outdb file (2) */
603  if (PG_ARGISNULL(2)) {
604  elog(NOTICE, "Out-db raster file not provided. Returning original raster");
605  if (pgraster != NULL) {
606  rt_raster_destroy(raster);
607  PG_RETURN_POINTER(pgraster);
608  }
609  else
610  PG_RETURN_NULL();
611  }
612  else {
613  outdbfile = text_to_cstring(PG_GETARG_TEXT_P(2));
614  if (!strlen(outdbfile)) {
615  elog(NOTICE, "Out-db raster file not provided. Returning original raster");
616  if (pgraster != NULL) {
617  rt_raster_destroy(raster);
618  PG_RETURN_POINTER(pgraster);
619  }
620  else
621  PG_RETURN_NULL();
622  }
623  }
624 
625  /* outdb band index (3) */
626  if (!PG_ARGISNULL(3)) {
627  ArrayType *array;
628  Oid etype;
629  Datum *e;
630  bool *nulls;
631 
632  int16 typlen;
633  bool typbyval;
634  char typalign;
635 
636  allbands = FALSE;
637 
638  array = PG_GETARG_ARRAYTYPE_P(3);
639  etype = ARR_ELEMTYPE(array);
640  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
641 
642  switch (etype) {
643  case INT2OID:
644  case INT4OID:
645  break;
646  default:
647  if (pgraster != NULL) {
648  rt_raster_destroy(raster);
649  PG_FREE_IF_COPY(pgraster, 0);
650  }
651  elog(ERROR, "RASTER_addBandOutDB: Invalid data type for band indexes");
652  PG_RETURN_NULL();
653  break;
654  }
655 
656  deconstruct_array(array, etype, typlen, typbyval, typalign, &e, &nulls, &numsrcnband);
657 
658  srcnband = palloc(sizeof(int) * numsrcnband);
659  if (srcnband == NULL) {
660  if (pgraster != NULL) {
661  rt_raster_destroy(raster);
662  PG_FREE_IF_COPY(pgraster, 0);
663  }
664  elog(ERROR, "RASTER_addBandOutDB: Could not allocate memory for band indexes");
665  PG_RETURN_NULL();
666  }
667 
668  for (i = 0, j = 0; i < numsrcnband; i++) {
669  if (nulls[i]) continue;
670 
671  switch (etype) {
672  case INT2OID:
673  srcnband[j] = DatumGetInt16(e[i]);
674  break;
675  case INT4OID:
676  srcnband[j] = DatumGetInt32(e[i]);
677  break;
678  }
679  j++;
680  }
681 
682  if (j < numsrcnband) {
683  srcnband = repalloc(srcnband, sizeof(int) * j);
684  if (srcnband == NULL) {
685  if (pgraster != NULL) {
686  rt_raster_destroy(raster);
687  PG_FREE_IF_COPY(pgraster, 0);
688  }
689  elog(ERROR, "RASTER_addBandOutDB: Could not reallocate memory for band indexes");
690  PG_RETURN_NULL();
691  }
692 
693  numsrcnband = j;
694  }
695  }
696  else
697  allbands = TRUE;
698 
699  /* nodataval (4) */
700  if (!PG_ARGISNULL(4)) {
701  hasnodata = TRUE;
702  nodataval = PG_GETARG_FLOAT8(4);
703  }
704  else
705  hasnodata = FALSE;
706 
707  /* validate input */
708 
709  /* make sure dstnband is valid */
710  if (raster != NULL) {
711  numbands = rt_raster_get_num_bands(raster);
712  if (!appendband) {
713  if (dstnband < 1) {
714  elog(NOTICE, "Invalid band index %d for adding bands. Using band index 1", dstnband);
715  dstnband = 1;
716  }
717  else if (numbands > 0 && dstnband > numbands) {
718  elog(NOTICE, "Invalid band index %d for adding bands. Using band index %d", dstnband, numbands);
719  dstnband = numbands + 1;
720  }
721  }
722  else
723  dstnband = numbands + 1;
724  }
725 
726  /* open outdb raster file */
727  rt_util_gdal_register_all(0);
728  hdsOut = rt_util_gdal_open(outdbfile, GA_ReadOnly, 0);
729  if (hdsOut == NULL) {
730  if (pgraster != NULL) {
731  rt_raster_destroy(raster);
732  PG_FREE_IF_COPY(pgraster, 0);
733  }
734  elog(ERROR, "RASTER_addBandOutDB: Could not open out-db file with GDAL");
735  PG_RETURN_NULL();
736  }
737 
738  /* get offline raster's geotransform */
739  if (GDALGetGeoTransform(hdsOut, ogt) != CE_None) {
740  ogt[0] = 0;
741  ogt[1] = 1;
742  ogt[2] = 0;
743  ogt[3] = 0;
744  ogt[4] = 0;
745  ogt[5] = -1;
746  }
747 
748  /* raster doesn't exist, create it now */
749  if (raster == NULL) {
750  raster = rt_raster_new(GDALGetRasterXSize(hdsOut), GDALGetRasterYSize(hdsOut));
751  if (rt_raster_is_empty(raster)) {
752  elog(ERROR, "RASTER_addBandOutDB: Could not create new raster");
753  PG_RETURN_NULL();
754  }
755  rt_raster_set_geotransform_matrix(raster, ogt);
756  rt_raster_get_geotransform_matrix(raster, gt);
757 
758  if (rt_util_gdal_sr_auth_info(hdsOut, &authname, &authcode) == ES_NONE) {
759  if (
760  authname != NULL &&
761  strcmp(authname, "EPSG") == 0 &&
762  authcode != NULL
763  ) {
764  rt_raster_set_srid(raster, atoi(authcode));
765  }
766  else
767  elog(INFO, "Unknown SRS auth name and code from out-db file. Defaulting SRID of new raster to %d", SRID_UNKNOWN);
768  }
769  else
770  elog(INFO, "Could not get SRS auth name and code from out-db file. Defaulting SRID of new raster to %d", SRID_UNKNOWN);
771  }
772 
773  /* some raster info */
774  width = rt_raster_get_width(raster);
775  height = rt_raster_get_height(raster);
776 
777  /* are rasters aligned? */
778  _rast = rt_raster_new(1, 1);
779  rt_raster_set_geotransform_matrix(_rast, ogt);
780  rt_raster_set_srid(_rast, rt_raster_get_srid(raster));
781  err = rt_raster_same_alignment(raster, _rast, &aligned, NULL);
782  rt_raster_destroy(_rast);
783 
784  if (err != ES_NONE) {
785  GDALClose(hdsOut);
786  if (raster != NULL)
787  rt_raster_destroy(raster);
788  if (pgraster != NULL)
789  PG_FREE_IF_COPY(pgraster, 0);
790  elog(ERROR, "RASTER_addBandOutDB: Could not test alignment of out-db file");
791  return ES_ERROR;
792  }
793  else if (!aligned)
794  elog(WARNING, "The in-db representation of the out-db raster is not aligned. Band data may be incorrect");
795 
796  numbands = GDALGetRasterCount(hdsOut);
797 
798  /* build up srcnband */
799  if (allbands) {
800  numsrcnband = numbands;
801  srcnband = palloc(sizeof(int) * numsrcnband);
802  if (srcnband == NULL) {
803  GDALClose(hdsOut);
804  if (raster != NULL)
805  rt_raster_destroy(raster);
806  if (pgraster != NULL)
807  PG_FREE_IF_COPY(pgraster, 0);
808  elog(ERROR, "RASTER_addBandOutDB: Could not allocate memory for band indexes");
809  PG_RETURN_NULL();
810  }
811 
812  for (i = 0, j = 1; i < numsrcnband; i++, j++)
813  srcnband[i] = j;
814  }
815 
816  /* check band properties and add band */
817  for (i = 0, j = dstnband - 1; i < numsrcnband; i++, j++) {
818  /* valid index? */
819  if (srcnband[i] < 1 || srcnband[i] > numbands) {
820  elog(NOTICE, "Out-db file does not have a band at index %d. Returning original raster", srcnband[i]);
821  GDALClose(hdsOut);
822  if (raster != NULL)
823  rt_raster_destroy(raster);
824  if (pgraster != NULL)
825  PG_RETURN_POINTER(pgraster);
826  else
827  PG_RETURN_NULL();
828  }
829 
830  /* get outdb band */
831  hbandOut = NULL;
832  hbandOut = GDALGetRasterBand(hdsOut, srcnband[i]);
833  if (NULL == hbandOut) {
834  GDALClose(hdsOut);
835  if (raster != NULL)
836  rt_raster_destroy(raster);
837  if (pgraster != NULL)
838  PG_FREE_IF_COPY(pgraster, 0);
839  elog(ERROR, "RASTER_addBandOutDB: Could not get band %d from GDAL dataset", srcnband[i]);
840  PG_RETURN_NULL();
841  }
842 
843  /* supported pixel type */
844  gdpixtype = GDALGetRasterDataType(hbandOut);
845  pt = rt_util_gdal_datatype_to_pixtype(gdpixtype);
846  if (pt == PT_END) {
847  elog(NOTICE, "Pixel type %s of band %d from GDAL dataset is not supported. Returning original raster", GDALGetDataTypeName(gdpixtype), srcnband[i]);
848  GDALClose(hdsOut);
849  if (raster != NULL)
850  rt_raster_destroy(raster);
851  if (pgraster != NULL)
852  PG_RETURN_POINTER(pgraster);
853  else
854  PG_RETURN_NULL();
855  }
856 
857  /* use out-db band's nodata value if nodataval not already set */
858  if (!hasnodata)
859  nodataval = GDALGetRasterNoDataValue(hbandOut, &hasnodata);
860 
861  /* add band */
862  band = rt_band_new_offline(
863  width, height,
864  pt,
865  hasnodata, nodataval,
866  srcnband[i] - 1, outdbfile
867  );
868  if (band == NULL) {
869  GDALClose(hdsOut);
870  if (raster != NULL)
871  rt_raster_destroy(raster);
872  if (pgraster != NULL)
873  PG_FREE_IF_COPY(pgraster, 0);
874  elog(ERROR, "RASTER_addBandOutDB: Could not create new out-db band");
875  PG_RETURN_NULL();
876  }
877 
878  if (rt_raster_add_band(raster, band, j) < 0) {
879  GDALClose(hdsOut);
880  if (raster != NULL)
881  rt_raster_destroy(raster);
882  if (pgraster != NULL)
883  PG_FREE_IF_COPY(pgraster, 0);
884  elog(ERROR, "RASTER_addBandOutDB: Could not add new out-db band to raster");
885  PG_RETURN_NULL();
886  }
887  }
888 
889  pgrtn = rt_raster_serialize(raster);
890  rt_raster_destroy(raster);
891  if (pgraster != NULL)
892  PG_FREE_IF_COPY(pgraster, 0);
893  if (!pgrtn)
894  PG_RETURN_NULL();
895 
896  SET_VARSIZE(pgrtn, pgrtn->size);
897  PG_RETURN_POINTER(pgrtn);
898 }
899 
903 PG_FUNCTION_INFO_V1(RASTER_copyBand);
904 Datum RASTER_copyBand(PG_FUNCTION_ARGS)
905 {
906  rt_pgraster *pgto = NULL;
907  rt_pgraster *pgfrom = NULL;
908  rt_pgraster *pgrtn = NULL;
909  rt_raster torast = NULL;
910  rt_raster fromrast = NULL;
911  int toindex = 0;
912  int fromband = 0;
913  int oldtorastnumbands = 0;
914  int newtorastnumbands = 0;
915  int newbandindex = 0;
916 
917  /* Deserialize torast */
918  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
919  pgto = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
920 
921  torast = rt_raster_deserialize(pgto, FALSE);
922  if (!torast) {
923  PG_FREE_IF_COPY(pgto, 0);
924  elog(ERROR, "RASTER_copyBand: Could not deserialize first raster");
925  PG_RETURN_NULL();
926  }
927 
928  /* Deserialize fromrast */
929  if (!PG_ARGISNULL(1)) {
930  pgfrom = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
931 
932  fromrast = rt_raster_deserialize(pgfrom, FALSE);
933  if (!fromrast) {
934  rt_raster_destroy(torast);
935  PG_FREE_IF_COPY(pgfrom, 1);
936  PG_FREE_IF_COPY(pgto, 0);
937  elog(ERROR, "RASTER_copyBand: Could not deserialize second raster");
938  PG_RETURN_NULL();
939  }
940 
941  oldtorastnumbands = rt_raster_get_num_bands(torast);
942 
943  if (PG_ARGISNULL(2))
944  fromband = 1;
945  else
946  fromband = PG_GETARG_INT32(2);
947 
948  if (PG_ARGISNULL(3))
949  toindex = oldtorastnumbands + 1;
950  else
951  toindex = PG_GETARG_INT32(3);
952 
953  /* Copy band fromrast torast */
954  newbandindex = rt_raster_copy_band(
955  torast, fromrast,
956  fromband - 1, toindex - 1
957  );
958 
959  newtorastnumbands = rt_raster_get_num_bands(torast);
960  if (newtorastnumbands == oldtorastnumbands || newbandindex == -1) {
961  elog(NOTICE, "RASTER_copyBand: Could not add band to raster. "
962  "Returning original raster."
963  );
964  }
965 
966  rt_raster_destroy(fromrast);
967  PG_FREE_IF_COPY(pgfrom, 1);
968  }
969 
970  /* Serialize and return torast */
971  pgrtn = rt_raster_serialize(torast);
972  rt_raster_destroy(torast);
973  PG_FREE_IF_COPY(pgto, 0);
974  if (!pgrtn) PG_RETURN_NULL();
975 
976  SET_VARSIZE(pgrtn, pgrtn->size);
977  PG_RETURN_POINTER(pgrtn);
978 }
979 
983 PG_FUNCTION_INFO_V1(RASTER_tile);
984 Datum RASTER_tile(PG_FUNCTION_ARGS)
985 {
986  FuncCallContext *funcctx;
987  int call_cntr;
988  int max_calls;
989  int i = 0;
990  int j = 0;
991 
992  struct tile_arg_t {
993 
994  struct {
995  rt_raster raster;
996  double gt[6];
997  int srid;
998  int width;
999  int height;
1000  } raster;
1001 
1002  struct {
1003  int width;
1004  int height;
1005 
1006  int nx;
1007  int ny;
1008  } tile;
1009 
1010  int numbands;
1011  int *nbands;
1012 
1013  struct {
1014  int pad;
1015  double hasnodata;
1016  double nodataval;
1017  } pad;
1018  };
1019  struct tile_arg_t *arg1 = NULL;
1020  struct tile_arg_t *arg2 = NULL;
1021 
1022  if (SRF_IS_FIRSTCALL()) {
1023  MemoryContext oldcontext;
1024  rt_pgraster *pgraster = NULL;
1025  int numbands;
1026 
1027  ArrayType *array;
1028  Oid etype;
1029  Datum *e;
1030  bool *nulls;
1031 
1032  int16 typlen;
1033  bool typbyval;
1034  char typalign;
1035 
1036  POSTGIS_RT_DEBUG(2, "RASTER_tile: first call");
1037 
1038  /* create a function context for cross-call persistence */
1039  funcctx = SRF_FIRSTCALL_INIT();
1040 
1041  /* switch to memory context appropriate for multiple function calls */
1042  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
1043 
1044  /* Get input arguments */
1045  if (PG_ARGISNULL(0)) {
1046  MemoryContextSwitchTo(oldcontext);
1047  SRF_RETURN_DONE(funcctx);
1048  }
1049 
1050  /* allocate arg1 */
1051  arg1 = palloc(sizeof(struct tile_arg_t));
1052  if (arg1 == NULL) {
1053  MemoryContextSwitchTo(oldcontext);
1054  elog(ERROR, "RASTER_tile: Could not allocate memory for arguments");
1055  SRF_RETURN_DONE(funcctx);
1056  }
1057 
1058  pgraster = (rt_pgraster *) PG_DETOAST_DATUM_COPY(PG_GETARG_DATUM(0));
1059  arg1->raster.raster = rt_raster_deserialize(pgraster, FALSE);
1060  if (!arg1->raster.raster) {
1061  ereport(ERROR, (
1062  errcode(ERRCODE_OUT_OF_MEMORY),
1063  errmsg("Could not deserialize raster")
1064  ));
1065  pfree(arg1);
1066  PG_FREE_IF_COPY(pgraster, 0);
1067  MemoryContextSwitchTo(oldcontext);
1068  SRF_RETURN_DONE(funcctx);
1069  }
1070 
1071  /* raster has bands */
1072  numbands = rt_raster_get_num_bands(arg1->raster.raster);
1073  /*
1074  if (!numbands) {
1075  elog(NOTICE, "Raster provided has no bands");
1076  rt_raster_destroy(arg1->raster.raster);
1077  pfree(arg1);
1078  PG_FREE_IF_COPY(pgraster, 0);
1079  MemoryContextSwitchTo(oldcontext);
1080  SRF_RETURN_DONE(funcctx);
1081  }
1082  */
1083 
1084  /* width (1) */
1085  if (PG_ARGISNULL(1)) {
1086  elog(NOTICE, "Width cannot be NULL. Returning NULL");
1087  rt_raster_destroy(arg1->raster.raster);
1088  pfree(arg1);
1089  PG_FREE_IF_COPY(pgraster, 0);
1090  MemoryContextSwitchTo(oldcontext);
1091  SRF_RETURN_DONE(funcctx);
1092  }
1093  arg1->tile.width = PG_GETARG_INT32(1);
1094  if (arg1->tile.width < 1) {
1095  elog(NOTICE, "Width must be greater than zero. Returning NULL");
1096  rt_raster_destroy(arg1->raster.raster);
1097  pfree(arg1);
1098  PG_FREE_IF_COPY(pgraster, 0);
1099  MemoryContextSwitchTo(oldcontext);
1100  SRF_RETURN_DONE(funcctx);
1101  }
1102 
1103  /* height (2) */
1104  if (PG_ARGISNULL(2)) {
1105  elog(NOTICE, "Height cannot be NULL. Returning NULL");
1106  rt_raster_destroy(arg1->raster.raster);
1107  pfree(arg1);
1108  PG_FREE_IF_COPY(pgraster, 0);
1109  MemoryContextSwitchTo(oldcontext);
1110  SRF_RETURN_DONE(funcctx);
1111  }
1112  arg1->tile.height = PG_GETARG_INT32(2);
1113  if (arg1->tile.height < 1) {
1114  elog(NOTICE, "Height must be greater than zero. Returning NULL");
1115  rt_raster_destroy(arg1->raster.raster);
1116  pfree(arg1);
1117  PG_FREE_IF_COPY(pgraster, 0);
1118  MemoryContextSwitchTo(oldcontext);
1119  SRF_RETURN_DONE(funcctx);
1120  }
1121 
1122  /* nband, array (3) */
1123  if (numbands && !PG_ARGISNULL(3)) {
1124  array = PG_GETARG_ARRAYTYPE_P(3);
1125  etype = ARR_ELEMTYPE(array);
1126  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1127 
1128  switch (etype) {
1129  case INT2OID:
1130  case INT4OID:
1131  break;
1132  default:
1133  rt_raster_destroy(arg1->raster.raster);
1134  pfree(arg1);
1135  PG_FREE_IF_COPY(pgraster, 0);
1136  MemoryContextSwitchTo(oldcontext);
1137  elog(ERROR, "RASTER_tile: Invalid data type for band indexes");
1138  SRF_RETURN_DONE(funcctx);
1139  break;
1140  }
1141 
1142  deconstruct_array(array, etype, typlen, typbyval, typalign, &e, &nulls, &(arg1->numbands));
1143 
1144  arg1->nbands = palloc(sizeof(int) * arg1->numbands);
1145  if (arg1->nbands == NULL) {
1146  rt_raster_destroy(arg1->raster.raster);
1147  pfree(arg1);
1148  PG_FREE_IF_COPY(pgraster, 0);
1149  MemoryContextSwitchTo(oldcontext);
1150  elog(ERROR, "RASTER_tile: Could not allocate memory for band indexes");
1151  SRF_RETURN_DONE(funcctx);
1152  }
1153 
1154  for (i = 0, j = 0; i < arg1->numbands; i++) {
1155  if (nulls[i]) continue;
1156 
1157  switch (etype) {
1158  case INT2OID:
1159  arg1->nbands[j] = DatumGetInt16(e[i]) - 1;
1160  break;
1161  case INT4OID:
1162  arg1->nbands[j] = DatumGetInt32(e[i]) - 1;
1163  break;
1164  }
1165 
1166  j++;
1167  }
1168 
1169  if (j < arg1->numbands) {
1170  arg1->nbands = repalloc(arg1->nbands, sizeof(int) * j);
1171  if (arg1->nbands == NULL) {
1172  rt_raster_destroy(arg1->raster.raster);
1173  pfree(arg1);
1174  PG_FREE_IF_COPY(pgraster, 0);
1175  MemoryContextSwitchTo(oldcontext);
1176  elog(ERROR, "RASTER_tile: Could not reallocate memory for band indexes");
1177  SRF_RETURN_DONE(funcctx);
1178  }
1179 
1180  arg1->numbands = j;
1181  }
1182 
1183  /* validate nbands */
1184  for (i = 0; i < arg1->numbands; i++) {
1185  if (!rt_raster_has_band(arg1->raster.raster, arg1->nbands[i])) {
1186  elog(NOTICE, "Band at index %d not found in raster", arg1->nbands[i] + 1);
1187  rt_raster_destroy(arg1->raster.raster);
1188  pfree(arg1->nbands);
1189  pfree(arg1);
1190  PG_FREE_IF_COPY(pgraster, 0);
1191  MemoryContextSwitchTo(oldcontext);
1192  SRF_RETURN_DONE(funcctx);
1193  }
1194  }
1195  }
1196  else {
1197  arg1->numbands = numbands;
1198 
1199  if (numbands) {
1200  arg1->nbands = palloc(sizeof(int) * arg1->numbands);
1201 
1202  if (arg1->nbands == NULL) {
1203  rt_raster_destroy(arg1->raster.raster);
1204  pfree(arg1);
1205  PG_FREE_IF_COPY(pgraster, 0);
1206  MemoryContextSwitchTo(oldcontext);
1207  elog(ERROR, "RASTER_dumpValues: Could not allocate memory for pixel values");
1208  SRF_RETURN_DONE(funcctx);
1209  }
1210 
1211  for (i = 0; i < arg1->numbands; i++) {
1212  arg1->nbands[i] = i;
1213  POSTGIS_RT_DEBUGF(4, "arg1->nbands[%d] = %d", arg1->nbands[i], i);
1214  }
1215  }
1216  }
1217 
1218  /* pad (4) and padnodata (5) */
1219  if (!PG_ARGISNULL(4)) {
1220  arg1->pad.pad = PG_GETARG_BOOL(4) ? 1 : 0;
1221 
1222  if (arg1->pad.pad && !PG_ARGISNULL(5)) {
1223  arg1->pad.hasnodata = 1;
1224  arg1->pad.nodataval = PG_GETARG_FLOAT8(5);
1225  }
1226  else {
1227  arg1->pad.hasnodata = 0;
1228  arg1->pad.nodataval = 0;
1229  }
1230  }
1231  else {
1232  arg1->pad.pad = 0;
1233  arg1->pad.hasnodata = 0;
1234  arg1->pad.nodataval = 0;
1235  }
1236 
1237  /* store some additional metadata */
1238  arg1->raster.srid = rt_raster_get_srid(arg1->raster.raster);
1239  arg1->raster.width = rt_raster_get_width(arg1->raster.raster);
1240  arg1->raster.height = rt_raster_get_height(arg1->raster.raster);
1241  rt_raster_get_geotransform_matrix(arg1->raster.raster, arg1->raster.gt);
1242 
1243  /* determine maximum number of tiles from raster */
1244  arg1->tile.nx = ceil(arg1->raster.width / (double) arg1->tile.width);
1245  arg1->tile.ny = ceil(arg1->raster.height / (double) arg1->tile.height);
1246  POSTGIS_RT_DEBUGF(4, "# of tiles (x, y) = (%d, %d)", arg1->tile.nx, arg1->tile.ny);
1247 
1248  /* Store needed information */
1249  funcctx->user_fctx = arg1;
1250 
1251  /* total number of tuples to be returned */
1252  funcctx->max_calls = (arg1->tile.nx * arg1->tile.ny);
1253 
1254  MemoryContextSwitchTo(oldcontext);
1255  }
1256 
1257  /* stuff done on every call of the function */
1258  funcctx = SRF_PERCALL_SETUP();
1259 
1260  call_cntr = funcctx->call_cntr;
1261  max_calls = funcctx->max_calls;
1262  arg2 = funcctx->user_fctx;
1263 
1264  /* do when there is more left to send */
1265  if (call_cntr < max_calls) {
1266  rt_pgraster *pgtile = NULL;
1267  rt_raster tile = NULL;
1268  rt_band _band = NULL;
1269  rt_band band = NULL;
1270  rt_pixtype pixtype = PT_END;
1271  int hasnodata = 0;
1272  double nodataval = 0;
1273  int width = 0;
1274  int height = 0;
1275 
1276  int k = 0;
1277  int tx = 0;
1278  int ty = 0;
1279  int rx = 0;
1280  int ry = 0;
1281  int ex = 0; /* edge tile on right */
1282  int ey = 0; /* edge tile on bottom */
1283  double ulx = 0;
1284  double uly = 0;
1285  uint16_t len = 0;
1286  void *vals = NULL;
1287  uint16_t nvals;
1288 
1289  POSTGIS_RT_DEBUGF(3, "call number %d", call_cntr);
1290 
1291  /*
1292  find offset based upon tile #
1293 
1294  0 1 2
1295  3 4 5
1296  6 7 8
1297  */
1298  ty = call_cntr / arg2->tile.nx;
1299  tx = call_cntr % arg2->tile.nx;
1300  POSTGIS_RT_DEBUGF(4, "tile (x, y) = (%d, %d)", tx, ty);
1301 
1302  /* edge tile? only important if padding is false */
1303  if (!arg2->pad.pad) {
1304  if (ty + 1 == arg2->tile.ny)
1305  ey = 1;
1306  if (tx + 1 == arg2->tile.nx)
1307  ex = 1;
1308  }
1309 
1310  /* upper-left of tile in raster coordinates */
1311  rx = tx * arg2->tile.width;
1312  ry = ty * arg2->tile.height;
1313  POSTGIS_RT_DEBUGF(4, "raster coordinates = %d, %d", rx, ry);
1314 
1315  /* determine tile width and height */
1316  /* default to user-defined */
1317  width = arg2->tile.width;
1318  height = arg2->tile.height;
1319 
1320  /* override user-defined if edge tile (only possible if padding is false */
1321  if (ex || ey) {
1322  /* right edge */
1323  if (ex)
1324  width = arg2->raster.width - rx;
1325  /* bottom edge */
1326  if (ey)
1327  height = arg2->raster.height - ry;
1328  }
1329 
1330  /* create empty raster */
1331  tile = rt_raster_new(width, height);
1332  rt_raster_set_geotransform_matrix(tile, arg2->raster.gt);
1333  rt_raster_set_srid(tile, arg2->raster.srid);
1334 
1335  /* upper-left of tile in spatial coordinates */
1336  if (rt_raster_cell_to_geopoint(arg2->raster.raster, rx, ry, &ulx, &uly, arg2->raster.gt) != ES_NONE) {
1337  rt_raster_destroy(tile);
1338  rt_raster_destroy(arg2->raster.raster);
1339  if (arg2->numbands) pfree(arg2->nbands);
1340  pfree(arg2);
1341  elog(ERROR, "RASTER_tile: Could not compute the coordinates of the upper-left corner of the output tile");
1342  SRF_RETURN_DONE(funcctx);
1343  }
1344  rt_raster_set_offsets(tile, ulx, uly);
1345  POSTGIS_RT_DEBUGF(4, "spatial coordinates = %f, %f", ulx, uly);
1346 
1347  /* compute length of pixel line to read */
1348  len = arg2->tile.width;
1349  if (rx + arg2->tile.width >= arg2->raster.width)
1350  len = arg2->raster.width - rx;
1351  POSTGIS_RT_DEBUGF(3, "read line len = %d", len);
1352 
1353  /* copy bands to tile */
1354  for (i = 0; i < arg2->numbands; i++) {
1355  POSTGIS_RT_DEBUGF(4, "copying band %d to tile %d", arg2->nbands[i], call_cntr);
1356 
1357  _band = rt_raster_get_band(arg2->raster.raster, arg2->nbands[i]);
1358  if (_band == NULL) {
1359  int nband = arg2->nbands[i] + 1;
1360  rt_raster_destroy(tile);
1361  rt_raster_destroy(arg2->raster.raster);
1362  if (arg2->numbands) pfree(arg2->nbands);
1363  pfree(arg2);
1364  elog(ERROR, "RASTER_tile: Could not get band %d from source raster", nband);
1365  SRF_RETURN_DONE(funcctx);
1366  }
1367 
1368  pixtype = rt_band_get_pixtype(_band);
1369  hasnodata = rt_band_get_hasnodata_flag(_band);
1370  if (hasnodata)
1371  rt_band_get_nodata(_band, &nodataval);
1372  else if (arg2->pad.pad && arg2->pad.hasnodata) {
1373  hasnodata = 1;
1374  nodataval = arg2->pad.nodataval;
1375  }
1376  else
1377  nodataval = rt_band_get_min_value(_band);
1378 
1379  /* inline band */
1380  if (!rt_band_is_offline(_band)) {
1381  if (rt_raster_generate_new_band(tile, pixtype, nodataval, hasnodata, nodataval, i) < 0) {
1382  rt_raster_destroy(tile);
1383  rt_raster_destroy(arg2->raster.raster);
1384  pfree(arg2->nbands);
1385  pfree(arg2);
1386  elog(ERROR, "RASTER_tile: Could not add new band to output tile");
1387  SRF_RETURN_DONE(funcctx);
1388  }
1389  band = rt_raster_get_band(tile, i);
1390  if (band == NULL) {
1391  rt_raster_destroy(tile);
1392  rt_raster_destroy(arg2->raster.raster);
1393  if (arg2->numbands) pfree(arg2->nbands);
1394  pfree(arg2);
1395  elog(ERROR, "RASTER_tile: Could not get newly added band from output tile");
1396  SRF_RETURN_DONE(funcctx);
1397  }
1398 
1399  /* if isnodata, set flag and continue */
1400  if (rt_band_get_isnodata_flag(_band)) {
1401  rt_band_set_isnodata_flag(band, 1);
1402  continue;
1403  }
1404 
1405  /* copy data */
1406  for (j = 0; j < arg2->tile.height; j++) {
1407  k = ry + j;
1408 
1409  if (k >= arg2->raster.height) {
1410  POSTGIS_RT_DEBUGF(4, "row %d is beyond extent of source raster. skipping", k);
1411  continue;
1412  }
1413 
1414  POSTGIS_RT_DEBUGF(4, "getting pixel line %d, %d for %d pixels", rx, k, len);
1415  if (rt_band_get_pixel_line(_band, rx, k, len, &vals, &nvals) != ES_NONE) {
1416  rt_raster_destroy(tile);
1417  rt_raster_destroy(arg2->raster.raster);
1418  if (arg2->numbands) pfree(arg2->nbands);
1419  pfree(arg2);
1420  elog(ERROR, "RASTER_tile: Could not get pixel line from source raster");
1421  SRF_RETURN_DONE(funcctx);
1422  }
1423 
1424  if (nvals && rt_band_set_pixel_line(band, 0, j, vals, nvals) != ES_NONE) {
1425  rt_raster_destroy(tile);
1426  rt_raster_destroy(arg2->raster.raster);
1427  if (arg2->numbands) pfree(arg2->nbands);
1428  pfree(arg2);
1429  elog(ERROR, "RASTER_tile: Could not set pixel line of output tile");
1430  SRF_RETURN_DONE(funcctx);
1431  }
1432  }
1433  }
1434  /* offline */
1435  else {
1436  uint8_t bandnum = 0;
1437  rt_band_get_ext_band_num(_band, &bandnum);
1438 
1439  band = rt_band_new_offline(
1440  width, height,
1441  pixtype,
1442  hasnodata, nodataval,
1443  bandnum, rt_band_get_ext_path(_band)
1444  );
1445 
1446  if (band == NULL) {
1447  rt_raster_destroy(tile);
1448  rt_raster_destroy(arg2->raster.raster);
1449  if (arg2->numbands) pfree(arg2->nbands);
1450  pfree(arg2);
1451  elog(ERROR, "RASTER_tile: Could not create new offline band for output tile");
1452  SRF_RETURN_DONE(funcctx);
1453  }
1454 
1455  if (rt_raster_add_band(tile, band, i) < 0) {
1456  rt_band_destroy(band);
1457  rt_raster_destroy(tile);
1458  rt_raster_destroy(arg2->raster.raster);
1459  if (arg2->numbands) pfree(arg2->nbands);
1460  pfree(arg2);
1461  elog(ERROR, "RASTER_tile: Could not add new offline band to output tile");
1462  SRF_RETURN_DONE(funcctx);
1463  }
1464  }
1465  }
1466 
1467  pgtile = rt_raster_serialize(tile);
1468  rt_raster_destroy(tile);
1469  if (!pgtile) {
1470  rt_raster_destroy(arg2->raster.raster);
1471  if (arg2->numbands) pfree(arg2->nbands);
1472  pfree(arg2);
1473  SRF_RETURN_DONE(funcctx);
1474  }
1475 
1476  SET_VARSIZE(pgtile, pgtile->size);
1477  SRF_RETURN_NEXT(funcctx, PointerGetDatum(pgtile));
1478  }
1479  /* do when there is no more left */
1480  else {
1481  rt_raster_destroy(arg2->raster.raster);
1482  if (arg2->numbands) pfree(arg2->nbands);
1483  pfree(arg2);
1484  SRF_RETURN_DONE(funcctx);
1485  }
1486 }
1487 
1492 PG_FUNCTION_INFO_V1(RASTER_band);
1493 Datum RASTER_band(PG_FUNCTION_ARGS)
1494 {
1495  rt_pgraster *pgraster;
1496  rt_pgraster *pgrast;
1497  rt_raster raster;
1498  rt_raster rast;
1499 
1500  bool skip = FALSE;
1501  ArrayType *array;
1502  Oid etype;
1503  Datum *e;
1504  bool *nulls;
1505  int16 typlen;
1506  bool typbyval;
1507  char typalign;
1508 
1509  uint32_t numBands;
1510  uint32_t *bandNums;
1511  uint32 idx = 0;
1512  int n;
1513  int i = 0;
1514  int j = 0;
1515 
1516  if (PG_ARGISNULL(0))
1517  PG_RETURN_NULL();
1518  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1519 
1520  raster = rt_raster_deserialize(pgraster, FALSE);
1521  if (!raster) {
1522  PG_FREE_IF_COPY(pgraster, 0);
1523  elog(ERROR, "RASTER_band: Could not deserialize raster");
1524  PG_RETURN_NULL();
1525  }
1526 
1527  /* process bandNums */
1528  if (PG_ARGISNULL(1)) {
1529  elog(NOTICE, "Band number(s) not provided. Returning original raster");
1530  skip = TRUE;
1531  }
1532  do {
1533  if (skip) break;
1534 
1535  numBands = rt_raster_get_num_bands(raster);
1536 
1537  array = PG_GETARG_ARRAYTYPE_P(1);
1538  etype = ARR_ELEMTYPE(array);
1539  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1540 
1541  switch (etype) {
1542  case INT2OID:
1543  case INT4OID:
1544  break;
1545  default:
1546  rt_raster_destroy(raster);
1547  PG_FREE_IF_COPY(pgraster, 0);
1548  elog(ERROR, "RASTER_band: Invalid data type for band number(s)");
1549  PG_RETURN_NULL();
1550  break;
1551  }
1552 
1553  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1554  &nulls, &n);
1555 
1556  bandNums = palloc(sizeof(uint32_t) * n);
1557  for (i = 0, j = 0; i < n; i++) {
1558  if (nulls[i]) continue;
1559 
1560  switch (etype) {
1561  case INT2OID:
1562  idx = (uint32_t) DatumGetInt16(e[i]);
1563  break;
1564  case INT4OID:
1565  idx = (uint32_t) DatumGetInt32(e[i]);
1566  break;
1567  }
1568 
1569  POSTGIS_RT_DEBUGF(3, "band idx (before): %d", idx);
1570  if (idx > numBands || idx < 1) {
1571  elog(NOTICE, "Invalid band index (must use 1-based). Returning original raster");
1572  skip = TRUE;
1573  break;
1574  }
1575 
1576  bandNums[j] = idx - 1;
1577  POSTGIS_RT_DEBUGF(3, "bandNums[%d] = %d", j, bandNums[j]);
1578  j++;
1579  }
1580 
1581  if (skip || j < 1) {
1582  pfree(bandNums);
1583  skip = TRUE;
1584  }
1585  }
1586  while (0);
1587 
1588  if (!skip) {
1589  rast = rt_raster_from_band(raster, bandNums, j);
1590  pfree(bandNums);
1591  rt_raster_destroy(raster);
1592  PG_FREE_IF_COPY(pgraster, 0);
1593  if (!rast) {
1594  elog(ERROR, "RASTER_band: Could not create new raster");
1595  PG_RETURN_NULL();
1596  }
1597 
1598  pgrast = rt_raster_serialize(rast);
1599  rt_raster_destroy(rast);
1600 
1601  if (!pgrast)
1602  PG_RETURN_NULL();
1603 
1604  SET_VARSIZE(pgrast, pgrast->size);
1605  PG_RETURN_POINTER(pgrast);
1606  }
1607 
1608  PG_RETURN_POINTER(pgraster);
1609 }
rt_band_t
Definition: librtcore.h:2264
rt_raster_serialize
void * rt_raster_serialize(rt_raster raster)
Return this raster in serialized form.
Definition: rt_serialize.c:521
rt_util_gdal_register_all
int rt_util_gdal_register_all(int force_register_all)
Definition: rt_util.c:334
rt_raster_get_num_bands
int rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:372
rt_raster_same_alignment
rt_errorstate rt_raster_same_alignment(rt_raster rast1, rt_raster rast2, int *aligned, char **reason)
Definition: rt_spatial_relationship.c:45
rtrowdump.raster
raster
Be careful!! Zeros function&#39;s input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
rt_raster_from_band
rt_raster rt_raster_from_band(rt_raster raster, uint32_t *bandNums, int count)
Construct a new rt_raster from an existing rt_raster and an array of band numbers.
Definition: rt_raster.c:1439
rt_raster_set_geotransform_matrix
void rt_raster_set_geotransform_matrix(rt_raster raster, double *gt)
Set raster&#39;s geotransform using 6-element array.
Definition: rt_raster.c:727
rt_raster_get_geotransform_matrix
void rt_raster_get_geotransform_matrix(rt_raster raster, double *gt)
Get 6-element array of raster geotransform matrix.
Definition: rt_raster.c:706
rt_raster_set_skews
void rt_raster_set_skews(rt_raster raster, double skewX, double skewY)
Set skews about the X and Y axis.
Definition: rt_raster.c:168
ovdump.band
band
Definition: ovdump.py:57
rt_band_is_offline
int rt_band_is_offline(rt_band band)
Return non-zero if the given band data is on the filesystem.
Definition: rt_band.c:228
window.gt
gt
Definition: window.py:77
rt_pixtype
rt_pixtype
Definition: librtcore.h:185
rt_band_destroy
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition: rt_band.c:242
RASTER_addBandRasterArray
Datum RASTER_addBandRasterArray(PG_FUNCTION_ARGS)
Definition: rtpg_create.c:344
rt_raster_is_empty
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
Definition: rt_raster.c:1338
RASTER_band
Datum RASTER_band(PG_FUNCTION_ARGS)
Definition: rtpg_create.c:1493
POSTGIS_RT_DEBUGF
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition: rtpostgis.h:65
rt_band_get_nodata
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1597
rt_band_get_ext_path
const char * rt_band_get_ext_path(rt_band band)
Return band&#39;s external path (only valid when rt_band_is_offline returns non-zero).
Definition: rt_band.c:265
uint32_t
unsigned int uint32_t
Definition: uthash.h:78
rt_raster_cell_to_geopoint
rt_errorstate rt_raster_cell_to_geopoint(rt_raster raster, double xr, double yr, double *xw, double *yw, double *gt)
Convert an xr, yr raster point to an xw, yw point on map.
Definition: rt_raster.c:755
pixval.nband
nband
Definition: pixval.py:52
RASTER_copyBand
Datum RASTER_copyBand(PG_FUNCTION_ARGS)
Definition: rtpg_create.c:904
rt_raster_set_scale
void rt_raster_set_scale(rt_raster raster, double scaleX, double scaleY)
Set scale in projection units.
Definition: rt_raster.c:137
rt_raster_generate_new_band
int rt_raster_generate_new_band(rt_raster raster, rt_pixtype pixtype, double initialvalue, uint32_t hasnodata, double nodatavalue, int index)
Generate a new inline band and add it to a raster.
Definition: rt_raster.c:485
rt_pixtype_index_from_name
rt_pixtype rt_pixtype_index_from_name(const char *pixname)
Definition: rt_pixel.c:80
rt_band_get_pixel_line
rt_errorstate rt_band_get_pixel_line(rt_band band, int x, int y, uint16_t len, void **vals, uint16_t *nvals)
Get values of multiple pixels.
Definition: rt_band.c:1004
rt_raster_set_offsets
void rt_raster_set_offsets(rt_raster raster, double x, double y)
Set insertion points in projection units.
Definition: rt_raster.c:199
SRID_UNKNOWN
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:188
rt_raster_serialized_t::size
uint32_t size
Definition: librtcore.h:2203
rt_raster_get_band
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:381
RASTER_addBand
Datum RASTER_addBand(PG_FUNCTION_ARGS)
Definition: rtpg_create.c:121
rt_raster_copy_band
int rt_raster_copy_band(rt_raster torast, rt_raster fromrast, int fromindex, int toindex)
Copy one band from one raster to another.
Definition: rt_raster.c:1374
rt_util_gdal_datatype_to_pixtype
rt_pixtype rt_util_gdal_datatype_to_pixtype(GDALDataType gdt)
Convert GDALDataType to rt_pixtype.
Definition: rt_util.c:153
rt_band_get_hasnodata_flag
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:541
rt_raster_has_band
int rt_raster_has_band(rt_raster raster, int nband)
Return TRUE if the raster has a band of this number.
Definition: rt_raster.c:1351
RASTER_makeEmpty
Datum RASTER_makeEmpty(PG_FUNCTION_ARGS)
Definition: rtpg_create.c:55
rt_raster_set_srid
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster&#39;s SRID.
Definition: rt_raster.c:363
rt_raster_get_srid
int32_t rt_raster_get_srid(rt_raster raster)
Get raster&#39;s SRID.
Definition: rt_raster.c:356
rt_raster_destroy
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:82
rt_raster_new
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_raster.c:48
rt_band_get_ext_band_num
rt_errorstate rt_band_get_ext_band_num(rt_band band, uint8_t *bandnum)
Return bands&#39; external band number (only valid when rt_band_is_offline returns non-zero).
Definition: rt_band.c:278
RASTER_addBandOutDB
Datum RASTER_addBandOutDB(PG_FUNCTION_ARGS)
Definition: rtpg_create.c:546
rt_raster_get_width
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_raster.c:121
rt_band_set_isnodata_flag
rt_errorstate rt_band_set_isnodata_flag(rt_band band, int flag)
Set isnodata flag value.
Definition: rt_band.c:562
FALSE
#define FALSE
Definition: dbfopen.c:168
rt_raster_serialized_t
Struct definitions.
Definition: librtcore.h:2201
rt_band_get_min_value
double rt_band_get_min_value(rt_band band)
Returns the minimal possible value for the band according to the pixel type.
Definition: rt_band.c:1612
rt_pixtype_name
const char * rt_pixtype_name(rt_pixtype pixtype)
Definition: rt_pixel.c:110
RASTER_tile
Datum RASTER_tile(PG_FUNCTION_ARGS)
Definition: rtpg_create.c:984
rt_band_get_pixtype
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition: rt_band.c:498
rt_band_set_pixel_line
rt_errorstate rt_band_set_pixel_line(rt_band band, int x, int y, void *vals, uint32_t len)
Set values of multiple pixels.
Definition: rt_band.c:720
rt_raster_add_band
int rt_raster_add_band(rt_raster raster, rt_band band, int index)
Add band data to a raster.
Definition: rt_raster.c:405
rt_util_gdal_sr_auth_info
rt_errorstate rt_util_gdal_sr_auth_info(GDALDatasetH hds, char **authname, char **authcode)
Get auth name and code.
Definition: rt_util.c:270
rt_band_get_isnodata_flag
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition: rt_band.c:581
rt_raster_t
Definition: librtcore.h:2235
rt_raster_get_height
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_raster.c:129
PG_FUNCTION_INFO_V1
PG_FUNCTION_INFO_V1(RASTER_makeEmpty)
Make a new raster with no bands.
POSTGIS_RT_DEBUG
#define POSTGIS_RT_DEBUG(level, msg)
Definition: rtpostgis.h:61
uint8_t
unsigned char uint8_t
Definition: uthash.h:79
TRUE
#define TRUE
Definition: dbfopen.c:169
rt_raster_deserialize
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition: rt_serialize.c:717
rt_util_gdal_open
GDALDatasetH rt_util_gdal_open(const char *fn, GDALAccess fn_access, int shared)
Definition: rt_util.c:379
rt_band_new_offline
rt_band rt_band_new_offline(uint16_t width, uint16_t height, rt_pixtype pixtype, uint32_t hasnodata, double nodataval, uint8_t bandNum, const char *path)
Create an out-db rt_band.
Definition: rt_band.c:119