PostGIS  2.1.10dev-r@@SVN_REVISION@@
rt_pg.c
Go to the documentation of this file.
1 /*
2  * $Id: rt_pg.c 14705 2016-02-26 11:37:39Z pramsey $
3  *
4  * WKTRaster - Raster Types for PostGIS
5  * http://trac.osgeo.org/postgis/wiki/WKTRaster
6  *
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@keybit.net>
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 <math.h>
32 #include <string.h>
33 #include <stdio.h>
34 #include <stdlib.h> /* for strtod in RASTER_reclass */
35 #include <errno.h>
36 #include <assert.h>
37 #include <ctype.h> /* for isspace */
38 
39 #include <postgres.h> /* for palloc */
40 #include <access/gist.h>
41 #include <access/itup.h>
42 #include <fmgr.h>
43 #include <utils/elog.h>
44 #include <utils/builtins.h>
45 #include <executor/spi.h>
46 #include <executor/executor.h> /* for GetAttributeByName in RASTER_reclass */
47 #include <funcapi.h>
48 
49 #include "../../postgis_config.h"
50 
51 #include "lwgeom_pg.h"
52 #include "rt_pg.h"
53 #include "pgsql_compat.h"
54 
55 #include "utils/lsyscache.h" /* for get_typlenbyvalalign */
56 #include "utils/array.h" /* for ArrayType */
57 #include "catalog/pg_type.h" /* for INT2OID, INT4OID, FLOAT4OID, FLOAT8OID and TEXTOID */
58 
59 #if POSTGIS_PGSQL_VERSION > 92
60 #include "access/htup_details.h"
61 #endif
62 
63 /* maximum char length required to hold any double or long long value */
64 #define MAX_DBL_CHARLEN (3 + DBL_MANT_DIG - DBL_MIN_EXP)
65 #define MAX_INT_CHARLEN 32
66 
67 /*
68  * This is required for builds against pgsql
69  */
71 
72 void _PG_init(void);
73 
74 /***************************************************************
75  * Internal functions must be prefixed with rtpg_. This is
76  * keeping inline with the use of pgis_ for ./postgis C utility
77  * functions.
78  ****************************************************************/
79 /* Internal funcs */
80 static char *rtpg_strreplace(
81  const char *str,
82  const char *oldstr, const char *newstr,
83  int *count
84 );
85 static char *rtpg_strtoupper(char *str);
86 static char *rtpg_chartrim(const char* input, char *remove);
87 static char **rtpg_strsplit(const char *str, const char *delimiter, int *n);
88 static char *rtpg_removespaces(char *str);
89 static char *rtpg_trim(const char* input);
90 static char *rtpg_getSR(int srid);
91 
92 extern char *gdal_enabled_drivers;
93 extern char enable_outdb_rasters;
94 
95 #define ENV_POSTGIS_GDAL_ENABLED_DRIVERS "POSTGIS_GDAL_ENABLED_DRIVERS"
96 #define ENV_POSTGIS_ENABLE_OUTDB_RASTERS "POSTGIS_ENABLE_OUTDB_RASTERS"
97 
98 
99 /* ---------------------------------------------------------------- */
100 /* Memory allocation / error reporting hooks */
101 /* ---------------------------------------------------------------- */
102 
103 static void *
104 rt_pg_alloc(size_t size)
105 {
106  void * result;
107 
108  POSTGIS_RT_DEBUGF(5, "rt_pgalloc(%ld) called", (long int) size);
109 
110  result = palloc(size);
111 
112  return result;
113 }
114 
115 static void *
116 rt_pg_realloc(void *mem, size_t size)
117 {
118  void * result;
119 
120  POSTGIS_RT_DEBUGF(5, "rt_pg_realloc(%ld) called", (long int) size);
121 
122  if (mem)
123  result = repalloc(mem, size);
124 
125  else
126  result = palloc(size);
127 
128  return result;
129 }
130 
131 static void
132 rt_pg_free(void *ptr)
133 {
134  POSTGIS_RT_DEBUG(5, "rt_pfree called");
135  pfree(ptr);
136 }
137 
138 static void
139 rt_pg_error(const char *fmt, va_list ap)
140 {
141 #define ERRMSG_MAXLEN 256
142 
143  char errmsg[ERRMSG_MAXLEN+1];
144 
145  vsnprintf (errmsg, ERRMSG_MAXLEN, fmt, ap);
146 
147  errmsg[ERRMSG_MAXLEN]='\0';
148  ereport(ERROR, (errmsg_internal("%s", errmsg)));
149 }
150 
151 static void
152 rt_pg_notice(const char *fmt, va_list ap)
153 {
154  char *msg;
155 
156  /*
157  * This is a GNU extension.
158  * Dunno how to handle errors here.
159  */
160  if (!lw_vasprintf (&msg, fmt, ap))
161  {
162  va_end (ap);
163  return;
164  }
165  ereport(NOTICE, (errmsg_internal("%s", msg)));
166  free(msg);
167 }
168 
169 
170 /* rtpg_assignHookGDALEnabledDrivers() should only be called by _PG_init */
171 static void
173  char *enabled_drivers = NULL;
174  int enable_all = 0;
175  int disable_all = 0;
176 
177  char **enabled_drivers_array = NULL;
178  int enabled_drivers_count = 0;
179  bool *enabled_drivers_found = NULL;
180  char *gdal_skip = NULL;
181 
182  uint32_t i;
183  uint32_t j;
184 
185  enabled_drivers = getenv(ENV_POSTGIS_GDAL_ENABLED_DRIVERS);
186 
187  POSTGIS_RT_DEBUGF(4, "GDAL_SKIP = \"%s\"", CPLGetConfigOption("GDAL_SKIP", ""));
188  POSTGIS_RT_DEBUGF(4, "enabled_drivers = \"%s\"", enabled_drivers);
189 
190  if (enabled_drivers != NULL) {
191  gdal_enabled_drivers = palloc(
192  sizeof(char) * (strlen(enabled_drivers) + 1)
193  );
194 
195  sprintf(gdal_enabled_drivers, "%s", enabled_drivers);
196 
197  enabled_drivers_array = rtpg_strsplit(enabled_drivers, " ", &enabled_drivers_count);
198 #if POSTGIS_DEBUG_LEVEL > 0
199  POSTGIS_RT_DEBUGF(4, "enabled_drivers_count = %d", enabled_drivers_count)
200  for (i = 0; i < enabled_drivers_count; i++) {
201  POSTGIS_RT_DEBUGF(4, "enabled_drivers_array[%d] = \"%s\"", i, enabled_drivers_array[i]);
202  }
203 #endif
204 
205  enabled_drivers_found = palloc(sizeof(bool) * enabled_drivers_count);
206  memset(enabled_drivers_found, FALSE, sizeof(bool) * enabled_drivers_count);
207  }
208  else {
209  gdal_enabled_drivers = palloc(sizeof(char));
210  gdal_enabled_drivers[0] = '\0';
211  }
212 
213  /* destroy the driver manager */
214  /* this is the only way to ensure GDAL_SKIP is recognized */
215  GDALDestroyDriverManager();
216  CPLSetConfigOption("GDAL_SKIP", NULL);
217 
218  /* force wrapper function to call GDALAllRegister() */
220 
221  /* scan for keywords DISABLE_ALL and ENABLE_ALL */
222  disable_all = 0;
223  enable_all = 0;
224  if (strstr(gdal_enabled_drivers, GDAL_DISABLE_ALL) != NULL) {
225  for (i = 0; i < enabled_drivers_count; i++) {
226  if (strstr(enabled_drivers_array[i], GDAL_DISABLE_ALL) != NULL) {
227  enabled_drivers_found[i] = TRUE;
228  disable_all = 1;
229  }
230  }
231  }
232  else if (strstr(gdal_enabled_drivers, GDAL_ENABLE_ALL) != NULL) {
233  for (i = 0; i < enabled_drivers_count; i++) {
234  if (strstr(enabled_drivers_array[i], GDAL_ENABLE_ALL) != NULL) {
235  enabled_drivers_found[i] = TRUE;
236  enable_all = 1;
237  }
238  }
239  }
240 
241  if (!enable_all) {
242  int found = 0;
243  uint32_t drv_count = 0;
244  rt_gdaldriver drv_set = rt_raster_gdal_drivers(&drv_count, 0);
245 
246  POSTGIS_RT_DEBUGF(4, "driver count = %d", drv_count);
247 
248  /* all other drivers than those in new drivers are added to GDAL_SKIP */
249  for (i = 0; i < drv_count; i++) {
250  POSTGIS_RT_DEBUGF(4, "drv_set[%d] = \"%s\"", i, drv_set[i].short_name);
251  found = 0;
252 
253  if (!disable_all) {
254  /* gdal driver found in gdal_enabled_drivers, continue to thorough search */
255  if (strstr(gdal_enabled_drivers, drv_set[i].short_name) != NULL) {
256  POSTGIS_RT_DEBUGF(4, "\"%s\" found in gdal_enabled_drivers", drv_set[i].short_name);
257  /* thorough search of enabled_drivers */
258  for (j = 0; j < enabled_drivers_count; j++) {
259  /* driver found */
260  if (strcmp(enabled_drivers_array[j], drv_set[i].short_name) == 0) {
261  POSTGIS_RT_DEBUGF(4, "\"%s\" found in enabled_drivers_array", drv_set[i].short_name);
262  enabled_drivers_found[j] = TRUE;
263  found = 1;
264  }
265  }
266  }
267  }
268 
269  /* driver found, continue */
270  if (found)
271  continue;
272 
273  /* driver not found, add to gdal_skip */
274  if (gdal_skip == NULL) {
275  gdal_skip = palloc(sizeof(char) * (strlen(drv_set[i].short_name) + 1));
276  gdal_skip[0] = '\0';
277  }
278  else {
279  gdal_skip = repalloc(
280  gdal_skip,
281  sizeof(char) * (
282  strlen(gdal_skip) + 1 + strlen(drv_set[i].short_name) + 1
283  )
284  );
285  strcat(gdal_skip, " ");
286  }
287  strcat(gdal_skip, drv_set[i].short_name);
288  }
289 
290  for (i = 0; i < drv_count; i++) {
291  pfree(drv_set[i].short_name);
292  pfree(drv_set[i].long_name);
293  pfree(drv_set[i].create_options);
294  }
295  if (drv_count) pfree(drv_set);
296 
297  }
298 
299  for (i = 0; i < enabled_drivers_count; i++) {
300  if (enabled_drivers_found[i])
301  continue;
302 
303  if (disable_all)
304  elog(WARNING, "%s set. Ignoring GDAL driver: %s", GDAL_DISABLE_ALL, enabled_drivers_array[i]);
305  else if (enable_all)
306  elog(WARNING, "%s set. Ignoring GDAL driver: %s", GDAL_ENABLE_ALL, enabled_drivers_array[i]);
307  else
308  elog(WARNING, "Unknown GDAL driver: %s", enabled_drivers_array[i]);
309  }
310 
311  /* destroy the driver manager */
312  /* this is the only way to ensure GDAL_SKIP is recognized */
313  GDALDestroyDriverManager();
314 
315  /* set GDAL_SKIP */
316  POSTGIS_RT_DEBUGF(4, "gdal_skip = \"%s\"", gdal_skip);
317  CPLSetConfigOption("GDAL_SKIP", gdal_skip);
318  if (gdal_skip != NULL) pfree(gdal_skip);
319 
320  /* force wrapper function to call GDALAllRegister() */
322 
323  if (enabled_drivers_count) {
324  pfree(enabled_drivers_array);
325  pfree(enabled_drivers_found);
326  }
327  POSTGIS_RT_DEBUGF(4, "GDAL_SKIP = \"%s\"", CPLGetConfigOption("GDAL_SKIP", ""));
328 }
329 
330 /*
331  * Module load callback
332  */
333 void _PG_init(void) {
334  char *env_postgis_enable_outdb_rasters = NULL;
335 
336  /* Install liblwgeom handlers */
337  pg_install_lwgeom_handlers();
338 
339  /* Install rt_api handlers */
341 
342  /*
343  * use POSTGIS_ENABLE_OUTDB_RASTERS to enable access to out-db rasters
344  */
346  env_postgis_enable_outdb_rasters = getenv(ENV_POSTGIS_ENABLE_OUTDB_RASTERS);
347  if (env_postgis_enable_outdb_rasters != NULL) {
348  char *env = rtpg_trim(env_postgis_enable_outdb_rasters);
349 
350  /* out of memory */
351  if (env == NULL) {
352  elog(
353  ERROR,
354  "_PG_init: Cannot process environmental variable: %s",
356  );
357  return;
358  }
359 
360  if (strcmp(env, "1") == 0)
362 
363  pfree(env);
364  }
365 
366  /*
367  * use POSTGIS_GDAL_ENABLED_DRIVERS to restrict drivers
368  */
370 
371 }
372 
373 /***************************************************************
374  * Some rules for returning NOTICE or ERROR...
375  *
376  * Send an ERROR like:
377  *
378  * elog(ERROR, "RASTER_out: Could not deserialize raster");
379  *
380  * only when:
381  *
382  * -something wrong happen with memory,
383  * -a function got an invalid argument ('3BUI' as pixel type) so that no row can
384  * be processed
385  *
386  * *** IMPORTANT: elog(ERROR, ...) does NOT return to calling function ***
387  *
388  * Send a NOTICE like:
389  *
390  * elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
391  *
392  * when arguments (e.g. x, y, band) are NULL or out of range so that some or
393  * most rows can be processed anyway
394  *
395  * in this case,
396  * for SET functions or function normally returning a modified raster, return
397  * the original raster
398  * for GET functions, return NULL
399  * try to deduce a valid parameter value if it makes sence (e.g. out of range
400  * index for addBand)
401  *
402  * Do not put the name of the faulty function for NOTICEs, only with ERRORs.
403  *
404  ****************************************************************/
405 
406 /******************************************************************************
407  * Some notes on memory management...
408  *
409  * Every time a SQL function is called, PostgreSQL creates a new memory context.
410  * So, all the memory allocated with palloc/repalloc in that context is
411  * automatically free'd at the end of the function. If you want some data to
412  * live between function calls, you have 2 options:
413  *
414  * - Use fcinfo->flinfo->fn_mcxt contex to store the data (by pointing the
415  * data you want to keep with fcinfo->flinfo->fn_extra)
416  * - Use SRF funcapi, and storing the data at multi_call_memory_ctx (by pointing
417  * the data you want to keep with funcctx->user_fctx. funcctx is created by
418  * funcctx = SPI_FIRSTCALL_INIT()). Recommended way in functions returning rows,
419  * like RASTER_dumpAsPolygons (see section 34.9.9 at
420  * http://www.postgresql.org/docs/8.4/static/xfunc-c.html).
421  *
422  * But raster code follows the same philosophy than the rest of PostGIS: keep
423  * memory as clean as possible. So, we free all allocated memory.
424  *
425  * TODO: In case of functions returning NULL, we should free the memory too.
426  *****************************************************************************/
427 
428 /******************************************************************************
429  * Notes for use of PG_DETOAST_DATUM(), PG_DETOAST_DATUM_SLICE()
430  * and PG_DETOAST_DATUM_COPY()
431  *
432  * When ONLY getting raster (not band) metadata, use PG_DETOAST_DATUM_SLICE()
433  * as it is generally quicker to get only the chunk of memory that contains
434  * the raster metadata.
435  *
436  * Example: PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(0), 0,
437  * sizeof(struct rt_raster_serialized_t))
438  *
439  * When ONLY setting raster or band(s) metadata OR reading band data, use
440  * PG_DETOAST_DATUM() as rt_raster_deserialize() allocates local memory
441  * for the raster and band(s) metadata.
442  *
443  * Example: PG_DETOAST_DATUM(PG_GETARG_DATUM(0))
444  *
445  * When SETTING band pixel values, use PG_DETOAST_DATUM_COPY(). This is
446  * because band data (not metadata) is just a pointer to the correct
447  * memory location in the detoasted datum. What is returned from
448  * PG_DETOAST_DATUM() may or may not be a copy of the input datum.
449  *
450  * From the comments in postgresql/src/include/fmgr.h...
451  *
452  * pg_detoast_datum() gives you either the input datum (if not toasted)
453  * or a detoasted copy allocated with palloc().
454  *
455  * From the mouth of Tom Lane...
456  * http://archives.postgresql.org/pgsql-hackers/2002-01/msg01289.php
457  *
458  * PG_DETOAST_DATUM_COPY guarantees to give you a copy, even if the
459  * original wasn't toasted. This allows you to scribble on the input,
460  * in case that happens to be a useful way of forming your result.
461  * Without a forced copy, a routine for a pass-by-ref datatype must
462  * NEVER, EVER scribble on its input ... because very possibly it'd
463  * be scribbling on a valid tuple in a disk buffer, or a valid entry
464  * in the syscache.
465  *
466  * The key detail above is that the raster datatype is a varlena, a
467  * passed by reference datatype.
468  *
469  * Example: PG_DETOAST_DATUM_COPY(PG_GETARG_DATUM(0))
470  *
471  * If in doubt, use PG_DETOAST_DATUM_COPY() as that guarantees that the input
472  * datum is copied for use.
473  *****************************************************************************/
474 
475 /* Prototypes */
476 
477 /* Utility functions */
478 Datum RASTER_lib_version(PG_FUNCTION_ARGS);
479 Datum RASTER_lib_build_date(PG_FUNCTION_ARGS);
480 Datum RASTER_gdal_version(PG_FUNCTION_ARGS);
481 Datum RASTER_minPossibleValue(PG_FUNCTION_ARGS);
482 
483 /* Input/output and format conversions */
484 Datum RASTER_in(PG_FUNCTION_ARGS);
485 Datum RASTER_out(PG_FUNCTION_ARGS);
486 
487 Datum RASTER_to_bytea(PG_FUNCTION_ARGS);
488 Datum RASTER_to_binary(PG_FUNCTION_ARGS);
489 
490 /* Raster as geometry operations */
491 Datum RASTER_convex_hull(PG_FUNCTION_ARGS);
492 Datum RASTER_dumpAsPolygons(PG_FUNCTION_ARGS);
493 
494 /* Get all the properties of a raster */
495 Datum RASTER_getSRID(PG_FUNCTION_ARGS);
496 Datum RASTER_getWidth(PG_FUNCTION_ARGS);
497 Datum RASTER_getHeight(PG_FUNCTION_ARGS);
498 Datum RASTER_getNumBands(PG_FUNCTION_ARGS);
499 Datum RASTER_getXScale(PG_FUNCTION_ARGS);
500 Datum RASTER_getYScale(PG_FUNCTION_ARGS);
501 Datum RASTER_getXSkew(PG_FUNCTION_ARGS);
502 Datum RASTER_getYSkew(PG_FUNCTION_ARGS);
503 Datum RASTER_getXUpperLeft(PG_FUNCTION_ARGS);
504 Datum RASTER_getYUpperLeft(PG_FUNCTION_ARGS);
505 Datum RASTER_getPixelWidth(PG_FUNCTION_ARGS);
506 Datum RASTER_getPixelHeight(PG_FUNCTION_ARGS);
507 Datum RASTER_getGeotransform(PG_FUNCTION_ARGS);
508 
509 /* Set all the properties of a raster */
510 Datum RASTER_setSRID(PG_FUNCTION_ARGS);
511 Datum RASTER_setScale(PG_FUNCTION_ARGS);
512 Datum RASTER_setScaleXY(PG_FUNCTION_ARGS);
513 Datum RASTER_setSkew(PG_FUNCTION_ARGS);
514 Datum RASTER_setSkewXY(PG_FUNCTION_ARGS);
515 Datum RASTER_setUpperLeftXY(PG_FUNCTION_ARGS);
516 Datum RASTER_setRotation(PG_FUNCTION_ARGS);
517 Datum RASTER_setGeotransform(PG_FUNCTION_ARGS);
518 
519 /* Get all the properties of a raster band */
520 Datum RASTER_getBandPixelType(PG_FUNCTION_ARGS);
521 Datum RASTER_getBandPixelTypeName(PG_FUNCTION_ARGS);
522 Datum RASTER_getBandNoDataValue(PG_FUNCTION_ARGS);
523 Datum RASTER_getBandPath(PG_FUNCTION_ARGS);
524 Datum RASTER_bandIsNoData(PG_FUNCTION_ARGS);
525 Datum RASTER_isEmpty(PG_FUNCTION_ARGS);
526 Datum RASTER_hasNoBand(PG_FUNCTION_ARGS);
527 
528 /* Set all the properties of a raster band */
529 Datum RASTER_setBandIsNoData(PG_FUNCTION_ARGS);
530 Datum RASTER_setBandNoDataValue(PG_FUNCTION_ARGS);
531 
532 /* Get pixel value */
533 Datum RASTER_getPixelValue(PG_FUNCTION_ARGS);
534 Datum RASTER_dumpValues(PG_FUNCTION_ARGS);
535 
536 /* Set pixel value(s) */
537 Datum RASTER_setPixelValue(PG_FUNCTION_ARGS);
538 Datum RASTER_setPixelValuesArray(PG_FUNCTION_ARGS);
539 Datum RASTER_setPixelValuesGeomval(PG_FUNCTION_ARGS);
540 
541 /* Get pixel geographical shape */
542 Datum RASTER_getPixelPolygons(PG_FUNCTION_ARGS);
543 
544 /* Get raster band's polygon */
545 Datum RASTER_getPolygon(PG_FUNCTION_ARGS);
546 
547 /* Get pixels of value */
548 Datum RASTER_pixelOfValue(PG_FUNCTION_ARGS);
549 
550 /* Get nearest value to a point */
551 Datum RASTER_nearestValue(PG_FUNCTION_ARGS);
552 
553 /* Get the neighborhood around a pixel */
554 Datum RASTER_neighborhood(PG_FUNCTION_ARGS);
555 
556 /* Raster and band creation */
557 Datum RASTER_makeEmpty(PG_FUNCTION_ARGS);
558 Datum RASTER_addBand(PG_FUNCTION_ARGS);
559 Datum RASTER_copyBand(PG_FUNCTION_ARGS);
560 Datum RASTER_addBandRasterArray(PG_FUNCTION_ARGS);
561 Datum RASTER_addBandOutDB(PG_FUNCTION_ARGS);
562 Datum RASTER_tile(PG_FUNCTION_ARGS);
563 
564 /* create new raster from existing raster's bands */
565 Datum RASTER_band(PG_FUNCTION_ARGS);
566 
567 /* Get summary stats */
568 Datum RASTER_summaryStats(PG_FUNCTION_ARGS);
569 Datum RASTER_summaryStatsCoverage(PG_FUNCTION_ARGS);
570 
571 /* get histogram */
572 Datum RASTER_histogram(PG_FUNCTION_ARGS);
573 Datum RASTER_histogramCoverage(PG_FUNCTION_ARGS);
574 
575 /* get quantiles */
576 Datum RASTER_quantile(PG_FUNCTION_ARGS);
577 Datum RASTER_quantileCoverage(PG_FUNCTION_ARGS);
578 
579 /* get counts of values */
580 Datum RASTER_valueCount(PG_FUNCTION_ARGS);
581 Datum RASTER_valueCountCoverage(PG_FUNCTION_ARGS);
582 
583 /* reclassify specified bands of a raster */
584 Datum RASTER_reclass(PG_FUNCTION_ARGS);
585 
586 /* apply colormap to specified band of a raster */
587 Datum RASTER_colorMap(PG_FUNCTION_ARGS);
588 
589 /* convert GDAL raster to raster */
590 Datum RASTER_fromGDALRaster(PG_FUNCTION_ARGS);
591 
592 /* convert raster to GDAL raster */
593 Datum RASTER_asGDALRaster(PG_FUNCTION_ARGS);
594 Datum RASTER_getGDALDrivers(PG_FUNCTION_ARGS);
595 
596 /* rasterize a geometry */
597 Datum RASTER_asRaster(PG_FUNCTION_ARGS);
598 
599 /* warp a raster using GDAL Warp API */
600 Datum RASTER_GDALWarp(PG_FUNCTION_ARGS);
601 
602 /* get raster's meta data */
603 Datum RASTER_metadata(PG_FUNCTION_ARGS);
604 
605 /* get raster band's meta data */
606 Datum RASTER_bandmetadata(PG_FUNCTION_ARGS);
607 
608 /* convert pixel/line to spatial coordinates */
609 Datum RASTER_rasterToWorldCoord(PG_FUNCTION_ARGS);
610 
611 /* convert spatial coordinates to pixel/line*/
612 Datum RASTER_worldToRasterCoord(PG_FUNCTION_ARGS);
613 
614 /* determine if two rasters intersect */
615 Datum RASTER_intersects(PG_FUNCTION_ARGS);
616 
617 /* determine if two rasters overlap */
618 Datum RASTER_overlaps(PG_FUNCTION_ARGS);
619 
620 /* determine if two rasters touch */
621 Datum RASTER_touches(PG_FUNCTION_ARGS);
622 
623 /* determine if the first raster contains the second raster */
624 Datum RASTER_contains(PG_FUNCTION_ARGS);
625 
626 /* determine if the first raster contains properly the second raster */
627 Datum RASTER_containsProperly(PG_FUNCTION_ARGS);
628 
629 /* determine if the first raster covers the second raster */
630 Datum RASTER_covers(PG_FUNCTION_ARGS);
631 
632 /* determine if the first raster is covered by the second raster */
633 Datum RASTER_coveredby(PG_FUNCTION_ARGS);
634 
635 /* determine if the two rasters are within the specified distance of each other */
636 Datum RASTER_dwithin(PG_FUNCTION_ARGS);
637 
638 /* determine if the two rasters are fully within the specified distance of each other */
639 Datum RASTER_dfullywithin(PG_FUNCTION_ARGS);
640 
641 /* determine if two rasters are aligned */
642 Datum RASTER_sameAlignment(PG_FUNCTION_ARGS);
643 Datum RASTER_notSameAlignmentReason(PG_FUNCTION_ARGS);
644 
645 /* one-raster MapAlgebra */
646 Datum RASTER_mapAlgebraExpr(PG_FUNCTION_ARGS);
647 Datum RASTER_mapAlgebraFct(PG_FUNCTION_ARGS);
648 
649 /* one-raster neighborhood MapAlgebra */
650 Datum RASTER_mapAlgebraFctNgb(PG_FUNCTION_ARGS);
651 
652 /* two-raster MapAlgebra */
653 Datum RASTER_mapAlgebra2(PG_FUNCTION_ARGS);
654 
655 /* n-raster MapAlgebra */
656 Datum RASTER_nMapAlgebra(PG_FUNCTION_ARGS);
657 Datum RASTER_nMapAlgebraExpr(PG_FUNCTION_ARGS);
658 
659 /* raster union aggregate */
660 Datum RASTER_union_transfn(PG_FUNCTION_ARGS);
661 Datum RASTER_union_finalfn(PG_FUNCTION_ARGS);
662 
663 /* raster clip */
664 Datum RASTER_clip(PG_FUNCTION_ARGS);
665 
666 /* string replacement function taken from
667  * http://ubuntuforums.org/showthread.php?s=aa6f015109fd7e4c7e30d2fd8b717497&t=141670&page=3
668  */
669 /* ---------------------------------------------------------------------------
670  Name : replace - Search & replace a substring by another one.
671  Creation : Thierry Husson, Sept 2010
672  Parameters :
673  str : Big string where we search
674  oldstr : Substring we are looking for
675  newstr : Substring we want to replace with
676  count : Optional pointer to int (input / output value). NULL to ignore.
677  Input: Maximum replacements to be done. NULL or < 1 to do all.
678  Output: Number of replacements done or -1 if not enough memory.
679  Returns : Pointer to the new string or NULL if error.
680  Notes :
681  - Case sensitive - Otherwise, replace functions "strstr" by "strcasestr"
682  - Always allocates memory for the result.
683 --------------------------------------------------------------------------- */
684 static char*
686  const char *str,
687  const char *oldstr, const char *newstr,
688  int *count
689 ) {
690  const char *tmp = str;
691  char *result;
692  int found = 0;
693  int length, reslen;
694  int oldlen = strlen(oldstr);
695  int newlen = strlen(newstr);
696  int limit = (count != NULL && *count > 0) ? *count : -1;
697 
698  tmp = str;
699  while ((tmp = strstr(tmp, oldstr)) != NULL && found != limit)
700  found++, tmp += oldlen;
701 
702  length = strlen(str) + found * (newlen - oldlen);
703  if ((result = (char *) palloc(length + 1)) == NULL) {
704  fprintf(stderr, "Not enough memory\n");
705  found = -1;
706  }
707  else {
708  tmp = str;
709  limit = found; /* Countdown */
710  reslen = 0; /* length of current result */
711 
712  /* Replace each old string found with new string */
713  while ((limit-- > 0) && (tmp = strstr(tmp, oldstr)) != NULL) {
714  length = (tmp - str); /* Number of chars to keep intouched */
715  strncpy(result + reslen, str, length); /* Original part keeped */
716  strcpy(result + (reslen += length), newstr); /* Insert new string */
717 
718  reslen += newlen;
719  tmp += oldlen;
720  str = tmp;
721  }
722  strcpy(result + reslen, str); /* Copies last part and ending null char */
723  }
724 
725  if (count != NULL) *count = found;
726  return result;
727 }
728 
729 static char *
730 rtpg_strtoupper(char * str) {
731  int j;
732 
733  for (j = strlen(str) - 1; j >= 0; j--)
734  str[j] = toupper(str[j]);
735 
736  return str;
737 }
738 
739 static char*
740 rtpg_chartrim(const char *input, char *remove) {
741  char *rtn = NULL;
742  char *ptr = NULL;
743  uint32_t offset = 0;
744 
745  if (!input)
746  return NULL;
747  else if (!*input)
748  return (char *) input;
749 
750  /* trim left */
751  while (strchr(remove, *input) != NULL)
752  input++;
753 
754  /* trim right */
755  ptr = ((char *) input) + strlen(input);
756  while (strchr(remove, *--ptr) != NULL)
757  offset++;
758 
759  rtn = palloc(sizeof(char) * (strlen(input) - offset + 1));
760  if (rtn == NULL) {
761  fprintf(stderr, "Not enough memory\n");
762  return NULL;
763  }
764  strncpy(rtn, input, strlen(input) - offset);
765  rtn[strlen(input) - offset] = '\0';
766 
767  return rtn;
768 }
769 
770 /* split a string based on a delimiter */
771 static char**
772 rtpg_strsplit(const char *str, const char *delimiter, int *n) {
773  char *tmp = NULL;
774  char **rtn = NULL;
775  char *token = NULL;
776 
777  *n = 0;
778  if (!str)
779  return NULL;
780 
781  /* copy str to tmp as strtok will mangle the string */
782  tmp = palloc(sizeof(char) * (strlen(str) + 1));
783  if (NULL == tmp) {
784  fprintf(stderr, "Not enough memory\n");
785  return NULL;
786  }
787  strcpy(tmp, str);
788 
789  if (!strlen(tmp) || !delimiter || !strlen(delimiter)) {
790  *n = 1;
791  rtn = (char **) palloc(*n * sizeof(char *));
792  if (NULL == rtn) {
793  fprintf(stderr, "Not enough memory\n");
794  return NULL;
795  }
796  rtn[0] = (char *) palloc(sizeof(char) * (strlen(tmp) + 1));
797  if (NULL == rtn[0]) {
798  fprintf(stderr, "Not enough memory\n");
799  return NULL;
800  }
801  strcpy(rtn[0], tmp);
802  pfree(tmp);
803  return rtn;
804  }
805 
806  token = strtok(tmp, delimiter);
807  while (token != NULL) {
808  if (*n < 1) {
809  rtn = (char **) palloc(sizeof(char *));
810  }
811  else {
812  rtn = (char **) repalloc(rtn, (*n + 1) * sizeof(char *));
813  }
814  if (NULL == rtn) {
815  fprintf(stderr, "Not enough memory\n");
816  return NULL;
817  }
818 
819  rtn[*n] = NULL;
820  rtn[*n] = (char *) palloc(sizeof(char) * (strlen(token) + 1));
821  if (NULL == rtn[*n]) {
822  fprintf(stderr, "Not enough memory\n");
823  return NULL;
824  }
825 
826  strcpy(rtn[*n], token);
827  *n = *n + 1;
828 
829  token = strtok(NULL, delimiter);
830  }
831 
832  pfree(tmp);
833  return rtn;
834 }
835 
836 static char *
837 rtpg_removespaces(char *str) {
838  char *rtn;
839  char *tmp;
840 
841  rtn = rtpg_strreplace(str, " ", "", NULL);
842 
843  tmp = rtpg_strreplace(rtn, "\n", "", NULL);
844  pfree(rtn);
845  rtn = rtpg_strreplace(tmp, "\t", "", NULL);
846  pfree(tmp);
847  tmp = rtpg_strreplace(rtn, "\f", "", NULL);
848  pfree(rtn);
849  rtn = rtpg_strreplace(tmp, "\r", "", NULL);
850  pfree(tmp);
851 
852  return rtn;
853 }
854 
855 static char*
856 rtpg_trim(const char *input) {
857  char *rtn;
858  char *ptr;
859  uint32_t offset = 0;
860  int inputlen = 0;
861 
862  if (!input)
863  return NULL;
864  else if (!*input)
865  return (char *) input;
866 
867  /* trim left */
868  while (isspace(*input) && *input != '\0')
869  input++;
870 
871  /* trim right */
872  inputlen = strlen(input);
873  if (inputlen) {
874  ptr = ((char *) input) + inputlen;
875  while (isspace(*--ptr))
876  offset++;
877  }
878 
879  rtn = palloc(sizeof(char) * (inputlen - offset + 1));
880  if (rtn == NULL) {
881  fprintf(stderr, "Not enough memory\n");
882  return NULL;
883  }
884  strncpy(rtn, input, inputlen - offset);
885  rtn[inputlen - offset] = '\0';
886 
887  return rtn;
888 }
889 
890 /*
891 * reverse string search function from
892 * http://stackoverflow.com/a/1634398
893 */
894 static char *
895 rtpg_strrstr(const char *s1, const char *s2) {
896  int s1len = strlen(s1);
897  int s2len = strlen(s2);
898  char *s;
899 
900  if (s2len > s1len)
901  return NULL;
902 
903  s = (char *) (s1 + s1len - s2len);
904  for (; s >= s1; --s)
905  if (strncmp(s, s2, s2len) == 0)
906  return s;
907 
908  return NULL;
909 }
910 
911 static char*
912 rtpg_getSR(int srid)
913 {
914  int i = 0;
915  int len = 0;
916  char *sql = NULL;
917  int spi_result;
918  TupleDesc tupdesc;
919  SPITupleTable *tuptable = NULL;
920  HeapTuple tuple;
921  char *tmp = NULL;
922  char *srs = NULL;
923 
924 /*
925 SELECT
926  CASE
927  WHEN (upper(auth_name) = 'EPSG' OR upper(auth_name) = 'EPSGA') AND length(COALESCE(auth_srid::text, '')) > 0
928  THEN upper(auth_name) || ':' || auth_srid
929  WHEN length(COALESCE(auth_name, '') || COALESCE(auth_srid::text, '')) > 0
930  THEN COALESCE(auth_name, '') || COALESCE(auth_srid::text, '')
931  ELSE ''
932  END,
933  proj4text,
934  srtext
935 FROM spatial_ref_sys
936 WHERE srid = X
937 LIMIT 1
938 */
939 
940  len = sizeof(char) * (strlen("SELECT CASE WHEN (upper(auth_name) = 'EPSG' OR upper(auth_name) = 'EPSGA') AND length(COALESCE(auth_srid::text, '')) > 0 THEN upper(auth_name) || ':' || auth_srid WHEN length(COALESCE(auth_name, '') || COALESCE(auth_srid::text, '')) > 0 THEN COALESCE(auth_name, '') || COALESCE(auth_srid::text, '') ELSE '' END, proj4text, srtext FROM spatial_ref_sys WHERE srid = LIMIT 1") + MAX_INT_CHARLEN + 1);
941  sql = (char *) palloc(len);
942  if (NULL == sql) {
943  elog(ERROR, "rtpg_getSR: Could not allocate memory for sql\n");
944  return NULL;
945  }
946 
947  spi_result = SPI_connect();
948  if (spi_result != SPI_OK_CONNECT) {
949  pfree(sql);
950  elog(ERROR, "rtpg_getSR: Could not connect to database using SPI\n");
951  return NULL;
952  }
953 
954  /* execute query */
955  snprintf(sql, len, "SELECT CASE WHEN (upper(auth_name) = 'EPSG' OR upper(auth_name) = 'EPSGA') AND length(COALESCE(auth_srid::text, '')) > 0 THEN upper(auth_name) || ':' || auth_srid WHEN length(COALESCE(auth_name, '') || COALESCE(auth_srid::text, '')) > 0 THEN COALESCE(auth_name, '') || COALESCE(auth_srid::text, '') ELSE '' END, proj4text, srtext FROM spatial_ref_sys WHERE srid = %d LIMIT 1", srid);
956  POSTGIS_RT_DEBUGF(4, "SRS query: %s", sql);
957  spi_result = SPI_execute(sql, TRUE, 0);
958  SPI_pfree(sql);
959  if (spi_result != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
960  if (SPI_tuptable) SPI_freetuptable(tuptable);
961  SPI_finish();
962  elog(ERROR, "rtpg_getSR: Cannot find SRID (%d) in spatial_ref_sys", srid);
963  return NULL;
964  }
965 
966  tupdesc = SPI_tuptable->tupdesc;
967  tuptable = SPI_tuptable;
968  tuple = tuptable->vals[0];
969 
970  /* which column to use? */
971  for (i = 1; i < 4; i++) {
972  tmp = SPI_getvalue(tuple, tupdesc, i);
973 
974  /* value AND GDAL supports this SR */
975  if (
976  SPI_result != SPI_ERROR_NOATTRIBUTE &&
977  SPI_result != SPI_ERROR_NOOUTFUNC &&
978  tmp != NULL &&
979  strlen(tmp) &&
981  ) {
982  POSTGIS_RT_DEBUGF(4, "Value for column %d is %s", i, tmp);
983 
984  len = strlen(tmp) + 1;
985  srs = SPI_palloc(sizeof(char) * len);
986  if (NULL == srs) {
987  pfree(tmp);
988  if (SPI_tuptable) SPI_freetuptable(tuptable);
989  SPI_finish();
990  elog(ERROR, "rtpg_getSR: Could not allocate memory for spatial reference text\n");
991  return NULL;
992  }
993  strncpy(srs, tmp, len);
994  pfree(tmp);
995 
996  break;
997  }
998 
999  if (tmp != NULL)
1000  pfree(tmp);
1001  continue;
1002  }
1003 
1004  if (SPI_tuptable) SPI_freetuptable(tuptable);
1005  SPI_finish();
1006 
1007  /* unable to get SR info */
1008  if (srs == NULL) {
1009  if (SPI_tuptable) SPI_freetuptable(tuptable);
1010  SPI_finish();
1011  elog(ERROR, "rtpg_getSR: Could not find a viable spatial reference for SRID (%d)", srid);
1012  return NULL;
1013  }
1014 
1015  return srs;
1016 }
1017 
1019 Datum RASTER_lib_version(PG_FUNCTION_ARGS)
1020 {
1021  char ver[64];
1022  text *result;
1023 
1024  snprintf(ver, 64, "%s r%d", POSTGIS_LIB_VERSION, POSTGIS_SVN_REVISION);
1025  ver[63] = '\0';
1026 
1027  result = cstring2text(ver);
1028  PG_RETURN_TEXT_P(result);
1029 }
1030 
1032 Datum RASTER_lib_build_date(PG_FUNCTION_ARGS)
1033 {
1034  char *ver = POSTGIS_BUILD_DATE;
1035  text *result;
1036  result = palloc(VARHDRSZ + strlen(ver));
1037  SET_VARSIZE(result, VARHDRSZ + strlen(ver));
1038  memcpy(VARDATA(result), ver, strlen(ver));
1039  PG_RETURN_POINTER(result);
1040 }
1041 
1043 Datum RASTER_gdal_version(PG_FUNCTION_ARGS)
1044 {
1045  const char *ver = rt_util_gdal_version("--version");
1046  text *result;
1047 
1048  /* add indicator if GDAL isn't configured right */
1049  if (!rt_util_gdal_configured()) {
1050  char *rtn = NULL;
1051  rtn = palloc(strlen(ver) + strlen(" GDAL_DATA not found") + 1);
1052  if (!rtn)
1053  result = cstring2text(ver);
1054  else {
1055  sprintf(rtn, "%s GDAL_DATA not found", ver);
1056  result = cstring2text(rtn);
1057  pfree(rtn);
1058  }
1059  }
1060  else
1061  result = cstring2text(ver);
1062 
1063  PG_RETURN_POINTER(result);
1064 }
1065 
1067 Datum RASTER_minPossibleValue(PG_FUNCTION_ARGS)
1068 {
1069  text *pixeltypetext = NULL;
1070  char *pixeltypechar = NULL;
1071  rt_pixtype pixtype = PT_END;
1072  double pixsize = 0;
1073 
1074  if (PG_ARGISNULL(0))
1075  PG_RETURN_NULL();
1076 
1077  pixeltypetext = PG_GETARG_TEXT_P(0);
1078  pixeltypechar = text_to_cstring(pixeltypetext);
1079 
1080  pixtype = rt_pixtype_index_from_name(pixeltypechar);
1081  if (pixtype == PT_END) {
1082  elog(ERROR, "RASTER_minPossibleValue: Invalid pixel type: %s", pixeltypechar);
1083  PG_RETURN_NULL();
1084  }
1085 
1086  pixsize = rt_pixtype_get_min_value(pixtype);
1087 
1088  /*
1089  correct pixsize of unsigned pixel types
1090  example: for PT_8BUI, the value is CHAR_MIN but if char is signed,
1091  the value returned is -127 instead of 0.
1092  */
1093  switch (pixtype) {
1094  case PT_1BB:
1095  case PT_2BUI:
1096  case PT_4BUI:
1097  case PT_8BUI:
1098  case PT_16BUI:
1099  case PT_32BUI:
1100  pixsize = 0;
1101  break;
1102  default:
1103  break;
1104  }
1105 
1106  PG_RETURN_FLOAT8(pixsize);
1107 }
1108 
1114 Datum RASTER_in(PG_FUNCTION_ARGS)
1115 {
1116  rt_raster raster;
1117  char *hexwkb = PG_GETARG_CSTRING(0);
1118  void *result = NULL;
1119 
1120  POSTGIS_RT_DEBUG(3, "Starting");
1121 
1122  raster = rt_raster_from_hexwkb(hexwkb, strlen(hexwkb));
1123  if (raster == NULL)
1124  PG_RETURN_NULL();
1125 
1126  result = rt_raster_serialize(raster);
1127  rt_raster_destroy(raster);
1128  if (result == NULL)
1129  PG_RETURN_NULL();
1130 
1131  SET_VARSIZE(result, ((rt_pgraster*)result)->size);
1132  PG_RETURN_POINTER(result);
1133 }
1134 
1139 Datum RASTER_out(PG_FUNCTION_ARGS)
1140 {
1141  rt_pgraster *pgraster = NULL;
1142  rt_raster raster = NULL;
1143  uint32_t hexwkbsize = 0;
1144  char *hexwkb = NULL;
1145 
1146  POSTGIS_RT_DEBUG(3, "Starting");
1147 
1148  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
1149  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1150 
1151  raster = rt_raster_deserialize(pgraster, FALSE);
1152  if (!raster) {
1153  PG_FREE_IF_COPY(pgraster, 0);
1154  elog(ERROR, "RASTER_out: Could not deserialize raster");
1155  PG_RETURN_NULL();
1156  }
1157 
1158  hexwkb = rt_raster_to_hexwkb(raster, FALSE, &hexwkbsize);
1159  if (!hexwkb) {
1160  rt_raster_destroy(raster);
1161  PG_FREE_IF_COPY(pgraster, 0);
1162  elog(ERROR, "RASTER_out: Could not HEX-WKBize raster");
1163  PG_RETURN_NULL();
1164  }
1165 
1166  /* Free the raster objects used */
1167  rt_raster_destroy(raster);
1168  PG_FREE_IF_COPY(pgraster, 0);
1169 
1170  PG_RETURN_CSTRING(hexwkb);
1171 }
1172 
1177 Datum RASTER_to_bytea(PG_FUNCTION_ARGS)
1178 {
1179  rt_pgraster *pgraster = NULL;
1180  rt_raster raster = NULL;
1181  uint8_t *wkb = NULL;
1182  uint32_t wkb_size = 0;
1183  bytea *result = NULL;
1184  int result_size = 0;
1185 
1186  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
1187  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1188 
1189  /* Get raster object */
1190  raster = rt_raster_deserialize(pgraster, FALSE);
1191  if (!raster) {
1192  PG_FREE_IF_COPY(pgraster, 0);
1193  elog(ERROR, "RASTER_to_bytea: Could not deserialize raster");
1194  PG_RETURN_NULL();
1195  }
1196 
1197  /* Parse raster to wkb object */
1198  wkb = rt_raster_to_wkb(raster, FALSE, &wkb_size);
1199  if (!wkb) {
1200  rt_raster_destroy(raster);
1201  PG_FREE_IF_COPY(pgraster, 0);
1202  elog(ERROR, "RASTER_to_bytea: Could not allocate and generate WKB data");
1203  PG_RETURN_NULL();
1204  }
1205 
1206  /* Create varlena object */
1207  result_size = wkb_size + VARHDRSZ;
1208  result = (bytea *)palloc(result_size);
1209  SET_VARSIZE(result, result_size);
1210  memcpy(VARDATA(result), wkb, VARSIZE(result) - VARHDRSZ);
1211 
1212  /* Free raster objects used */
1213  rt_raster_destroy(raster);
1214  pfree(wkb);
1215  PG_FREE_IF_COPY(pgraster, 0);
1216 
1217  PG_RETURN_POINTER(result);
1218 }
1219 
1224 Datum RASTER_to_binary(PG_FUNCTION_ARGS)
1225 {
1226  rt_pgraster *pgraster = NULL;
1227  rt_raster raster = NULL;
1228  uint8_t *wkb = NULL;
1229  uint32_t wkb_size = 0;
1230  char *result = NULL;
1231  int result_size = 0;
1232  int outasin = FALSE;
1233 
1234  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
1235  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1236 
1237  /* Get raster object */
1238  raster = rt_raster_deserialize(pgraster, FALSE);
1239  if (!raster) {
1240  PG_FREE_IF_COPY(pgraster, 0);
1241  elog(ERROR, "RASTER_to_binary: Could not deserialize raster");
1242  PG_RETURN_NULL();
1243  }
1244 
1245  if (!PG_ARGISNULL(1))
1246  outasin = PG_GETARG_BOOL(1);
1247 
1248  /* Parse raster to wkb object */
1249  wkb = rt_raster_to_wkb(raster, outasin, &wkb_size);
1250  if (!wkb) {
1251  rt_raster_destroy(raster);
1252  PG_FREE_IF_COPY(pgraster, 0);
1253  elog(ERROR, "RASTER_to_binary: Could not allocate and generate WKB data");
1254  PG_RETURN_NULL();
1255  }
1256 
1257  /* Create varlena object */
1258  result_size = wkb_size + VARHDRSZ;
1259  result = (char *)palloc(result_size);
1260  SET_VARSIZE(result, result_size);
1261  memcpy(VARDATA(result), wkb, VARSIZE(result) - VARHDRSZ);
1262 
1263  /* Free raster objects used */
1264  rt_raster_destroy(raster);
1265  pfree(wkb);
1266  PG_FREE_IF_COPY(pgraster, 0);
1267 
1268  PG_RETURN_POINTER(result);
1269 }
1270 
1275 Datum RASTER_convex_hull(PG_FUNCTION_ARGS)
1276 {
1277  rt_pgraster *pgraster;
1278  rt_raster raster;
1279  LWGEOM *geom = NULL;
1280  GSERIALIZED* gser = NULL;
1281  size_t gser_size;
1282  int err = ES_NONE;
1283 
1284  bool minhull = FALSE;
1285 
1286  if (PG_ARGISNULL(0))
1287  PG_RETURN_NULL();
1288 
1289  /* # of args */
1290  if (PG_NARGS() > 1)
1291  minhull = TRUE;
1292 
1293  if (!minhull) {
1294  pgraster = (rt_pgraster *) PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(0), 0, sizeof(struct rt_raster_serialized_t));
1295  raster = rt_raster_deserialize(pgraster, TRUE);
1296  }
1297  else {
1298  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1299  raster = rt_raster_deserialize(pgraster, FALSE);
1300  }
1301 
1302  if (!raster) {
1303  PG_FREE_IF_COPY(pgraster, 0);
1304  elog(ERROR, "RASTER_convex_hull: Could not deserialize raster");
1305  PG_RETURN_NULL();
1306  }
1307 
1308  if (!minhull)
1309  err = rt_raster_get_convex_hull(raster, &geom);
1310  else {
1311  int nband = -1;
1312 
1313  /* get arg 1 */
1314  if (!PG_ARGISNULL(1)) {
1315  nband = PG_GETARG_INT32(1);
1316  if (!rt_raster_has_band(raster, nband - 1)) {
1317  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
1318  rt_raster_destroy(raster);
1319  PG_FREE_IF_COPY(pgraster, 0);
1320  PG_RETURN_NULL();
1321  }
1322  nband = nband - 1;
1323  }
1324 
1325  err = rt_raster_get_perimeter(raster, nband, &geom);
1326  }
1327 
1328  rt_raster_destroy(raster);
1329  PG_FREE_IF_COPY(pgraster, 0);
1330 
1331  if (err != ES_NONE) {
1332  elog(ERROR, "RASTER_convex_hull: Could not get raster's convex hull");
1333  PG_RETURN_NULL();
1334  }
1335  else if (geom == NULL) {
1336  elog(NOTICE, "Raster's convex hull is NULL");
1337  PG_RETURN_NULL();
1338  }
1339 
1340  gser = gserialized_from_lwgeom(geom, 0, &gser_size);
1341  lwgeom_free(geom);
1342 
1343  SET_VARSIZE(gser, gser_size);
1344  PG_RETURN_POINTER(gser);
1345 }
1346 
1348 Datum RASTER_dumpAsPolygons(PG_FUNCTION_ARGS) {
1349  FuncCallContext *funcctx;
1350  TupleDesc tupdesc;
1351  rt_geomval geomval;
1352  rt_geomval geomval2;
1353  int call_cntr;
1354  int max_calls;
1355 
1356  /* stuff done only on the first call of the function */
1357  if (SRF_IS_FIRSTCALL()) {
1358  MemoryContext oldcontext;
1359  int numbands;
1360  rt_pgraster *pgraster = NULL;
1361  rt_raster raster = NULL;
1362  int nband;
1363  bool exclude_nodata_value = TRUE;
1364  int nElements;
1365 
1366  POSTGIS_RT_DEBUG(2, "RASTER_dumpAsPolygons first call");
1367 
1368  /* create a function context for cross-call persistence */
1369  funcctx = SRF_FIRSTCALL_INIT();
1370 
1371  /* switch to memory context appropriate for multiple function calls */
1372  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
1373 
1374  /* Get input arguments */
1375  if (PG_ARGISNULL(0)) {
1376  MemoryContextSwitchTo(oldcontext);
1377  SRF_RETURN_DONE(funcctx);
1378  }
1379  pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1380  raster = rt_raster_deserialize(pgraster, FALSE);
1381  if (!raster) {
1382  PG_FREE_IF_COPY(pgraster, 0);
1383  ereport(ERROR, (
1384  errcode(ERRCODE_OUT_OF_MEMORY),
1385  errmsg("Could not deserialize raster")
1386  ));
1387  MemoryContextSwitchTo(oldcontext);
1388  SRF_RETURN_DONE(funcctx);
1389  }
1390 
1391  if (!PG_ARGISNULL(1))
1392  nband = PG_GETARG_UINT32(1);
1393  else
1394  nband = 1; /* By default, first band */
1395 
1396  POSTGIS_RT_DEBUGF(3, "band %d", nband);
1397  numbands = rt_raster_get_num_bands(raster);
1398 
1399  if (nband < 1 || nband > numbands) {
1400  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
1401  rt_raster_destroy(raster);
1402  PG_FREE_IF_COPY(pgraster, 0);
1403  MemoryContextSwitchTo(oldcontext);
1404  SRF_RETURN_DONE(funcctx);
1405  }
1406 
1407  if (!PG_ARGISNULL(2))
1408  exclude_nodata_value = PG_GETARG_BOOL(2);
1409 
1410  /* see if band is NODATA */
1411  if (rt_band_get_isnodata_flag(rt_raster_get_band(raster, nband - 1))) {
1412  POSTGIS_RT_DEBUGF(3, "Band at index %d is NODATA. Returning NULL", nband);
1413  rt_raster_destroy(raster);
1414  PG_FREE_IF_COPY(pgraster, 0);
1415  MemoryContextSwitchTo(oldcontext);
1416  SRF_RETURN_DONE(funcctx);
1417  }
1418 
1419  /* Polygonize raster */
1420 
1424  geomval = rt_raster_gdal_polygonize(raster, nband - 1, exclude_nodata_value, &nElements);
1425  rt_raster_destroy(raster);
1426  PG_FREE_IF_COPY(pgraster, 0);
1427  if (NULL == geomval) {
1428  ereport(ERROR, (
1429  errcode(ERRCODE_NO_DATA_FOUND),
1430  errmsg("Could not polygonize raster")
1431  ));
1432  MemoryContextSwitchTo(oldcontext);
1433  SRF_RETURN_DONE(funcctx);
1434  }
1435 
1436  POSTGIS_RT_DEBUGF(3, "raster dump, %d elements returned", nElements);
1437 
1438  /* Store needed information */
1439  funcctx->user_fctx = geomval;
1440 
1441  /* total number of tuples to be returned */
1442  funcctx->max_calls = nElements;
1443 
1444  /* Build a tuple descriptor for our result type */
1445  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
1446  ereport(ERROR, (
1447  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1448  errmsg("function returning record called in context that cannot accept type record")
1449  ));
1450  }
1451 
1452  BlessTupleDesc(tupdesc);
1453  funcctx->tuple_desc = tupdesc;
1454 
1455  MemoryContextSwitchTo(oldcontext);
1456  }
1457 
1458  /* stuff done on every call of the function */
1459  funcctx = SRF_PERCALL_SETUP();
1460 
1461  call_cntr = funcctx->call_cntr;
1462  max_calls = funcctx->max_calls;
1463  tupdesc = funcctx->tuple_desc;
1464  geomval2 = funcctx->user_fctx;
1465 
1466  /* do when there is more left to send */
1467  if (call_cntr < max_calls) {
1468  int values_length = 2;
1469  Datum values[values_length];
1470  bool nulls[values_length];
1471  HeapTuple tuple;
1472  Datum result;
1473 
1474  GSERIALIZED *gser = NULL;
1475  size_t gser_size = 0;
1476 
1477  POSTGIS_RT_DEBUGF(3, "call number %d", call_cntr);
1478 
1479  memset(nulls, FALSE, sizeof(bool) * values_length);
1480 
1481  /* convert LWGEOM to GSERIALIZED */
1482  gser = gserialized_from_lwgeom(lwpoly_as_lwgeom(geomval2[call_cntr].geom), 0, &gser_size);
1483  lwgeom_free(lwpoly_as_lwgeom(geomval2[call_cntr].geom));
1484 
1485  values[0] = PointerGetDatum(gser);
1486  values[1] = Float8GetDatum(geomval2[call_cntr].val);
1487 
1488  /* build a tuple */
1489  tuple = heap_form_tuple(tupdesc, values, nulls);
1490 
1491  /* make the tuple into a datum */
1492  result = HeapTupleGetDatum(tuple);
1493 
1494  SRF_RETURN_NEXT(funcctx, result);
1495  }
1496  /* do when there is no more left */
1497  else {
1498  pfree(geomval2);
1499  SRF_RETURN_DONE(funcctx);
1500  }
1501 }
1502 
1507 Datum RASTER_makeEmpty(PG_FUNCTION_ARGS)
1508 {
1509  uint16 width = 0, height = 0;
1510  double ipx = 0, ipy = 0, scalex = 0, scaley = 0, skewx = 0, skewy = 0;
1511  int32_t srid = SRID_UNKNOWN;
1512  rt_pgraster *pgraster = NULL;
1513  rt_raster raster;
1514 
1515  if (PG_NARGS() < 9) {
1516  elog(ERROR, "RASTER_makeEmpty: ST_MakeEmptyRaster requires 9 args");
1517  PG_RETURN_NULL();
1518  }
1519 
1520  if (!PG_ARGISNULL(0))
1521  width = PG_GETARG_UINT16(0);
1522 
1523  if (!PG_ARGISNULL(1))
1524  height = PG_GETARG_UINT16(1);
1525 
1526  if (!PG_ARGISNULL(2))
1527  ipx = PG_GETARG_FLOAT8(2);
1528 
1529  if (!PG_ARGISNULL(3))
1530  ipy = PG_GETARG_FLOAT8(3);
1531 
1532  if (!PG_ARGISNULL(4))
1533  scalex = PG_GETARG_FLOAT8(4);
1534 
1535  if (!PG_ARGISNULL(5))
1536  scaley = PG_GETARG_FLOAT8(5);
1537 
1538  if (!PG_ARGISNULL(6))
1539  skewx = PG_GETARG_FLOAT8(6);
1540 
1541  if (!PG_ARGISNULL(7))
1542  skewy = PG_GETARG_FLOAT8(7);
1543 
1544  if (!PG_ARGISNULL(8))
1545  srid = PG_GETARG_INT32(8);
1546 
1547  POSTGIS_RT_DEBUGF(4, "%dx%d, ip:%g,%g, scale:%g,%g, skew:%g,%g srid:%d",
1548  width, height, ipx, ipy, scalex, scaley,
1549  skewx, skewy, srid);
1550 
1551  raster = rt_raster_new(width, height);
1552  if (raster == NULL)
1553  PG_RETURN_NULL(); /* error was supposedly printed already */
1554 
1555  rt_raster_set_scale(raster, scalex, scaley);
1556  rt_raster_set_offsets(raster, ipx, ipy);
1557  rt_raster_set_skews(raster, skewx, skewy);
1558  rt_raster_set_srid(raster, srid);
1559 
1560  pgraster = rt_raster_serialize(raster);
1561  rt_raster_destroy(raster);
1562  if (!pgraster)
1563  PG_RETURN_NULL();
1564 
1565  SET_VARSIZE(pgraster, pgraster->size);
1566  PG_RETURN_POINTER(pgraster);
1567 }
1568 
1573 Datum RASTER_getSRID(PG_FUNCTION_ARGS)
1574 {
1575  rt_pgraster *pgraster;
1576  rt_raster raster;
1577  int32_t srid;
1578 
1579  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
1580  pgraster = (rt_pgraster *) PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(0), 0, sizeof(struct rt_raster_serialized_t));
1581 
1582  raster = rt_raster_deserialize(pgraster, TRUE);
1583  if ( ! raster ) {
1584  PG_FREE_IF_COPY(pgraster, 0);
1585  elog(ERROR, "RASTER_getSRID: Could not deserialize raster");
1586  PG_RETURN_NULL();
1587  }
1588 
1589  srid = rt_raster_get_srid(raster);
1590 
1591  rt_raster_destroy(raster);
1592  PG_FREE_IF_COPY(pgraster, 0);
1593 
1594  PG_RETURN_INT32(srid);
1595 }
1596 
1601 Datum RASTER_setSRID(PG_FUNCTION_ARGS)
1602 {
1603  rt_pgraster *pgraster = NULL;
1604  rt_pgraster *pgrtn = NULL;
1605  rt_raster raster;
1606  int32_t newSRID = PG_GETARG_INT32(1);
1607 
1608  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
1609  pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1610 
1611  raster = rt_raster_deserialize(pgraster, FALSE);
1612  if (!raster) {
1613  PG_FREE_IF_COPY(pgraster, 0);
1614  elog(ERROR, "RASTER_setSRID: Could not deserialize raster");
1615  PG_RETURN_NULL();
1616  }
1617 
1618  rt_raster_set_srid(raster, newSRID);
1619 
1620  pgrtn = rt_raster_serialize(raster);
1621  rt_raster_destroy(raster);
1622  PG_FREE_IF_COPY(pgraster, 0);
1623  if (!pgrtn) PG_RETURN_NULL();
1624 
1625  SET_VARSIZE(pgrtn, pgrtn->size);
1626 
1627  PG_RETURN_POINTER(pgrtn);
1628 }
1629 
1634 Datum RASTER_getWidth(PG_FUNCTION_ARGS)
1635 {
1636  rt_pgraster *pgraster;
1637  rt_raster raster;
1638  uint16_t width;
1639 
1640  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
1641  pgraster = (rt_pgraster *) PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(0), 0, sizeof(struct rt_raster_serialized_t));
1642 
1643  raster = rt_raster_deserialize(pgraster, TRUE);
1644  if ( ! raster ) {
1645  PG_FREE_IF_COPY(pgraster, 0);
1646  elog(ERROR, "RASTER_getWidth: Could not deserialize raster");
1647  PG_RETURN_NULL();
1648  }
1649 
1650  width = rt_raster_get_width(raster);
1651 
1652  rt_raster_destroy(raster);
1653  PG_FREE_IF_COPY(pgraster, 0);
1654 
1655  PG_RETURN_INT32(width);
1656 }
1657 
1662 Datum RASTER_getHeight(PG_FUNCTION_ARGS)
1663 {
1664  rt_pgraster *pgraster;
1665  rt_raster raster;
1666  uint16_t height;
1667 
1668  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
1669  pgraster = (rt_pgraster *) PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(0), 0, sizeof(struct rt_raster_serialized_t));
1670 
1671  raster = rt_raster_deserialize(pgraster, TRUE);
1672  if ( ! raster ) {
1673  PG_FREE_IF_COPY(pgraster, 0);
1674  elog(ERROR, "RASTER_getHeight: Could not deserialize raster");
1675  PG_RETURN_NULL();
1676  }
1677 
1678  height = rt_raster_get_height(raster);
1679 
1680  rt_raster_destroy(raster);
1681  PG_FREE_IF_COPY(pgraster, 0);
1682 
1683  PG_RETURN_INT32(height);
1684 }
1685 
1690 Datum RASTER_getNumBands(PG_FUNCTION_ARGS)
1691 {
1692  rt_pgraster *pgraster;
1693  rt_raster raster;
1694  int32_t num_bands;
1695 
1696  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
1697  pgraster = (rt_pgraster *) PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(0), 0, sizeof(struct rt_raster_serialized_t));
1698 
1699  raster = rt_raster_deserialize(pgraster, TRUE);
1700  if ( ! raster ) {
1701  PG_FREE_IF_COPY(pgraster, 0);
1702  elog(ERROR, "RASTER_getNumBands: Could not deserialize raster");
1703  PG_RETURN_NULL();
1704  }
1705 
1706  num_bands = rt_raster_get_num_bands(raster);
1707 
1708  rt_raster_destroy(raster);
1709  PG_FREE_IF_COPY(pgraster, 0);
1710 
1711  PG_RETURN_INT32(num_bands);
1712 }
1713 
1718 Datum RASTER_getXScale(PG_FUNCTION_ARGS)
1719 {
1720  rt_pgraster *pgraster;
1721  rt_raster raster;
1722  double xsize;
1723 
1724  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
1725  pgraster = (rt_pgraster *) PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(0), 0, sizeof(struct rt_raster_serialized_t));
1726 
1727  raster = rt_raster_deserialize(pgraster, TRUE);
1728  if ( ! raster ) {
1729  PG_FREE_IF_COPY(pgraster, 0);
1730  elog(ERROR, "RASTER_getXScale: Could not deserialize raster");
1731  PG_RETURN_NULL();
1732  }
1733 
1734  xsize = rt_raster_get_x_scale(raster);
1735 
1736  rt_raster_destroy(raster);
1737  PG_FREE_IF_COPY(pgraster, 0);
1738 
1739  PG_RETURN_FLOAT8(xsize);
1740 }
1741 
1746 Datum RASTER_getYScale(PG_FUNCTION_ARGS)
1747 {
1748  rt_pgraster *pgraster;
1749  rt_raster raster;
1750  double ysize;
1751 
1752  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
1753  pgraster = (rt_pgraster *) PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(0), 0, sizeof(struct rt_raster_serialized_t));
1754 
1755  raster = rt_raster_deserialize(pgraster, TRUE);
1756  if ( ! raster ) {
1757  PG_FREE_IF_COPY(pgraster, 0);
1758  elog(ERROR, "RASTER_getYScale: Could not deserialize raster");
1759  PG_RETURN_NULL();
1760  }
1761 
1762  ysize = rt_raster_get_y_scale(raster);
1763 
1764  rt_raster_destroy(raster);
1765  PG_FREE_IF_COPY(pgraster, 0);
1766 
1767  PG_RETURN_FLOAT8(ysize);
1768 }
1769 
1774 Datum RASTER_setScale(PG_FUNCTION_ARGS)
1775 {
1776  rt_pgraster *pgraster = NULL;
1777  rt_pgraster *pgrtn = NULL;
1778  rt_raster raster;
1779  double size = PG_GETARG_FLOAT8(1);
1780 
1781  if (PG_ARGISNULL(0))
1782  PG_RETURN_NULL();
1783  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1784  raster = rt_raster_deserialize(pgraster, FALSE);
1785  if (!raster) {
1786  PG_FREE_IF_COPY(pgraster, 0);
1787  elog(ERROR, "RASTER_setScale: Could not deserialize raster");
1788  PG_RETURN_NULL();
1789  }
1790 
1791  rt_raster_set_scale(raster, size, size);
1792 
1793  pgrtn = rt_raster_serialize(raster);
1794  rt_raster_destroy(raster);
1795  PG_FREE_IF_COPY(pgraster, 0);
1796  if (!pgrtn)
1797  PG_RETURN_NULL();
1798 
1799  SET_VARSIZE(pgrtn, pgrtn->size);
1800  PG_RETURN_POINTER(pgrtn);
1801 }
1802 
1807 Datum RASTER_setScaleXY(PG_FUNCTION_ARGS)
1808 {
1809  rt_pgraster *pgraster = NULL;
1810  rt_pgraster *pgrtn = NULL;
1811  rt_raster raster;
1812  double xscale = PG_GETARG_FLOAT8(1);
1813  double yscale = PG_GETARG_FLOAT8(2);
1814 
1815  if (PG_ARGISNULL(0))
1816  PG_RETURN_NULL();
1817  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1818  raster = rt_raster_deserialize(pgraster, FALSE);
1819  if (!raster) {
1820  PG_FREE_IF_COPY(pgraster, 0);
1821  elog(ERROR, "RASTER_setScaleXY: Could not deserialize raster");
1822  PG_RETURN_NULL();
1823  }
1824 
1825  rt_raster_set_scale(raster, xscale, yscale);
1826  pgrtn = rt_raster_serialize(raster);
1827  rt_raster_destroy(raster);
1828  PG_FREE_IF_COPY(pgraster, 0);
1829  if (!pgrtn)
1830  PG_RETURN_NULL();
1831 
1832  SET_VARSIZE(pgrtn, pgrtn->size);
1833  PG_RETURN_POINTER(pgrtn);
1834 }
1835 
1840 Datum RASTER_getXSkew(PG_FUNCTION_ARGS)
1841 {
1842  rt_pgraster *pgraster;
1843  rt_raster raster;
1844  double xskew;
1845 
1846  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
1847  pgraster = (rt_pgraster *) PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(0), 0, sizeof(struct rt_raster_serialized_t));
1848 
1849  raster = rt_raster_deserialize(pgraster, TRUE);
1850  if ( ! raster ) {
1851  PG_FREE_IF_COPY(pgraster, 0);
1852  elog(ERROR, "RASTER_getXSkew: Could not deserialize raster");
1853  PG_RETURN_NULL();
1854  }
1855 
1856  xskew = rt_raster_get_x_skew(raster);
1857 
1858  rt_raster_destroy(raster);
1859  PG_FREE_IF_COPY(pgraster, 0);
1860 
1861  PG_RETURN_FLOAT8(xskew);
1862 }
1863 
1868 Datum RASTER_getYSkew(PG_FUNCTION_ARGS)
1869 {
1870  rt_pgraster *pgraster;
1871  rt_raster raster;
1872  double yskew;
1873 
1874  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
1875  pgraster = (rt_pgraster *) PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(0), 0, sizeof(struct rt_raster_serialized_t));
1876 
1877  raster = rt_raster_deserialize(pgraster, TRUE);
1878  if ( ! raster ) {
1879  PG_FREE_IF_COPY(pgraster, 0);
1880  elog(ERROR, "RASTER_getYSkew: Could not deserialize raster");
1881  PG_RETURN_NULL();
1882  }
1883 
1884  yskew = rt_raster_get_y_skew(raster);
1885 
1886  rt_raster_destroy(raster);
1887  PG_FREE_IF_COPY(pgraster, 0);
1888 
1889  PG_RETURN_FLOAT8(yskew);
1890 }
1891 
1896 Datum RASTER_setSkew(PG_FUNCTION_ARGS)
1897 {
1898  rt_pgraster *pgraster = NULL;
1899  rt_pgraster *pgrtn = NULL;
1900  rt_raster raster;
1901  double skew = PG_GETARG_FLOAT8(1);
1902 
1903  if (PG_ARGISNULL(0))
1904  PG_RETURN_NULL();
1905  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1906  raster = rt_raster_deserialize(pgraster, FALSE);
1907  if (!raster) {
1908  PG_FREE_IF_COPY(pgraster, 0);
1909  elog(ERROR, "RASTER_setSkew: Could not deserialize raster");
1910  PG_RETURN_NULL();
1911  }
1912 
1913  rt_raster_set_skews(raster, skew, skew);
1914 
1915  pgrtn = rt_raster_serialize(raster);
1916  rt_raster_destroy(raster);
1917  PG_FREE_IF_COPY(pgraster, 0);
1918  if (!pgrtn)
1919  PG_RETURN_NULL();
1920 
1921  SET_VARSIZE(pgrtn, pgrtn->size);
1922  PG_RETURN_POINTER(pgrtn);
1923 }
1924 
1929 Datum RASTER_setSkewXY(PG_FUNCTION_ARGS)
1930 {
1931  rt_pgraster *pgraster = NULL;
1932  rt_pgraster *pgrtn = NULL;
1933  rt_raster raster;
1934  double xskew = PG_GETARG_FLOAT8(1);
1935  double yskew = PG_GETARG_FLOAT8(2);
1936 
1937  if (PG_ARGISNULL(0))
1938  PG_RETURN_NULL();
1939  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1940  raster = rt_raster_deserialize(pgraster, FALSE);
1941  if (!raster) {
1942  PG_FREE_IF_COPY(pgraster, 0);
1943  elog(ERROR, "RASTER_setSkewXY: Could not deserialize raster");
1944  PG_RETURN_NULL();
1945  }
1946 
1947  rt_raster_set_skews(raster, xskew, yskew);
1948 
1949  pgrtn = rt_raster_serialize(raster);
1950  rt_raster_destroy(raster);
1951  PG_FREE_IF_COPY(pgraster, 0);
1952  if (!pgrtn)
1953  PG_RETURN_NULL();
1954 
1955  SET_VARSIZE(pgrtn, pgrtn->size);
1956  PG_RETURN_POINTER(pgrtn);
1957 }
1958 
1963 Datum RASTER_getXUpperLeft(PG_FUNCTION_ARGS)
1964 {
1965  rt_pgraster *pgraster;
1966  rt_raster raster;
1967  double xul;
1968 
1969  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
1970  pgraster = (rt_pgraster *) PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(0), 0, sizeof(struct rt_raster_serialized_t));
1971 
1972  raster = rt_raster_deserialize(pgraster, TRUE);
1973  if ( ! raster ) {
1974  PG_FREE_IF_COPY(pgraster, 0);
1975  elog(ERROR, "RASTER_getXUpperLeft: Could not deserialize raster");
1976  PG_RETURN_NULL();
1977  }
1978 
1979  xul = rt_raster_get_x_offset(raster);
1980 
1981  rt_raster_destroy(raster);
1982  PG_FREE_IF_COPY(pgraster, 0);
1983 
1984  PG_RETURN_FLOAT8(xul);
1985 }
1986 
1991 Datum RASTER_getYUpperLeft(PG_FUNCTION_ARGS)
1992 {
1993  rt_pgraster *pgraster;
1994  rt_raster raster;
1995  double yul;
1996 
1997  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
1998  pgraster = (rt_pgraster *) PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(0), 0, sizeof(struct rt_raster_serialized_t));
1999 
2000  raster = rt_raster_deserialize(pgraster, TRUE);
2001  if ( ! raster ) {
2002  PG_FREE_IF_COPY(pgraster, 0);
2003  elog(ERROR, "RASTER_getYUpperLeft: Could not deserialize raster");
2004  PG_RETURN_NULL();
2005  }
2006 
2007  yul = rt_raster_get_y_offset(raster);
2008 
2009  rt_raster_destroy(raster);
2010  PG_FREE_IF_COPY(pgraster, 0);
2011 
2012  PG_RETURN_FLOAT8(yul);
2013 }
2014 
2019 Datum RASTER_setUpperLeftXY(PG_FUNCTION_ARGS)
2020 {
2021  rt_pgraster *pgraster = NULL;
2022  rt_pgraster *pgrtn = NULL;
2023  rt_raster raster;
2024  double xoffset = PG_GETARG_FLOAT8(1);
2025  double yoffset = PG_GETARG_FLOAT8(2);
2026 
2027  if (PG_ARGISNULL(0))
2028  PG_RETURN_NULL();
2029  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2030  raster = rt_raster_deserialize(pgraster, FALSE);
2031  if (!raster) {
2032  PG_FREE_IF_COPY(pgraster, 0);
2033  elog(ERROR, "RASTER_setUpperLeftXY: Could not deserialize raster");
2034  PG_RETURN_NULL();
2035  }
2036 
2037  rt_raster_set_offsets(raster, xoffset, yoffset);
2038 
2039  pgrtn = rt_raster_serialize(raster);
2040  rt_raster_destroy(raster);
2041  PG_FREE_IF_COPY(pgraster, 0);
2042  if (!pgrtn)
2043  PG_RETURN_NULL();
2044 
2045  SET_VARSIZE(pgrtn, pgrtn->size);
2046  PG_RETURN_POINTER(pgrtn);
2047 }
2048 
2057 Datum RASTER_getPixelWidth(PG_FUNCTION_ARGS)
2058 {
2059  rt_pgraster *pgraster;
2060  rt_raster raster;
2061  double xscale;
2062  double yskew;
2063  double pwidth;
2064 
2065  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
2066  pgraster = (rt_pgraster *)PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(0), 0, sizeof(struct rt_raster_serialized_t));
2067 
2068  raster = rt_raster_deserialize(pgraster, TRUE);
2069  if (!raster) {
2070  PG_FREE_IF_COPY(pgraster, 0);
2071  elog(ERROR, "RASTER_getPixelWidth: Could not deserialize raster");
2072  PG_RETURN_NULL();
2073  }
2074 
2075  xscale = rt_raster_get_x_scale(raster);
2076  yskew = rt_raster_get_y_skew(raster);
2077  pwidth = sqrt(xscale*xscale + yskew*yskew);
2078 
2079  rt_raster_destroy(raster);
2080  PG_FREE_IF_COPY(pgraster, 0);
2081 
2082  PG_RETURN_FLOAT8(pwidth);
2083 }
2084 
2093 Datum RASTER_getPixelHeight(PG_FUNCTION_ARGS)
2094 {
2095  rt_pgraster *pgraster;
2096  rt_raster raster;
2097  double yscale;
2098  double xskew;
2099  double pheight;
2100 
2101  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
2102  pgraster = (rt_pgraster *)PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(0), 0, sizeof(struct rt_raster_serialized_t));
2103 
2104  raster = rt_raster_deserialize(pgraster, TRUE);
2105  if (!raster) {
2106  PG_FREE_IF_COPY(pgraster, 0);
2107  elog(ERROR, "RASTER_getPixelHeight: Could not deserialize raster");
2108  PG_RETURN_NULL();
2109  }
2110 
2111  yscale = rt_raster_get_y_scale(raster);
2112  xskew = rt_raster_get_x_skew(raster);
2113  pheight = sqrt(yscale*yscale + xskew*xskew);
2114 
2115  rt_raster_destroy(raster);
2116  PG_FREE_IF_COPY(pgraster, 0);
2117 
2118  PG_RETURN_FLOAT8(pheight);
2119 }
2120 
2125 Datum RASTER_setGeotransform(PG_FUNCTION_ARGS)
2126 {
2127  rt_pgraster *pgraster = NULL;
2128  rt_pgraster *pgrtn = NULL;
2129  rt_raster raster;
2130  float8 imag, jmag, theta_i, theta_ij, xoffset, yoffset;
2131 
2132  if (
2133  PG_ARGISNULL(0) || PG_ARGISNULL(1) || PG_ARGISNULL(2) ||
2134  PG_ARGISNULL(3) || PG_ARGISNULL(4) ||
2135  PG_ARGISNULL(5) || PG_ARGISNULL(6)
2136  ) {
2137  PG_RETURN_NULL();
2138  }
2139 
2140  /* get the inputs */
2141  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2142  imag = PG_GETARG_FLOAT8(1);
2143  jmag = PG_GETARG_FLOAT8(2);
2144  theta_i = PG_GETARG_FLOAT8(3);
2145  theta_ij = PG_GETARG_FLOAT8(4);
2146  xoffset = PG_GETARG_FLOAT8(5);
2147  yoffset = PG_GETARG_FLOAT8(6);
2148 
2149  raster = rt_raster_deserialize(pgraster, TRUE);
2150  if (!raster) {
2151  PG_FREE_IF_COPY(pgraster, 0);
2152  elog(ERROR, "RASTER_setGeotransform: Could not deserialize raster");
2153  PG_RETURN_NULL();
2154  }
2155 
2156  /* store the new geotransform */
2157  rt_raster_set_phys_params(raster, imag,jmag,theta_i,theta_ij);
2158  rt_raster_set_offsets(raster, xoffset, yoffset);
2159 
2160  /* prep the return value */
2161  pgrtn = rt_raster_serialize(raster);
2162  rt_raster_destroy(raster);
2163  PG_FREE_IF_COPY(pgraster, 0);
2164  if (!pgrtn)
2165  PG_RETURN_NULL();
2166 
2167  SET_VARSIZE(pgrtn, pgrtn->size);
2168  PG_RETURN_POINTER(pgrtn);
2169 }
2170 
2176 Datum RASTER_getGeotransform(PG_FUNCTION_ARGS)
2177 {
2178  rt_pgraster *pgraster = NULL;
2179  rt_raster raster = NULL;
2180 
2181  double imag;
2182  double jmag;
2183  double theta_i;
2184  double theta_ij;
2185  /*
2186  double xoffset;
2187  double yoffset;
2188  */
2189 
2190  TupleDesc result_tuple; /* for returning a composite */
2191  int values_length = 6;
2192  Datum values[values_length];
2193  bool nulls[values_length];
2194  HeapTuple heap_tuple ; /* instance of the tuple to return */
2195  Datum result;
2196 
2197  POSTGIS_RT_DEBUG(3, "RASTER_getGeotransform: Starting");
2198 
2199  /* get argument */
2200  if (PG_ARGISNULL(0))
2201  PG_RETURN_NULL();
2202  pgraster = (rt_pgraster *)PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(0), 0, sizeof(struct rt_raster_serialized_t));
2203 
2204  /* raster */
2205  raster = rt_raster_deserialize(pgraster, TRUE);
2206  if (!raster) {
2207  PG_FREE_IF_COPY(pgraster, 0);
2208  elog(ERROR, "RASTER_getGeotransform: Could not deserialize raster");
2209  PG_RETURN_NULL();
2210  }
2211 
2212  /* do the calculation */
2214  rt_raster_get_x_scale(raster),
2215  rt_raster_get_x_skew(raster),
2216  rt_raster_get_y_skew(raster),
2217  rt_raster_get_y_scale(raster),
2218  &imag, &jmag, &theta_i, &theta_ij) ;
2219 
2220  rt_raster_destroy(raster);
2221  PG_FREE_IF_COPY(pgraster, 0);
2222 
2223  /* setup the return value infrastructure */
2224  if (get_call_result_type(fcinfo, NULL, &result_tuple) != TYPEFUNC_COMPOSITE) {
2225  ereport(ERROR, (
2226  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
2227  errmsg("RASTER_getGeotransform(): function returning record called in context that cannot accept type record"
2228  )
2229  ));
2230  PG_RETURN_NULL();
2231  }
2232 
2233  BlessTupleDesc(result_tuple);
2234 
2235  /* get argument */
2236  /* prep the composite return value */
2237  /* construct datum array */
2238  values[0] = Float8GetDatum(imag);
2239  values[1] = Float8GetDatum(jmag);
2240  values[2] = Float8GetDatum(theta_i);
2241  values[3] = Float8GetDatum(theta_ij);
2242  values[4] = Float8GetDatum(rt_raster_get_x_offset(raster));
2243  values[5] = Float8GetDatum(rt_raster_get_y_offset(raster));
2244 
2245  memset(nulls, FALSE, sizeof(bool) * values_length);
2246 
2247  /* stick em on the heap */
2248  heap_tuple = heap_form_tuple(result_tuple, values, nulls);
2249 
2250  /* make the tuple into a datum */
2251  result = HeapTupleGetDatum(heap_tuple);
2252 
2253  PG_RETURN_DATUM(result);
2254 }
2255 
2256 
2269 Datum RASTER_setRotation(PG_FUNCTION_ARGS)
2270 {
2271  rt_pgraster *pgraster = NULL;
2272  rt_pgraster *pgrtn = NULL;
2273  rt_raster raster;
2274  double rotation = PG_GETARG_FLOAT8(1);
2275  double imag, jmag, theta_i, theta_ij;
2276 
2277  if (PG_ARGISNULL(0))
2278  PG_RETURN_NULL();
2279  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2280 
2281  raster = rt_raster_deserialize(pgraster, FALSE);
2282  if (! raster ) {
2283  PG_FREE_IF_COPY(pgraster, 0);
2284  elog(ERROR, "RASTER_setRotation: Could not deserialize raster");
2285  PG_RETURN_NULL();
2286  }
2287 
2288  /* preserve all defining characteristics of the grid except for rotation */
2289  rt_raster_get_phys_params(raster, &imag, &jmag, &theta_i, &theta_ij);
2290  rt_raster_set_phys_params(raster, imag, jmag, rotation, theta_ij);
2291 
2292  pgrtn = rt_raster_serialize(raster);
2293  rt_raster_destroy(raster);
2294  PG_FREE_IF_COPY(pgraster, 0);
2295  if (!pgrtn)
2296  PG_RETURN_NULL();
2297 
2298  SET_VARSIZE(pgrtn, pgrtn->size);
2299  PG_RETURN_POINTER(pgrtn);
2300 }
2301 
2307 Datum RASTER_getBandPixelType(PG_FUNCTION_ARGS)
2308 {
2309  rt_pgraster *pgraster = NULL;
2310  rt_raster raster = NULL;
2311  rt_band band = NULL;
2312  rt_pixtype pixtype;
2313  int32_t bandindex;
2314 
2315  /* Deserialize raster */
2316  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
2317  pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2318 
2319  /* Index is 1-based */
2320  bandindex = PG_GETARG_INT32(1);
2321  if ( bandindex < 1 ) {
2322  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
2323  PG_FREE_IF_COPY(pgraster, 0);
2324  PG_RETURN_NULL();
2325  }
2326 
2327  raster = rt_raster_deserialize(pgraster, FALSE);
2328  if ( ! raster ) {
2329  PG_FREE_IF_COPY(pgraster, 0);
2330  elog(ERROR, "RASTER_getBandPixelType: Could not deserialize raster");
2331  PG_RETURN_NULL();
2332  }
2333 
2334  /* Fetch requested band and its pixel type */
2335  band = rt_raster_get_band(raster, bandindex - 1);
2336  if ( ! band ) {
2337  elog(NOTICE, "Could not find raster band of index %d when getting pixel type. Returning NULL", bandindex);
2338  rt_raster_destroy(raster);
2339  PG_FREE_IF_COPY(pgraster, 0);
2340  PG_RETURN_NULL();
2341  }
2342 
2343  pixtype = rt_band_get_pixtype(band);
2344 
2345  rt_raster_destroy(raster);
2346  PG_FREE_IF_COPY(pgraster, 0);
2347 
2348  PG_RETURN_INT32(pixtype);
2349 }
2350 
2357 Datum RASTER_getBandPixelTypeName(PG_FUNCTION_ARGS)
2358 {
2359  rt_pgraster *pgraster = NULL;
2360  rt_raster raster = NULL;
2361  rt_band band = NULL;
2362  rt_pixtype pixtype;
2363  int32_t bandindex;
2364  const size_t name_size = 8; /* size of type name */
2365  size_t size = 0;
2366  char *ptr = NULL;
2367  text *result = NULL;
2368 
2369  /* Deserialize raster */
2370  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
2371  pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2372 
2373  /* Index is 1-based */
2374  bandindex = PG_GETARG_INT32(1);
2375  if ( bandindex < 1 ) {
2376  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
2377  PG_FREE_IF_COPY(pgraster, 0);
2378  PG_RETURN_NULL();
2379  }
2380 
2381  raster = rt_raster_deserialize(pgraster, FALSE);
2382  if ( ! raster ) {
2383  PG_FREE_IF_COPY(pgraster, 0);
2384  elog(ERROR, "RASTER_getBandPixelTypeName: Could not deserialize raster");
2385  PG_RETURN_NULL();
2386  }
2387 
2388  /* Fetch requested band and its pixel type */
2389  band = rt_raster_get_band(raster, bandindex - 1);
2390  if ( ! band ) {
2391  elog(NOTICE, "Could not find raster band of index %d when getting pixel type name. Returning NULL", bandindex);
2392  rt_raster_destroy(raster);
2393  PG_FREE_IF_COPY(pgraster, 0);
2394  PG_RETURN_NULL();
2395  }
2396 
2397  pixtype = rt_band_get_pixtype(band);
2398 
2399  result = palloc(VARHDRSZ + name_size);
2400  /* We don't need to check for NULL pointer, because if out of memory, palloc
2401  * exit via elog(ERROR). It never returns NULL.
2402  */
2403 
2404  memset(VARDATA(result), 0, name_size);
2405  ptr = (char *)result + VARHDRSZ;
2406  strcpy(ptr, rt_pixtype_name(pixtype));
2407 
2408  size = VARHDRSZ + strlen(ptr);
2409  SET_VARSIZE(result, size);
2410 
2411  rt_raster_destroy(raster);
2412  PG_FREE_IF_COPY(pgraster, 0);
2413 
2414  PG_RETURN_TEXT_P(result);
2415 }
2416 
2422 Datum RASTER_getBandNoDataValue(PG_FUNCTION_ARGS)
2423 {
2424  rt_pgraster *pgraster = NULL;
2425  rt_raster raster = NULL;
2426  rt_band band = NULL;
2427  int32_t bandindex;
2428  double nodata;
2429 
2430  /* Deserialize raster */
2431  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
2432  pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2433 
2434  /* Index is 1-based */
2435  bandindex = PG_GETARG_INT32(1);
2436  if ( bandindex < 1 ) {
2437  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
2438  PG_FREE_IF_COPY(pgraster, 0);
2439  PG_RETURN_NULL();
2440  }
2441 
2442  raster = rt_raster_deserialize(pgraster, FALSE);
2443  if ( ! raster ) {
2444  PG_FREE_IF_COPY(pgraster, 0);
2445  elog(ERROR, "RASTER_getBandNoDataValue: Could not deserialize raster");
2446  PG_RETURN_NULL();
2447  }
2448 
2449  /* Fetch requested band and its nodata value */
2450  band = rt_raster_get_band(raster, bandindex - 1);
2451  if ( ! band ) {
2452  elog(NOTICE, "Could not find raster band of index %d when getting band nodata value. Returning NULL", bandindex);
2453  rt_raster_destroy(raster);
2454  PG_FREE_IF_COPY(pgraster, 0);
2455  PG_RETURN_NULL();
2456  }
2457 
2458  if ( ! rt_band_get_hasnodata_flag(band) ) {
2459  /* Raster does not have a nodata value set so we return NULL */
2460  rt_raster_destroy(raster);
2461  PG_FREE_IF_COPY(pgraster, 0);
2462  PG_RETURN_NULL();
2463  }
2464 
2465  rt_band_get_nodata(band, &nodata);
2466 
2467  rt_raster_destroy(raster);
2468  PG_FREE_IF_COPY(pgraster, 0);
2469 
2470  PG_RETURN_FLOAT8(nodata);
2471 }
2472 
2473 
2478 Datum RASTER_setBandNoDataValue(PG_FUNCTION_ARGS)
2479 {
2480  rt_pgraster *pgraster = NULL;
2481  rt_pgraster *pgrtn = NULL;
2482  rt_raster raster = NULL;
2483  rt_band band = NULL;
2484  double nodata;
2485  int32_t bandindex;
2486  bool forcechecking = FALSE;
2487  bool skipset = FALSE;
2488 
2489  /* Deserialize raster */
2490  if (PG_ARGISNULL(0))
2491  PG_RETURN_NULL();
2492  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2493 
2494  /* Check index is not NULL or smaller than 1 */
2495  if (PG_ARGISNULL(1))
2496  bandindex = -1;
2497  else
2498  bandindex = PG_GETARG_INT32(1);
2499  if (bandindex < 1) {
2500  elog(NOTICE, "Invalid band index (must use 1-based). Nodata value not set. Returning original raster");
2501  skipset = TRUE;
2502  }
2503 
2504  raster = rt_raster_deserialize(pgraster, FALSE);
2505  if (!raster) {
2506  PG_FREE_IF_COPY(pgraster, 0);
2507  elog(ERROR, "RASTER_setBandNoDataValue: Could not deserialize raster");
2508  PG_RETURN_NULL();
2509  }
2510 
2511  if (!skipset) {
2512  /* Fetch requested band */
2513  band = rt_raster_get_band(raster, bandindex - 1);
2514  if (!band) {
2515  elog(NOTICE, "Could not find raster band of index %d when setting pixel value. Nodata value not set. Returning original raster", bandindex);
2516  }
2517  else {
2518  if (!PG_ARGISNULL(3))
2519  forcechecking = PG_GETARG_BOOL(3);
2520 
2521  if (PG_ARGISNULL(2)) {
2522  /* Set the hasnodata flag to FALSE */
2524  POSTGIS_RT_DEBUGF(3, "Raster band %d does not have a nodata value", bandindex);
2525  }
2526  else {
2527  /* Get the nodata value */
2528  nodata = PG_GETARG_FLOAT8(2);
2529 
2530  /* Set the band's nodata value */
2531  rt_band_set_nodata(band, nodata, NULL);
2532 
2533  /* Recheck all pixels if requested */
2534  if (forcechecking)
2536  }
2537  }
2538  }
2539 
2540  pgrtn = rt_raster_serialize(raster);
2541  rt_raster_destroy(raster);
2542  PG_FREE_IF_COPY(pgraster, 0);
2543  if (!pgrtn)
2544  PG_RETURN_NULL();
2545 
2546  SET_VARSIZE(pgrtn, pgrtn->size);
2547  PG_RETURN_POINTER(pgrtn);
2548 }
2549 
2551 Datum RASTER_setBandIsNoData(PG_FUNCTION_ARGS)
2552 {
2553  rt_pgraster *pgraster = NULL;
2554  rt_pgraster *pgrtn = NULL;
2555  rt_raster raster = NULL;
2556  rt_band band = NULL;
2557  int32_t bandindex;
2558 
2559  if (PG_ARGISNULL(0))
2560  PG_RETURN_NULL();
2561  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2562 
2563  raster = rt_raster_deserialize(pgraster, FALSE);
2564  if (!raster) {
2565  PG_FREE_IF_COPY(pgraster, 0);
2566  elog(ERROR, "RASTER_setBandIsNoData: Could not deserialize raster");
2567  PG_RETURN_NULL();
2568  }
2569 
2570  /* Check index is not NULL or smaller than 1 */
2571  if (PG_ARGISNULL(1))
2572  bandindex = -1;
2573  else
2574  bandindex = PG_GETARG_INT32(1);
2575 
2576  if (bandindex < 1)
2577  elog(NOTICE, "Invalid band index (must use 1-based). Isnodata flag not set. Returning original raster");
2578  else {
2579  /* Fetch requested band */
2580  band = rt_raster_get_band(raster, bandindex - 1);
2581 
2582  if (!band)
2583  elog(NOTICE, "Could not find raster band of index %d. Isnodata flag not set. Returning original raster", bandindex);
2584  else {
2585  if (!rt_band_get_hasnodata_flag(band)) {
2586  elog(NOTICE, "Band of index %d has no NODATA so cannot be NODATA. Returning original raster", bandindex);
2587  }
2588  /* Set the band's nodata value */
2589  else {
2590  rt_band_set_isnodata_flag(band, 1);
2591  }
2592  }
2593  }
2594 
2595  /* Serialize raster again */
2596  pgrtn = rt_raster_serialize(raster);
2597  rt_raster_destroy(raster);
2598  PG_FREE_IF_COPY(pgraster, 0);
2599  if (!pgrtn) PG_RETURN_NULL();
2600 
2601  SET_VARSIZE(pgrtn, pgrtn->size);
2602  PG_RETURN_POINTER(pgrtn);
2603 }
2604 
2606 Datum RASTER_bandIsNoData(PG_FUNCTION_ARGS)
2607 {
2608  rt_pgraster *pgraster = NULL;
2609  rt_raster raster = NULL;
2610  rt_band band = NULL;
2611  int32_t bandindex;
2612  bool forcechecking = FALSE;
2613  bool bandisnodata = FALSE;
2614 
2615  /* Index is 1-based */
2616  bandindex = PG_GETARG_INT32(1);
2617  if ( bandindex < 1 ) {
2618  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
2619  PG_RETURN_NULL();
2620  }
2621 
2622  /* Deserialize raster */
2623  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
2624  pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2625 
2626  raster = rt_raster_deserialize(pgraster, FALSE);
2627  if ( ! raster ) {
2628  PG_FREE_IF_COPY(pgraster, 0);
2629  elog(ERROR, "RASTER_bandIsNoData: Could not deserialize raster");
2630  PG_RETURN_NULL();
2631  }
2632 
2633  /* Fetch requested band and its nodata value */
2634  band = rt_raster_get_band(raster, bandindex - 1);
2635  if ( ! band ) {
2636  elog(NOTICE, "Could not find raster band of index %d when determining if band is nodata. Returning NULL", bandindex);
2637  rt_raster_destroy(raster);
2638  PG_FREE_IF_COPY(pgraster, 0);
2639  PG_RETURN_NULL();
2640  }
2641 
2642  forcechecking = PG_GETARG_BOOL(2);
2643 
2644  bandisnodata = (forcechecking) ?
2646 
2647  rt_raster_destroy(raster);
2648  PG_FREE_IF_COPY(pgraster, 0);
2649 
2650  PG_RETURN_BOOL(bandisnodata);
2651 }
2652 
2657 Datum RASTER_getBandPath(PG_FUNCTION_ARGS)
2658 {
2659  rt_pgraster *pgraster = NULL;
2660  rt_raster raster = NULL;
2661  rt_band band = NULL;
2662  int32_t bandindex;
2663  const char *bandpath;
2664  text *result;
2665 
2666  /* Index is 1-based */
2667  bandindex = PG_GETARG_INT32(1);
2668  if (bandindex < 1) {
2669  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
2670  PG_RETURN_NULL();
2671  }
2672 
2673  /* Deserialize raster */
2674  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
2675  pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2676 
2677  raster = rt_raster_deserialize(pgraster, FALSE);
2678  if (!raster) {
2679  PG_FREE_IF_COPY(pgraster, 0);
2680  elog(ERROR, "RASTER_getBandPath: Could not deserialize raster");
2681  PG_RETURN_NULL();
2682  }
2683 
2684  /* Fetch requested band and its nodata value */
2685  band = rt_raster_get_band(raster, bandindex - 1);
2686  if (!band) {
2687  elog(
2688  NOTICE,
2689  "Could not find raster band of index %d when getting band path. Returning NULL",
2690  bandindex
2691  );
2692  rt_raster_destroy(raster);
2693  PG_FREE_IF_COPY(pgraster, 0);
2694  PG_RETURN_NULL();
2695  }
2696 
2697  bandpath = rt_band_get_ext_path(band);
2698  if (!bandpath) {
2699  rt_band_destroy(band);
2700  rt_raster_destroy(raster);
2701  PG_FREE_IF_COPY(pgraster, 0);
2702  PG_RETURN_NULL();
2703  }
2704 
2705  result = (text *) palloc(VARHDRSZ + strlen(bandpath) + 1);
2706  SET_VARSIZE(result, VARHDRSZ + strlen(bandpath) + 1);
2707  strcpy((char *) VARDATA(result), bandpath);
2708 
2709  rt_band_destroy(band);
2710  rt_raster_destroy(raster);
2711  PG_FREE_IF_COPY(pgraster, 0);
2712 
2713  PG_RETURN_TEXT_P(result);
2714 }
2715 
2724 Datum RASTER_getPixelValue(PG_FUNCTION_ARGS)
2725 {
2726  rt_pgraster *pgraster = NULL;
2727  rt_raster raster = NULL;
2728  rt_band band = NULL;
2729  double pixvalue = 0;
2730  int32_t bandindex = 0;
2731  int32_t x = 0;
2732  int32_t y = 0;
2733  int result = 0;
2734  bool exclude_nodata_value = TRUE;
2735  int isnodata = 0;
2736 
2737  /* Index is 1-based */
2738  bandindex = PG_GETARG_INT32(1);
2739  if ( bandindex < 1 ) {
2740  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
2741  PG_RETURN_NULL();
2742  }
2743 
2744  x = PG_GETARG_INT32(2);
2745 
2746  y = PG_GETARG_INT32(3);
2747 
2748  exclude_nodata_value = PG_GETARG_BOOL(4);
2749 
2750  POSTGIS_RT_DEBUGF(3, "Pixel coordinates (%d, %d)", x, y);
2751 
2752  /* Deserialize raster */
2753  if (PG_ARGISNULL(0)) PG_RETURN_NULL();
2754  pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2755 
2756  raster = rt_raster_deserialize(pgraster, FALSE);
2757  if (!raster) {
2758  PG_FREE_IF_COPY(pgraster, 0);
2759  elog(ERROR, "RASTER_getPixelValue: Could not deserialize raster");
2760  PG_RETURN_NULL();
2761  }
2762 
2763  /* Fetch Nth band using 0-based internal index */
2764  band = rt_raster_get_band(raster, bandindex - 1);
2765  if (! band) {
2766  elog(NOTICE, "Could not find raster band of index %d when getting pixel "
2767  "value. Returning NULL", bandindex);
2768  rt_raster_destroy(raster);
2769  PG_FREE_IF_COPY(pgraster, 0);
2770  PG_RETURN_NULL();
2771  }
2772  /* Fetch pixel using 0-based coordinates */
2773  result = rt_band_get_pixel(band, x - 1, y - 1, &pixvalue, &isnodata);
2774 
2775  /* If the result is -1 or the value is nodata and we take nodata into account
2776  * then return nodata = NULL */
2777  if (result != ES_NONE || (exclude_nodata_value && isnodata)) {
2778  rt_raster_destroy(raster);
2779  PG_FREE_IF_COPY(pgraster, 0);
2780  PG_RETURN_NULL();
2781  }
2782 
2783  rt_raster_destroy(raster);
2784  PG_FREE_IF_COPY(pgraster, 0);
2785 
2786  PG_RETURN_FLOAT8(pixvalue);
2787 }
2788 
2789 /* ---------------------------------------------------------------- */
2790 /* ST_DumpValues function */
2791 /* ---------------------------------------------------------------- */
2792 
2796  int rows;
2797  int columns;
2798 
2799  int *nbands; /* 0-based */
2800  Datum **values;
2801  bool **nodata;
2802 };
2803 
2804 static rtpg_dumpvalues_arg rtpg_dumpvalues_arg_init() {
2805  rtpg_dumpvalues_arg arg = NULL;
2806 
2807  arg = palloc(sizeof(struct rtpg_dumpvalues_arg_t));
2808  if (arg == NULL) {
2809  elog(ERROR, "rtpg_dumpvalues_arg_init: Could not allocate memory for arguments");
2810  return NULL;
2811  }
2812 
2813  arg->numbands = 0;
2814  arg->rows = 0;
2815  arg->columns = 0;
2816 
2817  arg->nbands = NULL;
2818  arg->values = NULL;
2819  arg->nodata = NULL;
2820 
2821  return arg;
2822 }
2823 
2824 static void rtpg_dumpvalues_arg_destroy(rtpg_dumpvalues_arg arg) {
2825  int i = 0;
2826 
2827  if (arg->numbands) {
2828  if (arg->nbands != NULL)
2829  pfree(arg->nbands);
2830 
2831  if (arg->values != NULL) {
2832  for (i = 0; i < arg->numbands; i++) {
2833 
2834  if (arg->values[i] != NULL)
2835  pfree(arg->values[i]);
2836 
2837  if (arg->nodata[i] != NULL)
2838  pfree(arg->nodata[i]);
2839 
2840  }
2841 
2842  pfree(arg->values);
2843  }
2844 
2845  if (arg->nodata != NULL)
2846  pfree(arg->nodata);
2847  }
2848 
2849  pfree(arg);
2850 }
2851 
2853 Datum RASTER_dumpValues(PG_FUNCTION_ARGS)
2854 {
2855  FuncCallContext *funcctx;
2856  TupleDesc tupdesc;
2857  int call_cntr;
2858  int max_calls;
2859  int i = 0;
2860  int x = 0;
2861  int y = 0;
2862  int z = 0;
2863 
2864  int16 typlen;
2865  bool typbyval;
2866  char typalign;
2867 
2868  rtpg_dumpvalues_arg arg1 = NULL;
2869  rtpg_dumpvalues_arg arg2 = NULL;
2870 
2871  /* stuff done only on the first call of the function */
2872  if (SRF_IS_FIRSTCALL()) {
2873  MemoryContext oldcontext;
2874  rt_pgraster *pgraster = NULL;
2875  rt_raster raster = NULL;
2876  rt_band band = NULL;
2877  int numbands = 0;
2878  int j = 0;
2879  bool exclude_nodata_value = TRUE;
2880 
2881  ArrayType *array;
2882  Oid etype;
2883  Datum *e;
2884  bool *nulls;
2885 
2886  double val = 0;
2887  int isnodata = 0;
2888 
2889  POSTGIS_RT_DEBUG(2, "RASTER_dumpValues first call");
2890 
2891  /* create a function context for cross-call persistence */
2892  funcctx = SRF_FIRSTCALL_INIT();
2893 
2894  /* switch to memory context appropriate for multiple function calls */
2895  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
2896 
2897  /* Get input arguments */
2898  if (PG_ARGISNULL(0)) {
2899  MemoryContextSwitchTo(oldcontext);
2900  SRF_RETURN_DONE(funcctx);
2901  }
2902  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2903 
2904  raster = rt_raster_deserialize(pgraster, FALSE);
2905  if (!raster) {
2906  PG_FREE_IF_COPY(pgraster, 0);
2907  ereport(ERROR, (
2908  errcode(ERRCODE_OUT_OF_MEMORY),
2909  errmsg("Could not deserialize raster")
2910  ));
2911  MemoryContextSwitchTo(oldcontext);
2912  SRF_RETURN_DONE(funcctx);
2913  }
2914 
2915  /* check that raster is not empty */
2916  /*
2917  if (rt_raster_is_empty(raster)) {
2918  elog(NOTICE, "Raster provided is empty");
2919  rt_raster_destroy(raster);
2920  PG_FREE_IF_COPY(pgraster, 0);
2921  MemoryContextSwitchTo(oldcontext);
2922  SRF_RETURN_DONE(funcctx);
2923  }
2924  */
2925 
2926  /* raster has bands */
2927  numbands = rt_raster_get_num_bands(raster);
2928  if (!numbands) {
2929  elog(NOTICE, "Raster provided has no bands");
2930  rt_raster_destroy(raster);
2931  PG_FREE_IF_COPY(pgraster, 0);
2932  MemoryContextSwitchTo(oldcontext);
2933  SRF_RETURN_DONE(funcctx);
2934  }
2935 
2936  /* initialize arg1 */
2937  arg1 = rtpg_dumpvalues_arg_init();
2938  if (arg1 == NULL) {
2939  rt_raster_destroy(raster);
2940  PG_FREE_IF_COPY(pgraster, 0);
2941  MemoryContextSwitchTo(oldcontext);
2942  elog(ERROR, "RASTER_dumpValues: Could not initialize argument structure");
2943  SRF_RETURN_DONE(funcctx);
2944  }
2945 
2946  /* nband, array */
2947  if (!PG_ARGISNULL(1)) {
2948  array = PG_GETARG_ARRAYTYPE_P(1);
2949  etype = ARR_ELEMTYPE(array);
2950  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
2951 
2952  switch (etype) {
2953  case INT2OID:
2954  case INT4OID:
2955  break;
2956  default:
2958  rt_raster_destroy(raster);
2959  PG_FREE_IF_COPY(pgraster, 0);
2960  MemoryContextSwitchTo(oldcontext);
2961  elog(ERROR, "RASTER_dumpValues: Invalid data type for band indexes");
2962  SRF_RETURN_DONE(funcctx);
2963  break;
2964  }
2965 
2966  deconstruct_array(array, etype, typlen, typbyval, typalign, &e, &nulls, &(arg1->numbands));
2967 
2968  arg1->nbands = palloc(sizeof(int) * arg1->numbands);
2969  if (arg1->nbands == NULL) {
2971  rt_raster_destroy(raster);
2972  PG_FREE_IF_COPY(pgraster, 0);
2973  MemoryContextSwitchTo(oldcontext);
2974  elog(ERROR, "RASTER_dumpValues: Could not allocate memory for band indexes");
2975  SRF_RETURN_DONE(funcctx);
2976  }
2977 
2978  for (i = 0, j = 0; i < arg1->numbands; i++) {
2979  if (nulls[i]) continue;
2980 
2981  switch (etype) {
2982  case INT2OID:
2983  arg1->nbands[j] = DatumGetInt16(e[i]) - 1;
2984  break;
2985  case INT4OID:
2986  arg1->nbands[j] = DatumGetInt32(e[i]) - 1;
2987  break;
2988  }
2989 
2990  j++;
2991  }
2992 
2993  if (j < arg1->numbands) {
2994  arg1->nbands = repalloc(arg1->nbands, sizeof(int) * j);
2995  if (arg1->nbands == NULL) {
2997  rt_raster_destroy(raster);
2998  PG_FREE_IF_COPY(pgraster, 0);
2999  MemoryContextSwitchTo(oldcontext);
3000  elog(ERROR, "RASTER_dumpValues: Could not reallocate memory for band indexes");
3001  SRF_RETURN_DONE(funcctx);
3002  }
3003 
3004  arg1->numbands = j;
3005  }
3006 
3007  /* validate nbands */
3008  for (i = 0; i < arg1->numbands; i++) {
3009  if (!rt_raster_has_band(raster, arg1->nbands[i])) {
3010  elog(NOTICE, "Band at index %d not found in raster", arg1->nbands[i] + 1);
3012  rt_raster_destroy(raster);
3013  PG_FREE_IF_COPY(pgraster, 0);
3014  MemoryContextSwitchTo(oldcontext);
3015  SRF_RETURN_DONE(funcctx);
3016  }
3017  }
3018 
3019  }
3020  else {
3021  arg1->numbands = numbands;
3022  arg1->nbands = palloc(sizeof(int) * arg1->numbands);
3023 
3024  if (arg1->nbands == NULL) {
3026  rt_raster_destroy(raster);
3027  PG_FREE_IF_COPY(pgraster, 0);
3028  MemoryContextSwitchTo(oldcontext);
3029  elog(ERROR, "RASTER_dumpValues: Could not allocate memory for band indexes");
3030  SRF_RETURN_DONE(funcctx);
3031  }
3032 
3033  for (i = 0; i < arg1->numbands; i++) {
3034  arg1->nbands[i] = i;
3035  POSTGIS_RT_DEBUGF(4, "arg1->nbands[%d] = %d", arg1->nbands[i], i);
3036  }
3037  }
3038 
3039  arg1->rows = rt_raster_get_height(raster);
3040  arg1->columns = rt_raster_get_width(raster);
3041 
3042  /* exclude_nodata_value */
3043  if (!PG_ARGISNULL(2))
3044  exclude_nodata_value = PG_GETARG_BOOL(2);
3045  POSTGIS_RT_DEBUGF(4, "exclude_nodata_value = %d", exclude_nodata_value);
3046 
3047  /* allocate memory for each band's values and nodata flags */
3048  arg1->values = palloc(sizeof(Datum *) * arg1->numbands);
3049  arg1->nodata = palloc(sizeof(bool *) * arg1->numbands);
3050  if (arg1->values == NULL || arg1->nodata == NULL) {
3052  rt_raster_destroy(raster);
3053  PG_FREE_IF_COPY(pgraster, 0);
3054  MemoryContextSwitchTo(oldcontext);
3055  elog(ERROR, "RASTER_dumpValues: Could not allocate memory for pixel values");
3056  SRF_RETURN_DONE(funcctx);
3057  }
3058  memset(arg1->values, 0, sizeof(Datum *) * arg1->numbands);
3059  memset(arg1->nodata, 0, sizeof(bool *) * arg1->numbands);
3060 
3061  /* get each band and dump data */
3062  for (z = 0; z < arg1->numbands; z++) {
3063  /* shortcut if raster is empty */
3064  if (rt_raster_is_empty(raster))
3065  break;
3066 
3067  band = rt_raster_get_band(raster, arg1->nbands[z]);
3068  if (!band) {
3069  int nband = arg1->nbands[z] + 1;
3071  rt_raster_destroy(raster);
3072  PG_FREE_IF_COPY(pgraster, 0);
3073  MemoryContextSwitchTo(oldcontext);
3074  elog(ERROR, "RASTER_dumpValues: Could not get band at index %d", nband);
3075  SRF_RETURN_DONE(funcctx);
3076  }
3077 
3078  /* allocate memory for values and nodata flags */
3079  arg1->values[z] = palloc(sizeof(Datum) * arg1->rows * arg1->columns);
3080  arg1->nodata[z] = palloc(sizeof(bool) * arg1->rows * arg1->columns);
3081  if (arg1->values[z] == NULL || arg1->nodata[z] == NULL) {
3083  rt_raster_destroy(raster);
3084  PG_FREE_IF_COPY(pgraster, 0);
3085  MemoryContextSwitchTo(oldcontext);
3086  elog(ERROR, "RASTER_dumpValues: Could not allocate memory for pixel values");
3087  SRF_RETURN_DONE(funcctx);
3088  }
3089  memset(arg1->values[z], 0, sizeof(Datum) * arg1->rows * arg1->columns);
3090  memset(arg1->nodata[z], 0, sizeof(bool) * arg1->rows * arg1->columns);
3091 
3092  i = 0;
3093 
3094  /* shortcut if band is NODATA */
3095  if (rt_band_get_isnodata_flag(band)) {
3096  for (i = (arg1->rows * arg1->columns) - 1; i >= 0; i--)
3097  arg1->nodata[z][i] = TRUE;
3098  continue;
3099  }
3100 
3101  for (y = 0; y < arg1->rows; y++) {
3102  for (x = 0; x < arg1->columns; x++) {
3103  /* get pixel */
3104  if (rt_band_get_pixel(band, x, y, &val, &isnodata) != ES_NONE) {
3105  int nband = arg1->nbands[z] + 1;
3107  rt_raster_destroy(raster);
3108  PG_FREE_IF_COPY(pgraster, 0);
3109  MemoryContextSwitchTo(oldcontext);
3110  elog(ERROR, "RASTER_dumpValues: Could not pixel (%d, %d) of band %d", x, y, nband);
3111  SRF_RETURN_DONE(funcctx);
3112  }
3113 
3114  arg1->values[z][i] = Float8GetDatum(val);
3115  POSTGIS_RT_DEBUGF(5, "arg1->values[z][i] = %f", DatumGetFloat8(arg1->values[z][i]));
3116  POSTGIS_RT_DEBUGF(5, "clamped is?: %d", rt_band_clamped_value_is_nodata(band, val));
3117 
3118  if (exclude_nodata_value && isnodata) {
3119  arg1->nodata[z][i] = TRUE;
3120  POSTGIS_RT_DEBUG(5, "nodata = 1");
3121  }
3122  else
3123  POSTGIS_RT_DEBUG(5, "nodata = 0");
3124 
3125  i++;
3126  }
3127  }
3128  }
3129 
3130  /* cleanup */
3131  rt_raster_destroy(raster);
3132  PG_FREE_IF_COPY(pgraster, 0);
3133 
3134  /* Store needed information */
3135  funcctx->user_fctx = arg1;
3136 
3137  /* total number of tuples to be returned */
3138  funcctx->max_calls = arg1->numbands;
3139 
3140  /* Build a tuple descriptor for our result type */
3141  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
3142  MemoryContextSwitchTo(oldcontext);
3143  ereport(ERROR, (
3144  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
3145  errmsg(
3146  "function returning record called in context "
3147  "that cannot accept type record"
3148  )
3149  ));
3150  }
3151 
3152  BlessTupleDesc(tupdesc);
3153  funcctx->tuple_desc = tupdesc;
3154 
3155  MemoryContextSwitchTo(oldcontext);
3156  }
3157 
3158  /* stuff done on every call of the function */
3159  funcctx = SRF_PERCALL_SETUP();
3160 
3161  call_cntr = funcctx->call_cntr;
3162  max_calls = funcctx->max_calls;
3163  tupdesc = funcctx->tuple_desc;
3164  arg2 = funcctx->user_fctx;
3165 
3166  /* do when there is more left to send */
3167  if (call_cntr < max_calls) {
3168  int values_length = 2;
3169  Datum values[values_length];
3170  bool nulls[values_length];
3171  HeapTuple tuple;
3172  Datum result;
3173  ArrayType *mdValues = NULL;
3174  int ndim = 2;
3175  int dim[2] = {arg2->rows, arg2->columns};
3176  int lbound[2] = {1, 1};
3177 
3178  POSTGIS_RT_DEBUGF(3, "call number %d", call_cntr);
3179  POSTGIS_RT_DEBUGF(4, "dim = %d, %d", dim[0], dim[1]);
3180 
3181  memset(nulls, FALSE, sizeof(bool) * values_length);
3182 
3183  values[0] = Int32GetDatum(arg2->nbands[call_cntr] + 1);
3184 
3185  /* info about the type of item in the multi-dimensional array (float8). */
3186  get_typlenbyvalalign(FLOAT8OID, &typlen, &typbyval, &typalign);
3187 
3188  /* if values is NULL, return empty raster */
3189  if (arg2->values[call_cntr] == NULL)
3190  ndim = 0;
3191 
3192  /* assemble 3-dimension array of values */
3193  mdValues = construct_md_array(
3194  arg2->values[call_cntr], arg2->nodata[call_cntr],
3195  ndim, dim, lbound,
3196  FLOAT8OID,
3197  typlen, typbyval, typalign
3198  );
3199  values[1] = PointerGetDatum(mdValues);
3200 
3201  /* build a tuple and datum */
3202  tuple = heap_form_tuple(tupdesc, values, nulls);
3203  result = HeapTupleGetDatum(tuple);
3204 
3205  SRF_RETURN_NEXT(funcctx, result);
3206  }
3207  /* do when there is no more left */
3208  else {
3210  SRF_RETURN_DONE(funcctx);
3211  }
3212 }
3213 
3218 Datum RASTER_setPixelValue(PG_FUNCTION_ARGS)
3219 {
3220  rt_pgraster *pgraster = NULL;
3221  rt_pgraster *pgrtn = NULL;
3222  rt_raster raster = NULL;
3223  rt_band band = NULL;
3224  double pixvalue = 0;
3225  int32_t bandindex = 0;
3226  int32_t x = 0;
3227  int32_t y = 0;
3228  bool skipset = FALSE;
3229 
3230  if (PG_ARGISNULL(0))
3231  PG_RETURN_NULL();
3232 
3233  /* Check index is not NULL or < 1 */
3234  if (PG_ARGISNULL(1))
3235  bandindex = -1;
3236  else
3237  bandindex = PG_GETARG_INT32(1);
3238 
3239  if (bandindex < 1) {
3240  elog(NOTICE, "Invalid band index (must use 1-based). Value not set. Returning original raster");
3241  skipset = TRUE;
3242  }
3243 
3244  /* Validate pixel coordinates are not null */
3245  if (PG_ARGISNULL(2)) {
3246  elog(NOTICE, "X coordinate can not be NULL when setting pixel value. Value not set. Returning original raster");
3247  skipset = TRUE;
3248  }
3249  else
3250  x = PG_GETARG_INT32(2);
3251 
3252  if (PG_ARGISNULL(3)) {
3253  elog(NOTICE, "Y coordinate can not be NULL when setting pixel value. Value not set. Returning original raster");
3254  skipset = TRUE;
3255  }
3256  else
3257  y = PG_GETARG_INT32(3);
3258 
3259  POSTGIS_RT_DEBUGF(3, "Pixel coordinates (%d, %d)", x, y);
3260 
3261  /* Deserialize raster */
3262  pgraster = (rt_pgraster *) PG_DETOAST_DATUM_COPY(PG_GETARG_DATUM(0));
3263 
3264  raster = rt_raster_deserialize(pgraster, FALSE);
3265  if (!raster) {
3266  PG_FREE_IF_COPY(pgraster, 0);
3267  elog(ERROR, "RASTER_setPixelValue: Could not deserialize raster");
3268  PG_RETURN_NULL();
3269  }
3270 
3271  if (!skipset) {
3272  /* Fetch requested band */
3273  band = rt_raster_get_band(raster, bandindex - 1);
3274  if (!band) {
3275  elog(NOTICE, "Could not find raster band of index %d when setting "
3276  "pixel value. Value not set. Returning original raster",
3277  bandindex);
3278  PG_RETURN_POINTER(pgraster);
3279  }
3280  else {
3281  /* Set the pixel value */
3282  if (PG_ARGISNULL(4)) {
3283  if (!rt_band_get_hasnodata_flag(band)) {
3284  elog(NOTICE, "Raster do not have a nodata value defined. "
3285  "Set band nodata value first. Nodata value not set. "
3286  "Returning original raster");
3287  PG_RETURN_POINTER(pgraster);
3288  }
3289  else {
3290  rt_band_get_nodata(band, &pixvalue);
3291  rt_band_set_pixel(band, x - 1, y - 1, pixvalue, NULL);
3292  }
3293  }
3294  else {
3295  pixvalue = PG_GETARG_FLOAT8(4);
3296  rt_band_set_pixel(band, x - 1, y - 1, pixvalue, NULL);
3297  }
3298  }
3299  }
3300 
3301  pgrtn = rt_raster_serialize(raster);
3302  rt_raster_destroy(raster);
3303  PG_FREE_IF_COPY(pgraster, 0);
3304  if (!pgrtn)
3305  PG_RETURN_NULL();
3306 
3307  SET_VARSIZE(pgrtn, pgrtn->size);
3308  PG_RETURN_POINTER(pgrtn);
3309 }
3310 
3315 Datum RASTER_setPixelValuesArray(PG_FUNCTION_ARGS)
3316 {
3317  rt_pgraster *pgraster = NULL;
3318  rt_pgraster *pgrtn = NULL;
3319  rt_raster raster = NULL;
3320  rt_band band = NULL;
3321  int numbands = 0;
3322 
3323  int nband = 0;
3324  int width = 0;
3325  int height = 0;
3326 
3327  ArrayType *array;
3328  Oid etype;
3329  Datum *elements;
3330  bool *nulls;
3331  int16 typlen;
3332  bool typbyval;
3333  char typalign;
3334  int ndims = 1;
3335  int *dims;
3336  int num = 0;
3337 
3338  int ul[2] = {0};
3339  struct pixelvalue {
3340  int x;
3341  int y;
3342 
3343  bool noset;
3344  bool nodata;
3345  double value;
3346  };
3347  struct pixelvalue *pixval = NULL;
3348  int numpixval = 0;
3349  int dimpixval[2] = {1, 1};
3350  int dimnoset[2] = {1, 1};
3351  int hasnodata = FALSE;
3352  double nodataval = 0;
3353  bool keepnodata = FALSE;
3354  bool hasnosetval = FALSE;
3355  bool nosetvalisnull = FALSE;
3356  double nosetval = 0;
3357 
3358  int rtn = 0;
3359  double val = 0;
3360  int isnodata = 0;
3361 
3362  int i = 0;
3363  int j = 0;
3364  int x = 0;
3365  int y = 0;
3366 
3367  /* pgraster is null, return null */
3368  if (PG_ARGISNULL(0))
3369  PG_RETURN_NULL();
3370  pgraster = (rt_pgraster *) PG_DETOAST_DATUM_COPY(PG_GETARG_DATUM(0));
3371 
3372  /* raster */
3373  raster = rt_raster_deserialize(pgraster, FALSE);
3374  if (!raster) {
3375  PG_FREE_IF_COPY(pgraster, 0);
3376  elog(ERROR, "RASTER_setPixelValuesArray: Could not deserialize raster");
3377  PG_RETURN_NULL();
3378  }
3379 
3380  /* raster attributes */
3381  numbands = rt_raster_get_num_bands(raster);
3382  width = rt_raster_get_width(raster);
3383  height = rt_raster_get_height(raster);
3384 
3385  /* nband */
3386  if (PG_ARGISNULL(1)) {
3387  elog(NOTICE, "Band index cannot be NULL. Value must be 1-based. Returning original raster");
3388  rt_raster_destroy(raster);
3389  PG_RETURN_POINTER(pgraster);
3390  }
3391 
3392  nband = PG_GETARG_INT32(1);
3393  if (nband < 1 || nband > numbands) {
3394  elog(NOTICE, "Band index is invalid. Value must be 1-based. Returning original raster");
3395  rt_raster_destroy(raster);
3396  PG_RETURN_POINTER(pgraster);
3397  }
3398 
3399  /* x, y */
3400  for (i = 2, j = 0; i < 4; i++, j++) {
3401  if (PG_ARGISNULL(i)) {
3402  elog(NOTICE, "%s cannot be NULL. Value must be 1-based. Returning original raster", j < 1 ? "X" : "Y");
3403  rt_raster_destroy(raster);
3404  PG_RETURN_POINTER(pgraster);
3405  }
3406 
3407  ul[j] = PG_GETARG_INT32(i);
3408  if (
3409  (ul[j] < 1) || (
3410  (j < 1 && ul[j] > width) ||
3411  (j > 0 && ul[j] > height)
3412  )
3413  ) {
3414  elog(NOTICE, "%s is invalid. Value must be 1-based. Returning original raster", j < 1 ? "X" : "Y");
3415  rt_raster_destroy(raster);
3416  PG_RETURN_POINTER(pgraster);
3417  }
3418 
3419  /* force 0-based from 1-based */
3420  ul[j] -= 1;
3421  }
3422 
3423  /* new value set */
3424  if (PG_ARGISNULL(4)) {
3425  elog(NOTICE, "No values to set. Returning original raster");
3426  rt_raster_destroy(raster);
3427  PG_RETURN_POINTER(pgraster);
3428  }
3429 
3430  array = PG_GETARG_ARRAYTYPE_P(4);
3431  etype = ARR_ELEMTYPE(array);
3432  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3433 
3434  switch (etype) {
3435  case FLOAT4OID:
3436  case FLOAT8OID:
3437  break;
3438  default:
3439  rt_raster_destroy(raster);
3440  PG_FREE_IF_COPY(pgraster, 0);
3441  elog(ERROR, "RASTER_setPixelValuesArray: Invalid data type for new values");
3442  PG_RETURN_NULL();
3443  break;
3444  }
3445 
3446  ndims = ARR_NDIM(array);
3447  dims = ARR_DIMS(array);
3448  POSTGIS_RT_DEBUGF(4, "ndims = %d", ndims);
3449 
3450  if (ndims < 1 || ndims > 2) {
3451  elog(NOTICE, "New values array must be of 1 or 2 dimensions. Returning original raster");
3452  rt_raster_destroy(raster);
3453  PG_RETURN_POINTER(pgraster);
3454  }
3455  /* outer element, then inner element */
3456  /* i = 0, y */
3457  /* i = 1, x */
3458  if (ndims != 2)
3459  dimpixval[1] = dims[0];
3460  else {
3461  dimpixval[0] = dims[0];
3462  dimpixval[1] = dims[1];
3463  }
3464  POSTGIS_RT_DEBUGF(4, "dimpixval = (%d, %d)", dimpixval[0], dimpixval[1]);
3465 
3466  deconstruct_array(
3467  array,
3468  etype,
3469  typlen, typbyval, typalign,
3470  &elements, &nulls, &num
3471  );
3472 
3473  /* # of elements doesn't match dims */
3474  if (num < 1 || num != (dimpixval[0] * dimpixval[1])) {
3475  if (num) {
3476  pfree(elements);
3477  pfree(nulls);
3478  }
3479  rt_raster_destroy(raster);
3480  PG_FREE_IF_COPY(pgraster, 0);
3481  elog(ERROR, "RASTER_setPixelValuesArray: Could not deconstruct new values array");
3482  PG_RETURN_NULL();
3483  }
3484 
3485  /* allocate memory for pixval */
3486  numpixval = num;
3487  pixval = palloc(sizeof(struct pixelvalue) * numpixval);
3488  if (pixval == NULL) {
3489  pfree(elements);
3490  pfree(nulls);
3491  rt_raster_destroy(raster);
3492  PG_FREE_IF_COPY(pgraster, 0);
3493  elog(ERROR, "RASTER_setPixelValuesArray: Could not allocate memory for new pixel values");
3494  PG_RETURN_NULL();
3495  }
3496 
3497  /* load new values into pixval */
3498  i = 0;
3499  for (y = 0; y < dimpixval[0]; y++) {
3500  for (x = 0; x < dimpixval[1]; x++) {
3501  /* 0-based */
3502  pixval[i].x = ul[0] + x;
3503  pixval[i].y = ul[1] + y;
3504 
3505  pixval[i].noset = FALSE;
3506  pixval[i].nodata = FALSE;
3507  pixval[i].value = 0;
3508 
3509  if (nulls[i])
3510  pixval[i].nodata = TRUE;
3511  else {
3512  switch (etype) {
3513  case FLOAT4OID:
3514  pixval[i].value = DatumGetFloat4(elements[i]);
3515  break;
3516  case FLOAT8OID:
3517  pixval[i].value = DatumGetFloat8(elements[i]);
3518  break;
3519  }
3520  }
3521 
3522  i++;
3523  }
3524  }
3525 
3526  pfree(elements);
3527  pfree(nulls);
3528 
3529  /* now load noset flags */
3530  if (!PG_ARGISNULL(5)) {
3531  array = PG_GETARG_ARRAYTYPE_P(5);
3532  etype = ARR_ELEMTYPE(array);
3533  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3534 
3535  switch (etype) {
3536  case BOOLOID:
3537  break;
3538  default:
3539  pfree(pixval);
3540  rt_raster_destroy(raster);
3541  PG_FREE_IF_COPY(pgraster, 0);
3542  elog(ERROR, "RASTER_setPixelValuesArray: Invalid data type for noset flags");
3543  PG_RETURN_NULL();
3544  break;
3545  }
3546 
3547  ndims = ARR_NDIM(array);
3548  dims = ARR_DIMS(array);
3549  POSTGIS_RT_DEBUGF(4, "ndims = %d", ndims);
3550 
3551  if (ndims < 1 || ndims > 2) {
3552  elog(NOTICE, "Noset flags array must be of 1 or 2 dimensions. Returning original raster");
3553  pfree(pixval);
3554  rt_raster_destroy(raster);
3555  PG_RETURN_POINTER(pgraster);
3556  }
3557  /* outer element, then inner element */
3558  /* i = 0, y */
3559  /* i = 1, x */
3560  if (ndims != 2)
3561  dimnoset[1] = dims[0];
3562  else {
3563  dimnoset[0] = dims[0];
3564  dimnoset[1] = dims[1];
3565  }
3566  POSTGIS_RT_DEBUGF(4, "dimnoset = (%d, %d)", dimnoset[0], dimnoset[1]);
3567 
3568  deconstruct_array(
3569  array,
3570  etype,
3571  typlen, typbyval, typalign,
3572  &elements, &nulls, &num
3573  );
3574 
3575  /* # of elements doesn't match dims */
3576  if (num < 1 || num != (dimnoset[0] * dimnoset[1])) {
3577  pfree(pixval);
3578  if (num) {
3579  pfree(elements);
3580  pfree(nulls);
3581  }
3582  rt_raster_destroy(raster);
3583  PG_FREE_IF_COPY(pgraster, 0);
3584  elog(ERROR, "RASTER_setPixelValuesArray: Could not deconstruct noset flags array");
3585  PG_RETURN_NULL();
3586  }
3587 
3588  i = 0;
3589  j = 0;
3590  for (y = 0; y < dimnoset[0]; y++) {
3591  if (y >= dimpixval[0]) break;
3592 
3593  for (x = 0; x < dimnoset[1]; x++) {
3594  /* fast forward noset elements */
3595  if (x >= dimpixval[1]) {
3596  i += (dimnoset[1] - dimpixval[1]);
3597  break;
3598  }
3599 
3600  if (!nulls[i] && DatumGetBool(elements[i]))
3601  pixval[j].noset = TRUE;
3602 
3603  i++;
3604  j++;
3605  }
3606 
3607  /* fast forward pixval */
3608  if (x < dimpixval[1])
3609  j += (dimpixval[1] - dimnoset[1]);
3610  }
3611 
3612  pfree(elements);
3613  pfree(nulls);
3614  }
3615  /* hasnosetvalue and nosetvalue */
3616  else if (!PG_ARGISNULL(6) && PG_GETARG_BOOL(6)) {
3617  hasnosetval = TRUE;
3618  if (PG_ARGISNULL(7))
3619  nosetvalisnull = TRUE;
3620  else
3621  nosetval = PG_GETARG_FLOAT8(7);
3622  }
3623 
3624 #if POSTGIS_DEBUG_LEVEL > 0
3625  for (i = 0; i < numpixval; i++) {
3626  POSTGIS_RT_DEBUGF(4, "pixval[%d](x, y, noset, nodata, value) = (%d, %d, %d, %d, %f)",
3627  i,
3628  pixval[i].x,
3629  pixval[i].y,
3630  pixval[i].noset,
3631  pixval[i].nodata,
3632  pixval[i].value
3633  );
3634  }
3635 #endif
3636 
3637  /* keepnodata flag */
3638  if (!PG_ARGISNULL(8))
3639  keepnodata = PG_GETARG_BOOL(8);
3640 
3641  /* get band */
3642  band = rt_raster_get_band(raster, nband - 1);
3643  if (!band) {
3644  elog(NOTICE, "Could not find band at index %d. Returning original raster", nband);
3645  pfree(pixval);
3646  rt_raster_destroy(raster);
3647  PG_RETURN_POINTER(pgraster);
3648  }
3649 
3650  /* get band nodata info */
3651  /* has NODATA, use NODATA */
3652  hasnodata = rt_band_get_hasnodata_flag(band);
3653  if (hasnodata)
3654  rt_band_get_nodata(band, &nodataval);
3655  /* no NODATA, use min possible value */
3656  else
3657  nodataval = rt_band_get_min_value(band);
3658 
3659  /* set pixels */
3660  for (i = 0; i < numpixval; i++) {
3661  /* noset = true, skip */
3662  if (pixval[i].noset)
3663  continue;
3664  /* check against nosetval */
3665  else if (hasnosetval) {
3666  /* pixel = NULL AND nosetval = NULL */
3667  if (pixval[i].nodata && nosetvalisnull)
3668  continue;
3669  /* pixel value = nosetval */
3670  else if (!pixval[i].nodata && !nosetvalisnull && FLT_EQ(pixval[i].value, nosetval))
3671  continue;
3672  }
3673 
3674  /* if pixel is outside bounds, skip */
3675  if (
3676  (pixval[i].x < 0 || pixval[i].x >= width) ||
3677  (pixval[i].y < 0 || pixval[i].y >= height)
3678  ) {
3679  elog(NOTICE, "Cannot set value for pixel (%d, %d) outside raster bounds: %d x %d",
3680  pixval[i].x + 1, pixval[i].y + 1,
3681  width, height
3682  );
3683  continue;
3684  }
3685 
3686  /* if hasnodata = TRUE and keepnodata = TRUE, inspect pixel value */
3687  if (hasnodata && keepnodata) {
3688  rtn = rt_band_get_pixel(band, pixval[i].x, pixval[i].y, &val, &isnodata);
3689  if (rtn != ES_NONE) {
3690  pfree(pixval);
3691  rt_raster_destroy(raster);
3692  PG_FREE_IF_COPY(pgraster, 0);
3693  elog(ERROR, "Cannot get value of pixel");
3694  PG_RETURN_NULL();
3695  }
3696 
3697  /* pixel value = NODATA, skip */
3698  if (isnodata) {
3699  continue;
3700  }
3701  }
3702 
3703  if (pixval[i].nodata)
3704  rt_band_set_pixel(band, pixval[i].x, pixval[i].y, nodataval, NULL);
3705  else
3706  rt_band_set_pixel(band, pixval[i].x, pixval[i].y, pixval[i].value, NULL);
3707  }
3708 
3709  pfree(pixval);
3710 
3711  /* serialize new raster */
3712  pgrtn = rt_raster_serialize(raster);
3713  rt_raster_destroy(raster);
3714  PG_FREE_IF_COPY(pgraster, 0);
3715  if (!pgrtn)
3716  PG_RETURN_NULL();
3717 
3718  SET_VARSIZE(pgrtn, pgrtn->size);
3719  PG_RETURN_POINTER(pgrtn);
3720 }
3721 
3722 /* ---------------------------------------------------------------- */
3723 /* ST_SetValues using geomval array */
3724 /* ---------------------------------------------------------------- */
3725 
3728 
3730  int ngv;
3731  rtpg_setvaluesgv_geomval gv;
3732 
3734 };
3735 
3737  struct {
3738  int nodata;
3739  double value;
3740  } pixval;
3741 
3744 };
3745 
3746 static rtpg_setvaluesgv_arg rtpg_setvaluesgv_arg_init() {
3747  rtpg_setvaluesgv_arg arg = palloc(sizeof(struct rtpg_setvaluesgv_arg_t));
3748  if (arg == NULL) {
3749  elog(ERROR, "rtpg_setvaluesgv_arg_init: Could not allocate memory for function arguments");
3750  return NULL;
3751  }
3752 
3753  arg->ngv = 0;
3754  arg->gv = NULL;
3755  arg->keepnodata = 0;
3756 
3757  return arg;
3758 }
3759 
3760 static void rtpg_setvaluesgv_arg_destroy(rtpg_setvaluesgv_arg arg) {
3761  int i = 0;
3762 
3763  if (arg->gv != NULL) {
3764  for (i = 0; i < arg->ngv; i++) {
3765  if (arg->gv[i].geom != NULL)
3766  lwgeom_free(arg->gv[i].geom);
3767  if (arg->gv[i].mask != NULL)
3768  rt_raster_destroy(arg->gv[i].mask);
3769  }
3770 
3771  pfree(arg->gv);
3772  }
3773 
3774  pfree(arg);
3775 }
3776 
3778  rt_iterator_arg arg, void *userarg,
3779  double *value, int *nodata
3780 ) {
3781  rtpg_setvaluesgv_arg funcarg = (rtpg_setvaluesgv_arg) userarg;
3782  int i = 0;
3783  int j = 0;
3784 
3785  *value = 0;
3786  *nodata = 0;
3787 
3788  POSTGIS_RT_DEBUGF(4, "keepnodata = %d", funcarg->keepnodata);
3789 
3790  /* keepnodata = TRUE */
3791  if (funcarg->keepnodata && arg->nodata[0][0][0]) {
3792  POSTGIS_RT_DEBUG(3, "keepnodata = 1 AND source raster pixel is NODATA");
3793  *nodata = 1;
3794  return 1;
3795  }
3796 
3797  for (i = arg->rasters - 1, j = funcarg->ngv - 1; i > 0; i--, j--) {
3798  POSTGIS_RT_DEBUGF(4, "checking raster %d", i);
3799 
3800  /* mask is NODATA */
3801  if (arg->nodata[i][0][0])
3802  continue;
3803  /* mask is NOT NODATA */
3804  else {
3805  POSTGIS_RT_DEBUGF(4, "Using information from geometry %d", j);
3806 
3807  if (funcarg->gv[j].pixval.nodata)
3808  *nodata = 1;
3809  else
3810  *value = funcarg->gv[j].pixval.value;
3811 
3812  return 1;
3813  }
3814  }
3815 
3816  POSTGIS_RT_DEBUG(4, "Using information from source raster");
3817 
3818  /* still here */
3819  if (arg->nodata[0][0][0])
3820  *nodata = 1;
3821  else
3822  *value = arg->values[0][0][0];
3823 
3824  return 1;
3825 }
3826 
3828 Datum RASTER_setPixelValuesGeomval(PG_FUNCTION_ARGS)
3829 {
3830  rt_pgraster *pgraster = NULL;
3831  rt_pgraster *pgrtn = NULL;
3832  rt_raster raster = NULL;
3833  rt_band band = NULL;
3834  rt_raster _raster = NULL;
3835  rt_band _band = NULL;
3836  int nband = 0; /* 1-based */
3837 
3838  int numbands = 0;
3839  int width = 0;
3840  int height = 0;
3841  int srid = 0;
3842  double gt[6] = {0};
3843 
3844  rt_pixtype pixtype = PT_END;
3845  int hasnodata = 0;
3846  double nodataval = 0;
3847 
3848  rtpg_setvaluesgv_arg arg = NULL;
3849  int allpoint = 0;
3850 
3851  ArrayType *array;
3852  Oid etype;
3853  Datum *e;
3854  bool *nulls;
3855  int16 typlen;
3856  bool typbyval;
3857  char typalign;
3858  int n = 0;
3859 
3860  HeapTupleHeader tup;
3861  bool isnull;
3862  Datum tupv;
3863 
3864  GSERIALIZED *gser = NULL;
3865  uint8_t gtype;
3866  unsigned char *wkb = NULL;
3867  size_t wkb_len;
3868 
3869  int i = 0;
3870  int j = 0;
3871  int noerr = 1;
3872 
3873  /* pgraster is null, return null */
3874  if (PG_ARGISNULL(0))
3875  PG_RETURN_NULL();
3876  pgraster = (rt_pgraster *) PG_DETOAST_DATUM_COPY(PG_GETARG_DATUM(0));
3877 
3878  /* raster */
3879  raster = rt_raster_deserialize(pgraster, FALSE);
3880  if (!raster) {
3881  PG_FREE_IF_COPY(pgraster, 0);
3882  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not deserialize raster");
3883  PG_RETURN_NULL();
3884  }
3885 
3886  /* raster attributes */
3887  numbands = rt_raster_get_num_bands(raster);
3888  width = rt_raster_get_width(raster);
3889  height = rt_raster_get_height(raster);
3890  srid = clamp_srid(rt_raster_get_srid(raster));
3892 
3893  /* nband */
3894  if (PG_ARGISNULL(1)) {
3895  elog(NOTICE, "Band index cannot be NULL. Value must be 1-based. Returning original raster");
3896  rt_raster_destroy(raster);
3897  PG_RETURN_POINTER(pgraster);
3898  }
3899 
3900  nband = PG_GETARG_INT32(1);
3901  if (nband < 1 || nband > numbands) {
3902  elog(NOTICE, "Band index is invalid. Value must be 1-based. Returning original raster");
3903  rt_raster_destroy(raster);
3904  PG_RETURN_POINTER(pgraster);
3905  }
3906 
3907  /* get band attributes */
3908  band = rt_raster_get_band(raster, nband - 1);
3909  pixtype = rt_band_get_pixtype(band);
3910  hasnodata = rt_band_get_hasnodata_flag(band);
3911  if (hasnodata)
3912  rt_band_get_nodata(band, &nodataval);
3913 
3914  /* array of geomval (2) */
3915  if (PG_ARGISNULL(2)) {
3916  elog(NOTICE, "No values to set. Returning original raster");
3917  rt_raster_destroy(raster);
3918  PG_RETURN_POINTER(pgraster);
3919  }
3920 
3921  array = PG_GETARG_ARRAYTYPE_P(2);
3922  etype = ARR_ELEMTYPE(array);
3923  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3924 
3925  deconstruct_array(
3926  array,
3927  etype,
3928  typlen, typbyval, typalign,
3929  &e, &nulls, &n
3930  );
3931 
3932  if (!n) {
3933  elog(NOTICE, "No values to set. Returning original raster");
3934  rt_raster_destroy(raster);
3935  PG_RETURN_POINTER(pgraster);
3936  }
3937 
3938  /* init arg */
3939  arg = rtpg_setvaluesgv_arg_init();
3940  if (arg == NULL) {
3941  rt_raster_destroy(raster);
3942  PG_FREE_IF_COPY(pgraster, 0);
3943  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not intialize argument structure");
3944  PG_RETURN_NULL();
3945  }
3946 
3947  arg->gv = palloc(sizeof(struct rtpg_setvaluesgv_geomval_t) * n);
3948  if (arg->gv == NULL) {
3950  rt_raster_destroy(raster);
3951  PG_FREE_IF_COPY(pgraster, 0);
3952  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not allocate memory for geomval array");
3953  PG_RETURN_NULL();
3954  }
3955 
3956  /* process each element */
3957  arg->ngv = 0;
3958  for (i = 0; i < n; i++) {
3959  if (nulls[i])
3960  continue;
3961 
3962  arg->gv[arg->ngv].pixval.nodata = 0;
3963  arg->gv[arg->ngv].pixval.value = 0;
3964  arg->gv[arg->ngv].geom = NULL;
3965  arg->gv[arg->ngv].mask = NULL;
3966 
3967  /* each element is a tuple */
3968  tup = (HeapTupleHeader) DatumGetPointer(e[i]);
3969  if (NULL == tup) {
3971  rt_raster_destroy(raster);
3972  PG_FREE_IF_COPY(pgraster, 0);
3973  elog(ERROR, "RASTER_setPixelValuesGeomval: Invalid argument for geomval at index %d", i);
3974  PG_RETURN_NULL();
3975  }
3976 
3977  /* first element, geometry */
3978  POSTGIS_RT_DEBUG(4, "Processing first element (geometry)");
3979  tupv = GetAttributeByName(tup, "geom", &isnull);
3980  if (isnull) {
3981  elog(NOTICE, "First argument (geom) of geomval at index %d is NULL. Skipping", i);
3982  continue;
3983  }
3984 
3985  gser = (GSERIALIZED *) PG_DETOAST_DATUM(tupv);
3986  arg->gv[arg->ngv].geom = lwgeom_from_gserialized(gser);
3987  if (arg->gv[arg->ngv].geom == NULL) {
3989  rt_raster_destroy(raster);
3990  PG_FREE_IF_COPY(pgraster, 0);
3991  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not deserialize geometry of geomval at index %d", i);
3992  PG_RETURN_NULL();
3993  }
3994 
3995  /* empty geometry */
3996  if (lwgeom_is_empty(arg->gv[arg->ngv].geom)) {
3997  elog(NOTICE, "First argument (geom) of geomval at index %d is an empty geometry. Skipping", i);
3998  continue;
3999  }
4000 
4001  /* check SRID */
4002  if (clamp_srid(gserialized_get_srid(gser)) != srid) {
4003  elog(NOTICE, "Geometry provided for geomval at index %d does not have the same SRID as the raster: %d. Returning original raster", i, srid);
4005  rt_raster_destroy(raster);
4006  PG_RETURN_POINTER(pgraster);
4007  }
4008 
4009  /* Get a 2D version of the geometry if necessary */
4010  if (lwgeom_ndims(arg->gv[arg->ngv].geom) > 2) {
4011  LWGEOM *geom2d = lwgeom_force_2d(arg->gv[arg->ngv].geom);
4012  lwgeom_free(arg->gv[arg->ngv].geom);
4013  arg->gv[arg->ngv].geom = geom2d;
4014  }
4015 
4016  /* filter for types */
4017  gtype = gserialized_get_type(gser);
4018 
4019  /* shortcuts for POINT and MULTIPOINT */
4020  if (gtype == POINTTYPE || gtype == MULTIPOINTTYPE)
4021  allpoint++;
4022 
4023  /* get wkb of geometry */
4024  POSTGIS_RT_DEBUG(3, "getting wkb of geometry");
4025  wkb = lwgeom_to_wkb(arg->gv[arg->ngv].geom, WKB_SFSQL, &wkb_len);
4026 
4027  /* rasterize geometry */
4028  arg->gv[arg->ngv].mask = rt_raster_gdal_rasterize(
4029  wkb, wkb_len,
4030  NULL,
4031  0, NULL,
4032  NULL, NULL,
4033  NULL, NULL,
4034  NULL, NULL,
4035  &(gt[1]), &(gt[5]),
4036  NULL, NULL,
4037  &(gt[0]), &(gt[3]),
4038  &(gt[2]), &(gt[4]),
4039  NULL
4040  );
4041 
4042  pfree(wkb);
4043  if (gtype != POINTTYPE && gtype != MULTIPOINTTYPE) {
4044  lwgeom_free(arg->gv[arg->ngv].geom);
4045  arg->gv[arg->ngv].geom = NULL;
4046  }
4047 
4048  if (arg->gv[arg->ngv].mask == NULL) {
4050  rt_raster_destroy(raster);
4051  PG_FREE_IF_COPY(pgraster, 0);
4052  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not rasterize geometry of geomval at index %d", i);
4053  PG_RETURN_NULL();
4054  }
4055 
4056  /* set SRID */
4057  rt_raster_set_srid(arg->gv[arg->ngv].mask, srid);
4058 
4059  /* second element, value */
4060  POSTGIS_RT_DEBUG(4, "Processing second element (val)");
4061  tupv = GetAttributeByName(tup, "val", &isnull);
4062  if (isnull) {
4063  elog(NOTICE, "Second argument (val) of geomval at index %d is NULL. Treating as NODATA", i);
4064  arg->gv[arg->ngv].pixval.nodata = 1;
4065  }
4066  else
4067  arg->gv[arg->ngv].pixval.value = DatumGetFloat8(tupv);
4068 
4069  (arg->ngv)++;
4070  }
4071 
4072  /* redim arg->gv if needed */
4073  if (arg->ngv < n) {
4074  arg->gv = repalloc(arg->gv, sizeof(struct rtpg_setvaluesgv_geomval_t) * arg->ngv);
4075  if (arg->gv == NULL) {
4077  rt_raster_destroy(raster);
4078  PG_FREE_IF_COPY(pgraster, 0);
4079  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not reallocate memory for geomval array");
4080  PG_RETURN_NULL();
4081  }
4082  }
4083 
4084  /* keepnodata */
4085  if (!PG_ARGISNULL(3))
4086  arg->keepnodata = PG_GETARG_BOOL(3);
4087  POSTGIS_RT_DEBUGF(3, "keepnodata = %d", arg->keepnodata);
4088 
4089  /* keepnodata = TRUE and band is NODATA */
4090  if (arg->keepnodata && rt_band_get_isnodata_flag(band)) {
4091  POSTGIS_RT_DEBUG(3, "keepnodata = TRUE and band is NODATA. Not doing anything");
4092  }
4093  /* all elements are points */
4094  else if (allpoint == arg->ngv) {
4095  double igt[6] = {0};
4096  double xy[2] = {0};
4097  double value = 0;
4098  int isnodata = 0;
4099 
4100  LWCOLLECTION *coll = NULL;
4101  LWPOINT *point = NULL;
4102  POINT2D p;
4103 
4104  POSTGIS_RT_DEBUG(3, "all geometries are points, using direct to pixel method");
4105 
4106  /* cache inverse gretransform matrix */
4108 
4109  /* process each geometry */
4110  for (i = 0; i < arg->ngv; i++) {
4111  /* convert geometry to collection */
4112  coll = lwgeom_as_lwcollection(lwgeom_as_multi(arg->gv[i].geom));
4113 
4114  /* process each point in collection */
4115  for (j = 0; j < coll->ngeoms; j++) {
4116  point = lwgeom_as_lwpoint(coll->geoms[j]);
4117  getPoint2d_p(point->point, 0, &p);
4118 
4119  if (rt_raster_geopoint_to_cell(raster, p.x, p.y, &(xy[0]), &(xy[1]), igt) != ES_NONE) {
4121  rt_raster_destroy(raster);
4122  PG_FREE_IF_COPY(pgraster, 0);
4123  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not process coordinates of point");
4124  PG_RETURN_NULL();
4125  }
4126 
4127  /* skip point if outside raster */
4128  if (
4129  (xy[0] < 0 || xy[0] >= width) ||
4130  (xy[1] < 0 || xy[1] >= height)
4131  ) {
4132  elog(NOTICE, "Point is outside raster extent. Skipping");
4133  continue;
4134  }
4135 
4136  /* get pixel value */
4137  if (rt_band_get_pixel(band, xy[0], xy[1], &value, &isnodata) != ES_NONE) {
4139  rt_raster_destroy(raster);
4140  PG_FREE_IF_COPY(pgraster, 0);
4141  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not get pixel value");
4142  PG_RETURN_NULL();
4143  }
4144 
4145  /* keepnodata = TRUE AND pixel value is NODATA */
4146  if (arg->keepnodata && isnodata)
4147  continue;
4148 
4149  /* set pixel */
4150  if (arg->gv[i].pixval.nodata)
4151  noerr = rt_band_set_pixel(band, xy[0], xy[1], nodataval, NULL);
4152  else
4153  noerr = rt_band_set_pixel(band, xy[0], xy[1], arg->gv[i].pixval.value, NULL);
4154 
4155  if (noerr != ES_NONE) {
4157  rt_raster_destroy(raster);
4158  PG_FREE_IF_COPY(pgraster, 0);
4159  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not set pixel value");
4160  PG_RETURN_NULL();
4161  }
4162  }
4163  }
4164  }
4165  /* run iterator otherwise */
4166  else {
4167  rt_iterator itrset;
4168 
4169  POSTGIS_RT_DEBUG(3, "a mix of geometries, using iterator method");
4170 
4171  /* init itrset */
4172  itrset = palloc(sizeof(struct rt_iterator_t) * (arg->ngv + 1));
4173  if (itrset == NULL) {
4175  rt_raster_destroy(raster);
4176  PG_FREE_IF_COPY(pgraster, 0);
4177  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not allocate memory for iterator arguments");
4178  PG_RETURN_NULL();
4179  }
4180 
4181  /* set first raster's info */
4182  itrset[0].raster = raster;
4183  itrset[0].nband = nband - 1;
4184  itrset[0].nbnodata = 1;
4185 
4186  /* set other raster's info */
4187  for (i = 0, j = 1; i < arg->ngv; i++, j++) {
4188  itrset[j].raster = arg->gv[i].mask;
4189  itrset[j].nband = 0;
4190  itrset[j].nbnodata = 1;
4191  }
4192 
4193  /* pass to iterator */
4194  noerr = rt_raster_iterator(
4195  itrset, arg->ngv + 1,
4196  ET_FIRST, NULL,
4197  pixtype,
4198  hasnodata, nodataval,
4199  0, 0,
4200  arg,
4202  &_raster
4203  );
4204  pfree(itrset);
4205 
4206  if (noerr != ES_NONE) {
4208  rt_raster_destroy(raster);
4209  PG_FREE_IF_COPY(pgraster, 0);
4210  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not run raster iterator function");
4211  PG_RETURN_NULL();
4212  }
4213 
4214  /* copy band from _raster to raster */
4215  _band = rt_raster_get_band(_raster, 0);
4216  if (_band == NULL) {
4218  rt_raster_destroy(_raster);
4219  rt_raster_destroy(raster);
4220  PG_FREE_IF_COPY(pgraster, 0);
4221  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not get band from working raster");
4222  PG_RETURN_NULL();
4223  }
4224 
4225  _band = rt_raster_replace_band(raster, _band, nband - 1);
4226  if (_band == NULL) {
4228  rt_raster_destroy(_raster);
4229  rt_raster_destroy(raster);
4230  PG_FREE_IF_COPY(pgraster, 0);
4231  elog(ERROR, "RASTER_setPixelValuesGeomval: Could not replace band in output raster");
4232  PG_RETURN_NULL();
4233  }
4234 
4235  rt_band_destroy(_band);
4236  rt_raster_destroy(_raster);
4237  }
4238 
4240 
4241  pgrtn = rt_raster_serialize(raster);
4242  rt_raster_destroy(raster);
4243  PG_FREE_IF_COPY(pgraster, 0);
4244 
4245  POSTGIS_RT_DEBUG(3, "Finished");
4246 
4247  if (!pgrtn)
4248  PG_RETURN_NULL();
4249 
4250  SET_VARSIZE(pgrtn, pgrtn->size);
4251  PG_RETURN_POINTER(pgrtn);
4252 }
4253 
4258 Datum RASTER_getPixelPolygons(PG_FUNCTION_ARGS)
4259 {
4260  FuncCallContext *funcctx;
4261  TupleDesc tupdesc;
4262  rt_pixel pix = NULL;
4263  rt_pixel pix2;
4264  int call_cntr;
4265  int max_calls;
4266  int i = 0;
4267 
4268  /* stuff done only on the first call of the function */
4269  if (SRF_IS_FIRSTCALL()) {
4270  MemoryContext oldcontext;
4271  rt_pgraster *pgraster = NULL;
4272  rt_raster raster = NULL;
4273  rt_band band = NULL;
4274  int nband = 1;
4275  int numbands;
4276  bool hasband = TRUE;
4277  bool exclude_nodata_value = TRUE;
4278  bool nocolumnx = FALSE;
4279  bool norowy = FALSE;
4280  int x = 0;
4281  int y = 0;
4282  int bounds[4] = {0};
4283  int pixcount = 0;
4284  double value = 0;
4285  int isnodata = 0;
4286 
4287  LWPOLY *poly;
4288 
4289  POSTGIS_RT_DEBUG(3, "RASTER_getPixelPolygons first call");
4290 
4291  /* create a function context for cross-call persistence */
4292  funcctx = SRF_FIRSTCALL_INIT();
4293 
4294  /* switch to memory context appropriate for multiple function calls */
4295  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
4296 
4297  if (PG_ARGISNULL(0)) {
4298  MemoryContextSwitchTo(oldcontext);
4299  SRF_RETURN_DONE(funcctx);
4300  }
4301  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4302 
4303  /* band */
4304  if (PG_ARGISNULL(1))
4305  hasband = FALSE;
4306  else {
4307  nband = PG_GETARG_INT32(1);
4308  hasband = TRUE;
4309  }
4310 
4311  /* column */
4312  if (PG_ARGISNULL(2))
4313  nocolumnx = TRUE;
4314  else {
4315  bounds[0] = PG_GETARG_INT32(2);
4316  bounds[1] = bounds[0];
4317  }
4318 
4319  /* row */
4320  if (PG_ARGISNULL(3))
4321  norowy = TRUE;
4322  else {
4323  bounds[2] = PG_GETARG_INT32(3);
4324  bounds[3] = bounds[2];
4325  }
4326 
4327  /* exclude NODATA */
4328  if (!PG_ARGISNULL(4))
4329  exclude_nodata_value = PG_GETARG_BOOL(4);
4330 
4331  raster = rt_raster_deserialize(pgraster, FALSE);
4332  if (!raster) {
4333  PG_FREE_IF_COPY(pgraster, 0);
4334  ereport(ERROR, (
4335  errcode(ERRCODE_OUT_OF_MEMORY),
4336  errmsg("Could not deserialize raster")
4337  ));
4338  MemoryContextSwitchTo(oldcontext);
4339  SRF_RETURN_DONE(funcctx);
4340  }
4341 
4342  /* raster empty, return NULL */
4343  if (rt_raster_is_empty(raster)) {
4344  elog(NOTICE, "Raster is empty. Returning NULL");
4345  rt_raster_destroy(raster);
4346  PG_FREE_IF_COPY(pgraster, 0);
4347  MemoryContextSwitchTo(oldcontext);
4348  SRF_RETURN_DONE(funcctx);
4349  }
4350 
4351  /* band specified, load band and info */
4352  if (hasband) {
4353  numbands = rt_raster_get_num_bands(raster);
4354  POSTGIS_RT_DEBUGF(3, "band %d", nband);
4355  POSTGIS_RT_DEBUGF(3, "# of bands %d", numbands);
4356 
4357  if (nband < 1 || nband > numbands) {
4358  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
4359  rt_raster_destroy(raster);
4360  PG_FREE_IF_COPY(pgraster, 0);
4361  MemoryContextSwitchTo(oldcontext);
4362  SRF_RETURN_DONE(funcctx);
4363  }
4364 
4365  band = rt_raster_get_band(raster, nband - 1);
4366  if (!band) {
4367  elog(NOTICE, "Could not find band at index %d. Returning NULL", nband);
4368  rt_raster_destroy(raster);
4369  PG_FREE_IF_COPY(pgraster, 0);
4370  MemoryContextSwitchTo(oldcontext);
4371  SRF_RETURN_DONE(funcctx);
4372  }
4373 
4374  if (!rt_band_get_hasnodata_flag(band))
4375  exclude_nodata_value = FALSE;
4376  }
4377 
4378  /* set bounds if columnx, rowy not set */
4379  if (nocolumnx) {
4380  bounds[0] = 1;
4381  bounds[1] = rt_raster_get_width(raster);
4382  }
4383  if (norowy) {
4384  bounds[2] = 1;
4385  bounds[3] = rt_raster_get_height(raster);
4386  }
4387  POSTGIS_RT_DEBUGF(3, "bounds (min x, max x, min y, max y) = (%d, %d, %d, %d)",
4388  bounds[0], bounds[1], bounds[2], bounds[3]);
4389 
4390  /* rowy */
4391  pixcount = 0;
4392  for (y = bounds[2]; y <= bounds[3]; y++) {
4393  /* columnx */
4394  for (x = bounds[0]; x <= bounds[1]; x++) {
4395 
4396  value = 0;
4397  isnodata = TRUE;
4398 
4399  if (hasband) {
4400  if (rt_band_get_pixel(band, x - 1, y - 1, &value, &isnodata) != ES_NONE) {
4401 
4402  for (i = 0; i < pixcount; i++)
4403  lwgeom_free(pix[i].geom);
4404  if (pixcount) pfree(pix);
4405 
4406  rt_band_destroy(band);
4407  rt_raster_destroy(raster);
4408  PG_FREE_IF_COPY(pgraster, 0);
4409 
4410  MemoryContextSwitchTo(oldcontext);
4411  elog(ERROR, "RASTER_getPixelPolygons: Could not get pixel value");
4412  SRF_RETURN_DONE(funcctx);
4413  }
4414 
4415  /* don't continue if pixel is NODATA and to exclude NODATA */
4416  if (isnodata && exclude_nodata_value) {
4417  POSTGIS_RT_DEBUG(5, "pixel value is NODATA and exclude_nodata_value = TRUE");
4418  continue;
4419  }
4420  }
4421 
4422  /* geometry */
4423  poly = rt_raster_pixel_as_polygon(raster, x - 1, y - 1);
4424  if (!poly) {
4425  for (i = 0; i < pixcount; i++)
4426  lwgeom_free(pix[i].geom);
4427  if (pixcount) pfree(pix);
4428 
4429  if (hasband) rt_band_destroy(band);
4430  rt_raster_destroy(raster);
4431  PG_FREE_IF_COPY(pgraster, 0);
4432 
4433  MemoryContextSwitchTo(oldcontext);
4434  elog(ERROR, "RASTER_getPixelPolygons: Could not get pixel polygon");
4435  SRF_RETURN_DONE(funcctx);
4436  }
4437 
4438  if (!pixcount)
4439  pix = palloc(sizeof(struct rt_pixel_t) * (pixcount + 1));
4440  else
4441  pix = repalloc(pix, sizeof(struct rt_pixel_t) * (pixcount + 1));
4442  if (pix == NULL) {
4443 
4444  lwpoly_free(poly);
4445  if (hasband) rt_band_destroy(band);
4446  rt_raster_destroy(raster);
4447  PG_FREE_IF_COPY(pgraster, 0);
4448 
4449  MemoryContextSwitchTo(oldcontext);
4450  elog(ERROR, "RASTER_getPixelPolygons: Could not allocate memory for storing pixel polygons");
4451  SRF_RETURN_DONE(funcctx);
4452  }
4453  pix[pixcount].geom = (LWGEOM *) poly;
4454  POSTGIS_RT_DEBUGF(5, "poly @ %p", poly);
4455  POSTGIS_RT_DEBUGF(5, "geom @ %p", pix[pixcount].geom);
4456 
4457  /* x, y */
4458  pix[pixcount].x = x;
4459  pix[pixcount].y = y;
4460 
4461  /* value */
4462  pix[pixcount].value = value;
4463 
4464  /* NODATA */
4465  if (hasband) {
4466  if (exclude_nodata_value)
4467  pix[pixcount].nodata = isnodata;
4468  else
4469  pix[pixcount].nodata = FALSE;
4470  }
4471  else {
4472  pix[pixcount].nodata = isnodata;
4473  }
4474 
4475  pixcount++;
4476  }
4477  }
4478 
4479  if (hasband) rt_band_destroy(band);
4480  rt_raster_destroy(raster);
4481  PG_FREE_IF_COPY(pgraster, 0);
4482 
4483  /* shortcut if no pixcount */
4484  if (pixcount < 1) {
4485  elog(NOTICE, "No pixels found for band %d", nband);
4486  MemoryContextSwitchTo(oldcontext);
4487  SRF_RETURN_DONE(funcctx);
4488  }
4489 
4490  /* Store needed information */
4491  funcctx->user_fctx = pix;
4492 
4493  /* total number of tuples to be returned */
4494  funcctx->max_calls = pixcount;
4495  POSTGIS_RT_DEBUGF(3, "pixcount = %d", pixcount);
4496 
4497  /* Build a tuple descriptor for our result type */
4498  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
4499  ereport(ERROR, (
4500  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
4501  errmsg("function returning record called in context that cannot accept type record")
4502  ));
4503  }
4504 
4505  BlessTupleDesc(tupdesc);
4506  funcctx->tuple_desc = tupdesc;
4507 
4508  MemoryContextSwitchTo(oldcontext);
4509  }
4510 
4511  /* stuff done on every call of the function */
4512  funcctx = SRF_PERCALL_SETUP();
4513 
4514  call_cntr = funcctx->call_cntr;
4515  max_calls = funcctx->max_calls;
4516  tupdesc = funcctx->tuple_desc;
4517  pix2 = funcctx->user_fctx;
4518 
4519  /* do when there is more left to send */
4520  if (call_cntr < max_calls) {
4521  int values_length = 4;
4522  Datum values[values_length];
4523  bool nulls[values_length];
4524  HeapTuple tuple;
4525  Datum result;
4526 
4527  GSERIALIZED *gser = NULL;
4528  size_t gser_size = 0;
4529 
4530  POSTGIS_RT_DEBUGF(3, "call number %d", call_cntr);
4531 
4532  memset(nulls, FALSE, sizeof(bool) * values_length);
4533 
4534  /* convert LWGEOM to GSERIALIZED */
4535  gser = gserialized_from_lwgeom(pix2[call_cntr].geom, 0, &gser_size);
4536  lwgeom_free(pix2[call_cntr].geom);
4537 
4538  values[0] = PointerGetDatum(gser);
4539  if (pix2[call_cntr].nodata)
4540  nulls[1] = TRUE;
4541  else
4542  values[1] = Float8GetDatum(pix2[call_cntr].value);
4543  values[2] = Int32GetDatum(pix2[call_cntr].x);
4544  values[3] = Int32GetDatum(pix2[call_cntr].y);
4545 
4546  /* build a tuple */
4547  tuple = heap_form_tuple(tupdesc, values, nulls);
4548 
4549  /* make the tuple into a datum */
4550  result = HeapTupleGetDatum(tuple);
4551 
4552  SRF_RETURN_NEXT(funcctx, result);
4553  }
4554  /* do when there is no more left */
4555  else {
4556  pfree(pix2);
4557  SRF_RETURN_DONE(funcctx);
4558  }
4559 }
4560 
4565 Datum RASTER_getPolygon(PG_FUNCTION_ARGS)
4566 {
4567  rt_pgraster *pgraster = NULL;
4568  rt_raster raster = NULL;
4569  int num_bands = 0;
4570  int nband = 1;
4571  int err;
4572  LWMPOLY *surface = NULL;
4573  GSERIALIZED *rtn = NULL;
4574 
4575  /* raster */
4576  if (PG_ARGISNULL(0))
4577  PG_RETURN_NULL();
4578  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4579 
4580  raster = rt_raster_deserialize(pgraster, FALSE);
4581  if (!raster) {
4582  PG_FREE_IF_COPY(pgraster, 0);
4583  elog(ERROR, "RASTER_getPolygon: Could not deserialize raster");
4584  PG_RETURN_NULL();
4585  }
4586 
4587  /* num_bands */
4588  num_bands = rt_raster_get_num_bands(raster);
4589  if (num_bands < 1) {
4590  elog(NOTICE, "Raster provided has no bands");
4591  rt_raster_destroy(raster);
4592  PG_FREE_IF_COPY(pgraster, 0);
4593  PG_RETURN_NULL();
4594  }
4595 
4596  /* band index is 1-based */
4597  if (!PG_ARGISNULL(1))
4598  nband = PG_GETARG_INT32(1);
4599  if (nband < 1 || nband > num_bands) {
4600  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
4601  rt_raster_destroy(raster);
4602  PG_FREE_IF_COPY(pgraster, 0);
4603  PG_RETURN_NULL();
4604  }
4605 
4606  /* get band surface */
4607  err = rt_raster_surface(raster, nband - 1, &surface);
4608  rt_raster_destroy(raster);
4609  PG_FREE_IF_COPY(pgraster, 0);
4610 
4611  if (err != ES_NONE) {
4612  elog(ERROR, "RASTER_getPolygon: Could not get raster band's surface");
4613  PG_RETURN_NULL();
4614  }
4615  else if (surface == NULL) {
4616  elog(NOTICE, "Raster is empty or all pixels of band are NODATA. Returning NULL");
4617  PG_RETURN_NULL();
4618  }
4619 
4620  rtn = geometry_serialize(lwmpoly_as_lwgeom(surface));
4621  lwmpoly_free(surface);
4622 
4623  PG_RETURN_POINTER(rtn);
4624 }
4625 
4630 Datum RASTER_pixelOfValue(PG_FUNCTION_ARGS)
4631 {
4632  FuncCallContext *funcctx;
4633  TupleDesc tupdesc;
4634 
4635  rt_pixel pixels = NULL;
4636  rt_pixel pixels2 = NULL;
4637  int count = 0;
4638  int i = 0;
4639  int n = 0;
4640  int call_cntr;
4641  int max_calls;
4642 
4643  if (SRF_IS_FIRSTCALL()) {
4644  MemoryContext oldcontext;
4645 
4646  rt_pgraster *pgraster = NULL;
4647  rt_raster raster = NULL;
4648  rt_band band = NULL;
4649  int nband = 1;
4650  int num_bands = 0;
4651  double *search = NULL;
4652  int nsearch = 0;
4653  double val;
4654  bool exclude_nodata_value = TRUE;
4655 
4656  ArrayType *array;
4657  Oid etype;
4658  Datum *e;
4659  bool *nulls;
4660  int16 typlen;
4661  bool typbyval;
4662  char typalign;
4663 
4664  /* create a function context for cross-call persistence */
4665  funcctx = SRF_FIRSTCALL_INIT();
4666 
4667  /* switch to memory context appropriate for multiple function calls */
4668  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
4669 
4670  if (PG_ARGISNULL(0)) {
4671  MemoryContextSwitchTo(oldcontext);
4672  SRF_RETURN_DONE(funcctx);
4673  }
4674  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4675  raster = rt_raster_deserialize(pgraster, FALSE);
4676  if (!raster) {
4677  PG_FREE_IF_COPY(pgraster, 0);
4678  MemoryContextSwitchTo(oldcontext);
4679  elog(ERROR, "RASTER_pixelOfValue: Could not deserialize raster");
4680  SRF_RETURN_DONE(funcctx);
4681  }
4682 
4683  /* num_bands */
4684  num_bands = rt_raster_get_num_bands(raster);
4685  if (num_bands < 1) {
4686  elog(NOTICE, "Raster provided has no bands");
4687  rt_raster_destroy(raster);
4688  PG_FREE_IF_COPY(pgraster, 0);
4689  MemoryContextSwitchTo(oldcontext);
4690  SRF_RETURN_DONE(funcctx);
4691  }
4692 
4693  /* band index is 1-based */
4694  if (!PG_ARGISNULL(1))
4695  nband = PG_GETARG_INT32(1);
4696  if (nband < 1 || nband > num_bands) {
4697  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
4698  rt_raster_destroy(raster);
4699  PG_FREE_IF_COPY(pgraster, 0);
4700  MemoryContextSwitchTo(oldcontext);
4701  SRF_RETURN_DONE(funcctx);
4702  }
4703 
4704  /* search values */
4705  array = PG_GETARG_ARRAYTYPE_P(2);
4706  etype = ARR_ELEMTYPE(array);
4707  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
4708 
4709  switch (etype) {
4710  case FLOAT4OID:
4711  case FLOAT8OID:
4712  break;
4713  default:
4714  rt_raster_destroy(raster);
4715  PG_FREE_IF_COPY(pgraster, 0);
4716  MemoryContextSwitchTo(oldcontext);
4717  elog(ERROR, "RASTER_pixelOfValue: Invalid data type for pixel values");
4718  SRF_RETURN_DONE(funcctx);
4719  break;
4720  }
4721 
4722  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
4723  &nulls, &n);
4724 
4725  search = palloc(sizeof(double) * n);
4726  for (i = 0, nsearch = 0; i < n; i++) {
4727  if (nulls[i]) continue;
4728 
4729  switch (etype) {
4730  case FLOAT4OID:
4731  val = (double) DatumGetFloat4(e[i]);
4732  break;
4733  case FLOAT8OID:
4734  val = (double) DatumGetFloat8(e[i]);
4735  break;
4736  }
4737 
4738  search[nsearch] = val;
4739  POSTGIS_RT_DEBUGF(3, "search[%d] = %f", nsearch, search[nsearch]);
4740  nsearch++;
4741  }
4742 
4743  /* not searching for anything */
4744  if (nsearch < 1) {
4745  elog(NOTICE, "No search values provided. Returning NULL");
4746  pfree(search);
4747  rt_raster_destroy(raster);
4748  PG_FREE_IF_COPY(pgraster, 0);
4749  MemoryContextSwitchTo(oldcontext);
4750  SRF_RETURN_DONE(funcctx);
4751  }
4752  else if (nsearch < n)
4753  search = repalloc(search, sizeof(double) * nsearch);
4754 
4755  /* exclude_nodata_value flag */
4756  if (!PG_ARGISNULL(3))
4757  exclude_nodata_value = PG_GETARG_BOOL(3);
4758 
4759  /* get band */
4760  band = rt_raster_get_band(raster, nband - 1);
4761  if (!band) {
4762  elog(NOTICE, "Could not find band at index %d. Returning NULL", nband);
4763  rt_raster_destroy(raster);
4764  PG_FREE_IF_COPY(pgraster, 0);
4765  MemoryContextSwitchTo(oldcontext);
4766  SRF_RETURN_DONE(funcctx);
4767  }
4768 
4769  /* get pixels of values */
4771  band, exclude_nodata_value,
4772  search, nsearch,
4773  &pixels
4774  );
4775  pfree(search);
4776  rt_band_destroy(band);
4777  rt_raster_destroy(raster);
4778  PG_FREE_IF_COPY(pgraster, 0);
4779  if (count < 1) {
4780  /* error */
4781  if (count < 0)
4782  elog(NOTICE, "Could not get the pixels of search values for band at index %d", nband);
4783  /* no nearest pixel */
4784  else
4785  elog(NOTICE, "No pixels of search values found for band at index %d", nband);
4786 
4787  MemoryContextSwitchTo(oldcontext);
4788  SRF_RETURN_DONE(funcctx);
4789  }
4790 
4791  /* Store needed information */
4792  funcctx->user_fctx = pixels;
4793 
4794  /* total number of tuples to be returned */
4795  funcctx->max_calls = count;
4796 
4797  /* Build a tuple descriptor for our result type */
4798  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
4799  ereport(ERROR, (
4800  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
4801  errmsg(
4802  "function returning record called in context "
4803  "that cannot accept type record"
4804  )
4805  ));
4806  }
4807 
4808  BlessTupleDesc(tupdesc);
4809  funcctx->tuple_desc = tupdesc;
4810 
4811  MemoryContextSwitchTo(oldcontext);
4812  }
4813 
4814  /* stuff done on every call of the function */
4815  funcctx = SRF_PERCALL_SETUP();
4816 
4817  call_cntr = funcctx->call_cntr;
4818  max_calls = funcctx->max_calls;
4819  tupdesc = funcctx->tuple_desc;
4820  pixels2 = funcctx->user_fctx;
4821 
4822  /* do when there is more left to send */
4823  if (call_cntr < max_calls) {
4824  int values_length = 3;
4825  Datum values[values_length];
4826  bool nulls[values_length];
4827  HeapTuple tuple;
4828  Datum result;
4829 
4830  memset(nulls, FALSE, sizeof(bool) * values_length);
4831 
4832  /* 0-based to 1-based */
4833  pixels2[call_cntr].x += 1;
4834  pixels2[call_cntr].y += 1;
4835 
4836  values[0] = Float8GetDatum(pixels2[call_cntr].value);
4837  values[1] = Int32GetDatum(pixels2[call_cntr].x);
4838  values[2] = Int32GetDatum(pixels2[call_cntr].y);
4839 
4840  /* build a tuple */
4841  tuple = heap_form_tuple(tupdesc, values, nulls);
4842 
4843  /* make the tuple into a datum */
4844  result = HeapTupleGetDatum(tuple);
4845 
4846  SRF_RETURN_NEXT(funcctx, result);
4847  }
4848  else {
4849  pfree(pixels2);
4850  SRF_RETURN_DONE(funcctx);
4851  }
4852 }
4853 
4858 Datum RASTER_nearestValue(PG_FUNCTION_ARGS)
4859 {
4860  rt_pgraster *pgraster = NULL;
4861  rt_raster raster = NULL;
4862  rt_band band = NULL;
4863  int bandindex = 1;
4864  int num_bands = 0;
4865  GSERIALIZED *geom;
4866  bool exclude_nodata_value = TRUE;
4867  LWGEOM *lwgeom;
4868  LWPOINT *point = NULL;
4869  POINT2D p;
4870 
4871  double x;
4872  double y;
4873  int count;
4874  rt_pixel npixels = NULL;
4875  double value = 0;
4876  int hasvalue = 0;
4877  int isnodata = 0;
4878 
4879  if (PG_ARGISNULL(0))
4880  PG_RETURN_NULL();
4881  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4882  raster = rt_raster_deserialize(pgraster, FALSE);
4883  if (!raster) {
4884  PG_FREE_IF_COPY(pgraster, 0);
4885  elog(ERROR, "RASTER_nearestValue: Could not deserialize raster");
4886  PG_RETURN_NULL();
4887  }
4888 
4889  /* band index is 1-based */
4890  if (!PG_ARGISNULL(1))
4891  bandindex = PG_GETARG_INT32(1);
4892  num_bands = rt_raster_get_num_bands(raster);
4893  if (bandindex < 1 || bandindex > num_bands) {
4894  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
4895  rt_raster_destroy(raster);
4896  PG_FREE_IF_COPY(pgraster, 0);
4897  PG_RETURN_NULL();
4898  }
4899 
4900  /* point */
4901  geom = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(2));
4902  if (gserialized_get_type(geom) != POINTTYPE) {
4903  elog(NOTICE, "Geometry provided must be a point");
4904  rt_raster_destroy(raster);
4905  PG_FREE_IF_COPY(pgraster, 0);
4906  PG_FREE_IF_COPY(geom, 2);
4907  PG_RETURN_NULL();
4908  }
4909 
4910  /* exclude_nodata_value flag */
4911  if (!PG_ARGISNULL(3))
4912  exclude_nodata_value = PG_GETARG_BOOL(3);
4913 
4914  /* SRIDs of raster and geometry must match */
4916  elog(NOTICE, "SRIDs of geometry and raster do not match");
4917  rt_raster_destroy(raster);
4918  PG_FREE_IF_COPY(pgraster, 0);
4919  PG_FREE_IF_COPY(geom, 2);
4920  PG_RETURN_NULL();
4921  }
4922 
4923  /* get band */
4924  band = rt_raster_get_band(raster, bandindex - 1);
4925  if (!band) {
4926  elog(NOTICE, "Could not find band at index %d. Returning NULL", bandindex);
4927  rt_raster_destroy(raster);
4928  PG_FREE_IF_COPY(pgraster, 0);
4929  PG_FREE_IF_COPY(geom, 2);
4930  PG_RETURN_NULL();
4931  }
4932 
4933  /* process geometry */
4934  lwgeom = lwgeom_from_gserialized(geom);
4935 
4936  if (lwgeom_is_empty(lwgeom)) {
4937  elog(NOTICE, "Geometry provided cannot be empty");
4938  rt_raster_destroy(raster);
4939  PG_FREE_IF_COPY(pgraster, 0);
4940  PG_FREE_IF_COPY(geom, 2);
4941  PG_RETURN_NULL();
4942  }
4943 
4944  /* Get a 2D version of the geometry if necessary */
4945  if (lwgeom_ndims(lwgeom) > 2) {
4946  LWGEOM *lwgeom2d = lwgeom_force_2d(lwgeom);
4947  lwgeom_free(lwgeom);
4948  lwgeom = lwgeom2d;
4949  }
4950 
4951  point = lwgeom_as_lwpoint(lwgeom);
4952  getPoint2d_p(point->point, 0, &p);
4953 
4955  raster,
4956  p.x, p.y,
4957  &x, &y,
4958  NULL
4959  ) != ES_NONE) {
4960  rt_raster_destroy(raster);
4961  PG_FREE_IF_COPY(pgraster, 0);
4962  lwgeom_free(lwgeom);
4963  PG_FREE_IF_COPY(geom, 2);
4964  elog(ERROR, "RASTER_nearestValue: Could not compute pixel coordinates from spatial coordinates");
4965  PG_RETURN_NULL();
4966  }
4967 
4968  /* get value at point */
4969  if (
4970  (x >= 0 && x < rt_raster_get_width(raster)) &&
4971  (y >= 0 && y < rt_raster_get_height(raster))
4972  ) {
4973  if (rt_band_get_pixel(band, x, y, &value, &isnodata) != ES_NONE) {
4974  rt_raster_destroy(raster);
4975  PG_FREE_IF_COPY(pgraster, 0);
4976  lwgeom_free(lwgeom);
4977  PG_FREE_IF_COPY(geom, 2);
4978  elog(ERROR, "RASTER_nearestValue: Could not get pixel value for band at index %d", bandindex);
4979  PG_RETURN_NULL();
4980  }
4981 
4982  /* value at point, return value */
4983  if (!exclude_nodata_value || !isnodata) {
4984  rt_raster_destroy(raster);
4985  PG_FREE_IF_COPY(pgraster, 0);
4986  lwgeom_free(lwgeom);
4987  PG_FREE_IF_COPY(geom, 2);
4988 
4989  PG_RETURN_FLOAT8(value);
4990  }
4991  }
4992 
4993  /* get neighborhood */
4994  count = rt_band_get_nearest_pixel(
4995  band,
4996  x, y,
4997  0, 0,
4998  exclude_nodata_value,
4999  &npixels
5000  );
5001  rt_band_destroy(band);
5002  /* error or no neighbors */
5003  if (count < 1) {
5004  /* error */
5005  if (count < 0)
5006  elog(NOTICE, "Could not get the nearest value for band at index %d", bandindex);
5007  /* no nearest pixel */
5008  else
5009  elog(NOTICE, "No nearest value found for band at index %d", bandindex);
5010 
5011  lwgeom_free(lwgeom);
5012  PG_FREE_IF_COPY(geom, 2);
5013  rt_raster_destroy(raster);
5014  PG_FREE_IF_COPY(pgraster, 0);
5015  PG_RETURN_NULL();
5016  }
5017 
5018  /* more than one nearest value, see which one is closest */
5019  if (count > 1) {
5020  int i = 0;
5021  LWPOLY *poly = NULL;
5022  double lastdist = -1;
5023  double dist;
5024 
5025  for (i = 0; i < count; i++) {
5026  /* convex-hull of pixel */
5027  poly = rt_raster_pixel_as_polygon(raster, npixels[i].x, npixels[i].y);
5028  if (!poly) {
5029  lwgeom_free(lwgeom);
5030  PG_FREE_IF_COPY(geom, 2);
5031  rt_raster_destroy(raster);
5032  PG_FREE_IF_COPY(pgraster, 0);
5033  elog(ERROR, "RASTER_nearestValue: Could not get polygon of neighboring pixel");
5034  PG_RETURN_NULL();
5035  }
5036 
5037  /* distance between convex-hull and point */
5038  dist = lwgeom_mindistance2d(lwpoly_as_lwgeom(poly), lwgeom);
5039  if (lastdist < 0 || dist < lastdist) {
5040  value = npixels[i].value;
5041  hasvalue = 1;
5042  }
5043  lastdist = dist;
5044 
5045  lwpoly_free(poly);
5046  }
5047  }
5048  else {
5049  value = npixels[0].value;
5050  hasvalue = 1;
5051  }
5052 
5053  pfree(npixels);
5054  lwgeom_free(lwgeom);
5055  PG_FREE_IF_COPY(geom, 2);
5056  rt_raster_destroy(raster);
5057  PG_FREE_IF_COPY(pgraster, 0);
5058 
5059  if (hasvalue)
5060  PG_RETURN_FLOAT8(value);
5061  else
5062  PG_RETURN_NULL();
5063 }
5064 
5069 Datum RASTER_neighborhood(PG_FUNCTION_ARGS)
5070 {
5071  rt_pgraster *pgraster = NULL;
5072  rt_raster raster = NULL;
5073  rt_band band = NULL;
5074  int bandindex = 1;
5075  int num_bands = 0;
5076  int x = 0;
5077  int y = 0;
5078  int _x = 0;
5079  int _y = 0;
5080  int distance[2] = {0};
5081  bool exclude_nodata_value = TRUE;
5082  double pixval;
5083  int isnodata = 0;
5084 
5085  rt_pixel npixels = NULL;
5086  int count;
5087  double **value2D = NULL;
5088  int **nodata2D = NULL;
5089 
5090  int i = 0;
5091  int j = 0;
5092  int k = 0;
5093  Datum *value1D = NULL;
5094  bool *nodata1D = NULL;
5095  int dim[2] = {0};
5096  int lbound[2] = {1, 1};
5097  ArrayType *mdArray = NULL;
5098 
5099  int16 typlen;
5100  bool typbyval;
5101  char typalign;
5102 
5103  /* pgraster is null, return nothing */
5104  if (PG_ARGISNULL(0))
5105  PG_RETURN_NULL();
5106  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5107 
5108  raster = rt_raster_deserialize(pgraster, FALSE);
5109  if (!raster) {
5110  PG_FREE_IF_COPY(pgraster, 0);
5111  elog(ERROR, "RASTER_neighborhood: Could not deserialize raster");
5112  PG_RETURN_NULL();
5113  }
5114 
5115  /* band index is 1-based */
5116  if (!PG_ARGISNULL(1))
5117  bandindex = PG_GETARG_INT32(1);
5118  num_bands = rt_raster_get_num_bands(raster);
5119  if (bandindex < 1 || bandindex > num_bands) {
5120  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
5121  rt_raster_destroy(raster);
5122  PG_FREE_IF_COPY(pgraster, 0);
5123  PG_RETURN_NULL();
5124  }
5125 
5126  /* pixel column, 1-based */
5127  x = PG_GETARG_INT32(2);
5128  _x = x - 1;
5129 
5130  /* pixel row, 1-based */
5131  y = PG_GETARG_INT32(3);
5132  _y = y - 1;
5133 
5134  /* distance X axis */
5135  distance[0] = PG_GETARG_INT32(4);
5136  if (distance[0] < 0) {
5137  elog(NOTICE, "Invalid value for distancex (must be >= zero). Returning NULL");
5138  rt_raster_destroy(raster);
5139  PG_FREE_IF_COPY(pgraster, 0);
5140  PG_RETURN_NULL();
5141  }
5142  distance[0] = (uint16_t) distance[0];
5143 
5144  /* distance Y axis */
5145  distance[1] = PG_GETARG_INT32(5);
5146  if (distance[1] < 0) {
5147  elog(NOTICE, "Invalid value for distancey (must be >= zero). Returning NULL");
5148  rt_raster_destroy(raster);
5149  PG_FREE_IF_COPY(pgraster, 0);
5150  PG_RETURN_NULL();
5151  }
5152  distance[1] = (uint16_t) distance[1];
5153 
5154  /* exclude_nodata_value flag */
5155  if (!PG_ARGISNULL(6))
5156  exclude_nodata_value = PG_GETARG_BOOL(6);
5157 
5158  /* get band */
5159  band = rt_raster_get_band(raster, bandindex - 1);
5160  if (!band) {
5161  elog(NOTICE, "Could not find band at index %d. Returning NULL", bandindex);
5162  rt_raster_destroy(raster);
5163  PG_FREE_IF_COPY(pgraster, 0);
5164  PG_RETURN_NULL();
5165  }
5166 
5167  /* get neighborhood */
5168  count = 0;
5169  npixels = NULL;
5170  if (distance[0] > 0 || distance[1] > 0) {
5171  count = rt_band_get_nearest_pixel(
5172  band,
5173  _x, _y,
5174  distance[0], distance[1],
5175  exclude_nodata_value,
5176  &npixels
5177  );
5178  /* error */
5179  if (count < 0) {
5180  elog(NOTICE, "Could not get the pixel's neighborhood for band at index %d", bandindex);
5181 
5182  rt_band_destroy(band);
5183  rt_raster_destroy(raster);
5184  PG_FREE_IF_COPY(pgraster, 0);
5185 
5186  PG_RETURN_NULL();
5187  }
5188  }
5189 
5190  /* get pixel's value */
5191  if (
5192  (_x >= 0 && _x < rt_band_get_width(band)) &&
5193  (_y >= 0 && _y < rt_band_get_height(band))
5194  ) {
5195  if (rt_band_get_pixel(
5196  band,
5197  _x, _y,
5198  &pixval,
5199  &isnodata
5200  ) != ES_NONE) {
5201  elog(NOTICE, "Could not get the pixel of band at index %d. Returning NULL", bandindex);
5202  rt_band_destroy(band);
5203  rt_raster_destroy(raster);
5204  PG_FREE_IF_COPY(pgraster, 0);
5205  PG_RETURN_NULL();
5206  }
5207  }
5208  /* outside band extent, set to NODATA */
5209  else {
5210  /* has NODATA, use NODATA */
5211  if (rt_band_get_hasnodata_flag(band))
5212  rt_band_get_nodata(band, &pixval);
5213  /* no NODATA, use min possible value */
5214  else
5215  pixval = rt_band_get_min_value(band);
5216  isnodata = 1;
5217  }
5218  POSTGIS_RT_DEBUGF(4, "pixval: %f", pixval);
5219 
5220 
5221  /* add pixel to neighborhood */
5222  count++;
5223  if (count > 1)
5224  npixels = (rt_pixel) repalloc(npixels, sizeof(struct rt_pixel_t) * count);
5225  else
5226  npixels = (rt_pixel) palloc(sizeof(struct rt_pixel_t));
5227  if (npixels == NULL) {
5228 
5229  rt_band_destroy(band);
5230  rt_raster_destroy(raster);
5231  PG_FREE_IF_COPY(pgraster, 0);
5232 
5233  elog(ERROR, "RASTER_neighborhood: Could not reallocate memory for neighborhood");
5234  PG_RETURN_NULL();
5235  }
5236  npixels[count - 1].x = _x;
5237  npixels[count - 1].y = _y;
5238  npixels[count - 1].nodata = 1;
5239  npixels[count - 1].value = pixval;
5240 
5241  /* set NODATA */
5242  if (!exclude_nodata_value || !isnodata) {
5243  npixels[count - 1].nodata = 0;
5244  }
5245 
5246  /* free unnecessary stuff */
5247  rt_band_destroy(band);
5248  rt_raster_destroy(raster);
5249  PG_FREE_IF_COPY(pgraster, 0);
5250 
5251  /* convert set of rt_pixel to 2D array */
5252  /* dim is passed with element 0 being Y-axis and element 1 being X-axis */
5253  count = rt_pixel_set_to_array(
5254  npixels, count,
5255  _x, _y,
5256  distance[0], distance[1],
5257  &value2D,
5258  &nodata2D,
5259  &(dim[1]), &(dim[0])
5260  );
5261  pfree(npixels);
5262  if (count != ES_NONE) {
5263  elog(NOTICE, "Could not create 2D array of neighborhood");
5264  PG_RETURN_NULL();
5265  }
5266 
5267  /* 1D arrays for values and nodata from 2D arrays */
5268  value1D = palloc(sizeof(Datum) * dim[0] * dim[1]);
5269  nodata1D = palloc(sizeof(bool) * dim[0] * dim[1]);
5270 
5271  if (value1D == NULL || nodata1D == NULL) {
5272 
5273  for (i = 0; i < dim[0]; i++) {
5274  pfree(value2D[i]);
5275  pfree(nodata2D[i]);
5276  }
5277  pfree(value2D);
5278  pfree(nodata2D);
5279 
5280  elog(ERROR, "RASTER_neighborhood: Could not allocate memory for return 2D array");
5281  PG_RETURN_NULL();
5282  }
5283 
5284  /* copy values from 2D arrays to 1D arrays */
5285  k = 0;
5286  /* Y-axis */
5287  for (i = 0; i < dim[0]; i++) {
5288  /* X-axis */
5289  for (j = 0; j < dim[1]; j++) {
5290  nodata1D[k] = (bool) nodata2D[i][j];
5291  if (!nodata1D[k])
5292  value1D[k] = Float8GetDatum(value2D[i][j]);
5293  else
5294  value1D[k] = PointerGetDatum(NULL);
5295 
5296  k++;
5297  }
5298  }
5299 
5300  /* no more need for 2D arrays */
5301  for (i = 0; i < dim[0]; i++) {
5302  pfree(value2D[i]);
5303  pfree(nodata2D[i]);
5304  }
5305  pfree(value2D);
5306  pfree(nodata2D);
5307 
5308  /* info about the type of item in the multi-dimensional array (float8). */
5309  get_typlenbyvalalign(FLOAT8OID, &typlen, &typbyval, &typalign);
5310 
5311  mdArray = construct_md_array(
5312  value1D, nodata1D,
5313  2, dim, lbound,
5314  FLOAT8OID,
5315  typlen, typbyval, typalign
5316  );
5317 
5318  pfree(value1D);
5319  pfree(nodata1D);
5320 
5321  PG_RETURN_ARRAYTYPE_P(mdArray);
5322 }
5323 
5328 Datum RASTER_addBand(PG_FUNCTION_ARGS)
5329 {
5330  rt_pgraster *pgraster = NULL;
5331  rt_pgraster *pgrtn = NULL;
5332  rt_raster raster = NULL;
5333  int bandindex = 0;
5334  int maxbandindex = 0;
5335  int numbands = 0;
5336  int lastnumbands = 0;
5337 
5338  text *text_pixtype = NULL;
5339  char *char_pixtype = NULL;
5340 
5341  struct addbandarg {
5342  int index;
5343  bool append;
5344  rt_pixtype pixtype;
5345  double initialvalue;
5346  bool hasnodata;
5347  double nodatavalue;
5348  };
5349  struct addbandarg *arg = NULL;
5350 
5351  ArrayType *array;
5352  Oid etype;
5353  Datum *e;
5354  bool *nulls;
5355  int16 typlen;
5356  bool typbyval;
5357  char typalign;
5358  int n = 0;
5359 
5360  HeapTupleHeader tup;
5361  bool isnull;
5362  Datum tupv;
5363 
5364  int i = 0;
5365 
5366  /* pgraster is null, return null */
5367  if (PG_ARGISNULL(0))
5368  PG_RETURN_NULL();
5369  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5370 
5371  /* raster */
5372  raster = rt_raster_deserialize(pgraster, FALSE);
5373  if (!raster) {
5374  PG_FREE_IF_COPY(pgraster, 0);
5375  elog(ERROR, "RASTER_addBand: Could not deserialize raster");
5376  PG_RETURN_NULL();
5377  }
5378 
5379  /* process set of addbandarg */
5380  POSTGIS_RT_DEBUG(3, "Processing Arg 1 (addbandargset)");
5381  array = PG_GETARG_ARRAYTYPE_P(1);
5382  etype = ARR_ELEMTYPE(array);
5383  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
5384 
5385  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
5386  &nulls, &n);
5387 
5388  if (!n) {
5389  PG_FREE_IF_COPY(pgraster, 0);
5390  elog(ERROR, "RASTER_addBand: Invalid argument for addbandargset");
5391  PG_RETURN_NULL();
5392  }
5393 
5394  /* allocate addbandarg */
5395  arg = (struct addbandarg *) palloc(sizeof(struct addbandarg) * n);
5396  if (arg == NULL) {
5397  rt_raster_destroy(raster);
5398  PG_FREE_IF_COPY(pgraster, 0);
5399  elog(ERROR, "RASTER_addBand: Could not allocate memory for addbandarg");
5400  PG_RETURN_NULL();
5401  }
5402 
5403  /*
5404  process each element of addbandargset
5405  each element is the index of where to add the new band,
5406  new band's pixeltype, the new band's initial value and
5407  the new band's NODATA value if NOT NULL
5408  */
5409  for (i = 0; i < n; i++) {
5410  if (nulls[i]) continue;
5411 
5412  POSTGIS_RT_DEBUGF(4, "Processing addbandarg at index %d", i);
5413 
5414  /* each element is a tuple */
5415  tup = (HeapTupleHeader) DatumGetPointer(e[i]);
5416  if (NULL == tup) {
5417  pfree(arg);
5418  rt_raster_destroy(raster);
5419  PG_FREE_IF_COPY(pgraster, 0);
5420  elog(ERROR, "RASTER_addBand: Invalid argument for addbandargset");
5421  PG_RETURN_NULL();
5422  }
5423 
5424  /* new band index, 1-based */
5425  arg[i].index = 0;
5426  arg[i].append = TRUE;
5427  tupv = GetAttributeByName(tup, "index", &isnull);
5428  if (!isnull) {
5429  arg[i].index = DatumGetInt32(tupv);
5430  arg[i].append = FALSE;
5431  }
5432 
5433  /* for now, only check that band index is 1-based */
5434  if (!arg[i].append && arg[i].index < 1) {
5435  pfree(arg);
5436  rt_raster_destroy(raster);
5437  PG_FREE_IF_COPY(pgraster, 0);
5438  elog(ERROR, "RASTER_addBand: Invalid argument for addbandargset. Invalid band index (must be 1-based) for addbandarg of index %d", i);
5439  PG_RETURN_NULL();
5440  }
5441 
5442  /* new band pixeltype */
5443  arg[i].pixtype = PT_END;
5444  tupv = GetAttributeByName(tup, "pixeltype", &isnull);
5445  if (isnull) {
5446  pfree(arg);
5447  rt_raster_destroy(raster);
5448  PG_FREE_IF_COPY(pgraster, 0);
5449  elog(ERROR, "RASTER_addBand: Invalid argument for addbandargset. Pixel type cannot be NULL for addbandarg of index %d", i);
5450  PG_RETURN_NULL();
5451  }
5452  text_pixtype = (text *) DatumGetPointer(tupv);
5453  if (text_pixtype == NULL) {
5454  pfree(arg);
5455  rt_raster_destroy(raster);
5456  PG_FREE_IF_COPY(pgraster, 0);
5457  elog(ERROR, "RASTER_addBand: Invalid argument for addbandargset. Pixel type cannot be NULL for addbandarg of index %d", i);
5458  PG_RETURN_NULL();
5459  }
5460  char_pixtype = text_to_cstring(text_pixtype);
5461 
5462  arg[i].pixtype = rt_pixtype_index_from_name(char_pixtype);
5463  pfree(char_pixtype);
5464  if (arg[i].pixtype == PT_END) {
5465  pfree(arg);
5466  rt_raster_destroy(raster);
5467  PG_FREE_IF_COPY(pgraster, 0);
5468  elog(ERROR, "RASTER_addBand: Invalid argument for addbandargset. Invalid pixel type for addbandarg of index %d", i);
5469  PG_RETURN_NULL();
5470  }
5471 
5472  /* new band initialvalue */
5473  arg[i].initialvalue = 0;
5474  tupv = GetAttributeByName(tup, "initialvalue", &isnull);
5475  if (!isnull)
5476  arg[i].initialvalue = DatumGetFloat8(tupv);
5477 
5478  /* new band NODATA value */
5479  arg[i].hasnodata = FALSE;
5480  arg[i].nodatavalue = 0;
5481  tupv = GetAttributeByName(tup, "nodataval", &isnull);
5482  if (!isnull) {
5483  arg[i].hasnodata = TRUE;
5484  arg[i].nodatavalue = DatumGetFloat8(tupv);
5485  }
5486  }
5487 
5488  /* add new bands to raster */
5489  lastnumbands = rt_raster_get_num_bands(raster);
5490  for (i = 0; i < n; i++) {
5491  if (nulls[i]) continue;
5492 
5493  POSTGIS_RT_DEBUGF(3, "%d bands in old raster", lastnumbands);
5494  maxbandindex = lastnumbands + 1;
5495 
5496  /* check that new band's index doesn't exceed maxbandindex */
5497  if (!arg[i].append) {
5498  if (arg[i].index > maxbandindex) {
5499  elog(NOTICE, "Band index for addbandarg of index %d exceeds possible value. Adding band at index %d", i, maxbandindex);
5500  arg[i].index = maxbandindex;
5501  }
5502  }
5503  /* append, so use maxbandindex */
5504  else
5505  arg[i].index = maxbandindex;
5506 
5507  POSTGIS_RT_DEBUGF(4, "new band (index, pixtype, initialvalue, hasnodata, nodatavalue) = (%d, %s, %f, %s, %f)",
5508  arg[i].index,
5509  rt_pixtype_name(arg[i].pixtype),
5510  arg[i].initialvalue,
5511  arg[i].hasnodata ? "TRUE" : "FALSE",
5512  arg[i].nodatavalue
5513  );
5514 
5515  bandindex = rt_raster_generate_new_band(
5516  raster,
5517  arg[i].pixtype, arg[i].initialvalue,
5518  arg[i].hasnodata, arg[i].nodatavalue,
5519  arg[i].index - 1
5520  );
5521 
5522  numbands = rt_raster_get_num_bands(raster);
5523  if (numbands == lastnumbands || bandindex == -1) {
5524  pfree(arg);
5525  rt_raster_destroy(raster);
5526  PG_FREE_IF_COPY(pgraster, 0);
5527  elog(ERROR, "RASTER_addBand: Could not add band defined by addbandarg of index %d to raster", i);
5528  PG_RETURN_NULL();
5529  }
5530 
5531  lastnumbands = numbands;
5532  POSTGIS_RT_DEBUGF(3, "%d bands in new raster", lastnumbands);
5533  }
5534 
5535  pfree(arg);
5536 
5537  pgrtn = rt_raster_serialize(raster);
5538  rt_raster_destroy(raster);
5539  PG_FREE_IF_COPY(pgraster, 0);
5540  if (!pgrtn)
5541  PG_RETURN_NULL();
5542 
5543  SET_VARSIZE(pgrtn, pgrtn->size);
5544  PG_RETURN_POINTER(pgrtn);
5545 }
5546 
5551 Datum RASTER_addBandRasterArray(PG_FUNCTION_ARGS)
5552 {
5553  rt_pgraster *pgraster = NULL;
5554  rt_pgraster *pgsrc = NULL;
5555  rt_pgraster *pgrtn = NULL;
5556 
5557  rt_raster raster = NULL;
5558  rt_raster src = NULL;
5559 
5560  int srcnband = 1;
5561  bool appendband = FALSE;
5562  int dstnband = 1;
5563  int srcnumbands = 0;
5564  int dstnumbands = 0;
5565 
5566  ArrayType *array;
5567  Oid etype;
5568  Datum *e;
5569  bool *nulls;
5570  int16 typlen;
5571  bool typbyval;
5572  char typalign;
5573  int n = 0;
5574 
5575  int rtn = 0;
5576  int i = 0;
5577 
5578  /* destination raster */
5579  if (!PG_ARGISNULL(0)) {
5580  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5581 
5582  /* raster */
5583  raster = rt_raster_deserialize(pgraster, FALSE);
5584  if (!raster) {
5585  PG_FREE_IF_COPY(pgraster, 0);
5586  elog(ERROR, "RASTER_addBandRasterArray: Could not deserialize destination raster");
5587  PG_RETURN_NULL();
5588  }
5589 
5590  POSTGIS_RT_DEBUG(4, "destination raster isn't NULL");
5591  }
5592 
5593  /* source rasters' band index, 1-based */
5594  if (!PG_ARGISNULL(2))
5595  srcnband = PG_GETARG_INT32(2);
5596  if (srcnband < 1) {
5597  elog(NOTICE, "Invalid band index for source rasters (must be 1-based). Returning original raster");
5598  if (raster != NULL) {
5599  rt_raster_destroy(raster);
5600  PG_RETURN_POINTER(pgraster);
5601  }
5602  else
5603  PG_RETURN_NULL();
5604  }
5605  POSTGIS_RT_DEBUGF(4, "srcnband = %d", srcnband);
5606 
5607  /* destination raster's band index, 1-based */
5608  if (!PG_ARGISNULL(3)) {
5609  dstnband = PG_GETARG_INT32(3);
5610  appendband = FALSE;
5611 
5612  if (dstnband < 1) {
5613  elog(NOTICE, "Invalid band index for destination raster (must be 1-based). Returning original raster");
5614  if (raster != NULL) {
5615  rt_raster_destroy(raster);
5616  PG_RETURN_POINTER(pgraster);
5617  }
5618  else
5619  PG_RETURN_NULL();
5620  }
5621  }
5622  else
5623  appendband = TRUE;
5624 
5625  /* additional processing of dstnband */
5626  if (raster != NULL) {
5627  dstnumbands = rt_raster_get_num_bands(raster);
5628 
5629  if (dstnumbands < 1) {
5630  appendband = TRUE;
5631  dstnband = 1;
5632  }
5633  else if (appendband)
5634  dstnband = dstnumbands + 1;
5635  else if (dstnband > dstnumbands) {
5636  elog(NOTICE, "Band index provided for destination raster is greater than the number of bands in the raster. Bands will be appended");
5637  appendband = TRUE;
5638  dstnband = dstnumbands + 1;
5639  }
5640  }
5641  POSTGIS_RT_DEBUGF(4, "appendband = %d", appendband);
5642  POSTGIS_RT_DEBUGF(4, "dstnband = %d", dstnband);
5643 
5644  /* process set of source rasters */
5645  POSTGIS_RT_DEBUG(3, "Processing array of source rasters");
5646  array = PG_GETARG_ARRAYTYPE_P(1);
5647  etype = ARR_ELEMTYPE(array);
5648  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
5649 
5650  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
5651  &nulls, &n);
5652 
5653  /* decrement srcnband and dstnband by 1, now 0-based */
5654  srcnband--;
5655  dstnband--;
5656  POSTGIS_RT_DEBUGF(4, "0-based nband (src, dst) = (%d, %d)", srcnband, dstnband);
5657 
5658  /* time to copy bands */
5659  for (i = 0; i < n; i++) {
5660  if (nulls[i]) continue;
5661  src = NULL;
5662 
5663  pgsrc = (rt_pgraster *) PG_DETOAST_DATUM(e[i]);
5664  src = rt_raster_deserialize(pgsrc, FALSE);
5665  if (src == NULL) {
5666  pfree(nulls);
5667  pfree(e);
5668  if (raster != NULL)
5669  rt_raster_destroy(raster);
5670  if (pgraster != NULL)
5671  PG_FREE_IF_COPY(pgraster, 0);
5672  elog(ERROR, "RASTER_addBandRasterArray: Could not deserialize source raster at index %d", i + 1);
5673  PG_RETURN_NULL();
5674  }
5675 
5676  srcnumbands = rt_raster_get_num_bands(src);
5677  POSTGIS_RT_DEBUGF(4, "source raster %d has %d bands", i + 1, srcnumbands);
5678 
5679  /* band index isn't valid */
5680  if (srcnband > srcnumbands - 1) {
5681  elog(NOTICE, "Invalid band index for source raster at index %d. Returning original raster", i + 1);
5682  pfree(nulls);
5683  pfree(e);
5684  rt_raster_destroy(src);
5685  if (raster != NULL) {
5686  rt_raster_destroy(raster);
5687  PG_RETURN_POINTER(pgraster);
5688  }
5689  else
5690  PG_RETURN_NULL();
5691  }
5692 
5693  /* destination raster is empty, new raster */
5694  if (raster == NULL) {
5695  uint32_t srcnbands[1] = {srcnband};
5696 
5697  POSTGIS_RT_DEBUG(4, "empty destination raster, using rt_raster_from_band");
5698 
5699  raster = rt_raster_from_band(src, srcnbands, 1);
5700  rt_raster_destroy(src);
5701  if (raster == NULL) {
5702  pfree(nulls);
5703  pfree(e);
5704  if (pgraster != NULL)
5705  PG_FREE_IF_COPY(pgraster, 0);
5706  elog(ERROR, "RASTER_addBandRasterArray: Could not create raster from source raster at index %d", i + 1);
5707  PG_RETURN_NULL();
5708  }
5709  }
5710  /* copy band */
5711  else {
5712  rtn = rt_raster_copy_band(
5713  raster, src,
5714  srcnband, dstnband
5715  );
5716  rt_raster_destroy(src);
5717 
5718  if (rtn == -1 || rt_raster_get_num_bands(raster) == dstnumbands) {
5719  elog(NOTICE, "Could not add band from source raster at index %d to destination raster. Returning original raster", i + 1);
5720  rt_raster_destroy(raster);
5721  pfree(nulls);
5722  pfree(e);
5723  if (pgraster != NULL)
5724  PG_RETURN_POINTER(pgraster);
5725  else
5726  PG_RETURN_NULL();
5727  }
5728  }
5729 
5730  dstnband++;
5731  dstnumbands++;
5732  }
5733 
5734  if (raster != NULL) {
5735  pgrtn = rt_raster_serialize(raster);
5736  rt_raster_destroy(raster);
5737  if (pgraster != NULL)
5738  PG_FREE_IF_COPY(pgraster, 0);
5739  if (!pgrtn)
5740  PG_RETURN_NULL();
5741 
5742  SET_VARSIZE(pgrtn, pgrtn->size);
5743  PG_RETURN_POINTER(pgrtn);
5744  }
5745 
5746  PG_RETURN_NULL();
5747 }
5748 
5753 Datum RASTER_addBandOutDB(PG_FUNCTION_ARGS)
5754 {
5755  rt_pgraster *pgraster = NULL;
5756  rt_pgraster *pgrtn = NULL;
5757 
5758  rt_raster raster = NULL;
5759  rt_band band = NULL;
5760  int numbands = 0;
5761  int dstnband = 1; /* 1-based */
5762  int appendband = FALSE;
5763  char *outdbfile = NULL;
5764  int *srcnband = NULL; /* 1-based */
5765  int numsrcnband = 0;
5766  int allbands = FALSE;
5767  int hasnodata = FALSE;
5768  double nodataval = 0.;
5769  uint16_t width = 0;
5770  uint16_t height = 0;
5771  char *authname = NULL;
5772  char *authcode = NULL;
5773 
5774  int i = 0;
5775  int j = 0;
5776 
5777  GDALDatasetH hdsOut;
5778  GDALRasterBandH hbandOut;
5779  GDALDataType gdpixtype;
5780 
5781  rt_pixtype pt = PT_END;
5782  double gt[6] = {0.};
5783  double ogt[6] = {0.};
5784  rt_raster _rast = NULL;
5785  int aligned = 0;
5786  int err = 0;
5787 
5788  /* destination raster */
5789  if (!PG_ARGISNULL(0)) {
5790  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5791 
5792  /* raster */
5793  raster = rt_raster_deserialize(pgraster, FALSE);
5794  if (!raster) {
5795  PG_FREE_IF_COPY(pgraster, 0);
5796  elog(ERROR, "RASTER_addBandOutDB: Could not deserialize destination raster");
5797  PG_RETURN_NULL();
5798  }
5799 
5800  POSTGIS_RT_DEBUG(4, "destination raster isn't NULL");
5801  }
5802 
5803  /* destination band index (1) */
5804  if (!PG_ARGISNULL(1))
5805  dstnband = PG_GETARG_INT32(1);
5806  else
5807  appendband = TRUE;
5808 
5809  /* outdb file (2) */
5810  if (PG_ARGISNULL(2)) {
5811  elog(NOTICE, "Out-db raster file not provided. Returning original raster");
5812  if (pgraster != NULL) {
5813  rt_raster_destroy(raster);
5814  PG_RETURN_POINTER(pgraster);
5815  }
5816  else
5817  PG_RETURN_NULL();
5818  }
5819  else {
5820  outdbfile = text_to_cstring(PG_GETARG_TEXT_P(2));
5821  if (!strlen(outdbfile)) {
5822  elog(NOTICE, "Out-db raster file not provided. Returning original raster");
5823  if (pgraster != NULL) {
5824  rt_raster_destroy(raster);
5825  PG_RETURN_POINTER(pgraster);
5826  }
5827  else
5828  PG_RETURN_NULL();
5829  }
5830  }
5831 
5832  /* outdb band index (3) */
5833  if (!PG_ARGISNULL(3)) {
5834  ArrayType *array;
5835  Oid etype;
5836  Datum *e;
5837  bool *nulls;
5838 
5839  int16 typlen;
5840  bool typbyval;
5841  char typalign;
5842 
5843  allbands = FALSE;
5844 
5845  array = PG_GETARG_ARRAYTYPE_P(3);
5846  etype = ARR_ELEMTYPE(array);
5847  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
5848 
5849  switch (etype) {
5850  case INT2OID:
5851  case INT4OID:
5852  break;
5853  default:
5854  if (pgraster != NULL) {
5855  rt_raster_destroy(raster);
5856  PG_FREE_IF_COPY(pgraster, 0);
5857  }
5858  elog(ERROR, "RASTER_addBandOutDB: Invalid data type for band indexes");
5859  PG_RETURN_NULL();
5860  break;
5861  }
5862 
5863  deconstruct_array(array, etype, typlen, typbyval, typalign, &e, &nulls, &numsrcnband);
5864 
5865  srcnband = palloc(sizeof(int) * numsrcnband);
5866  if (srcnband == NULL) {
5867  if (pgraster != NULL) {
5868  rt_raster_destroy(raster);
5869  PG_FREE_IF_COPY(pgraster, 0);
5870  }
5871  elog(ERROR, "RASTER_addBandOutDB: Could not allocate memory for band indexes");
5872  PG_RETURN_NULL();
5873  }
5874 
5875  for (i = 0, j = 0; i < numsrcnband; i++) {
5876  if (nulls[i]) continue;
5877 
5878  switch (etype) {
5879  case INT2OID:
5880  srcnband[j] = DatumGetInt16(e[i]);
5881  break;
5882  case INT4OID:
5883  srcnband[j] = DatumGetInt32(e[i]);
5884  break;
5885  }
5886  j++;
5887  }
5888 
5889  if (j < numsrcnband) {
5890  srcnband = repalloc(srcnband, sizeof(int) * j);
5891  if (srcnband == NULL) {
5892  if (pgraster != NULL) {
5893  rt_raster_destroy(raster);
5894  PG_FREE_IF_COPY(pgraster, 0);
5895  }
5896  elog(ERROR, "RASTER_addBandOutDB: Could not reallocate memory for band indexes");
5897  PG_RETURN_NULL();
5898  }
5899 
5900  numsrcnband = j;
5901  }
5902  }
5903  else
5904  allbands = TRUE;
5905 
5906  /* nodataval (4) */
5907  if (!PG_ARGISNULL(4)) {
5908  hasnodata = TRUE;
5909  nodataval = PG_GETARG_FLOAT8(4);
5910  }
5911  else
5912  hasnodata = FALSE;
5913 
5914  /* validate input */
5915 
5916  /* make sure dstnband is valid */
5917  if (raster != NULL) {
5918  numbands = rt_raster_get_num_bands(raster);
5919  if (!appendband) {
5920  if (dstnband < 1) {
5921  elog(NOTICE, "Invalid band index %d for adding bands. Using band index 1", dstnband);
5922  dstnband = 1;
5923  }
5924  else if (dstnband > numbands) {
5925  elog(NOTICE, "Invalid band index %d for adding bands. Using band index %d", dstnband, numbands);
5926  dstnband = numbands + 1;
5927  }
5928  }
5929  else
5930  dstnband = numbands + 1;
5931  }
5932 
5933  /* open outdb raster file */
5935  hdsOut = rt_util_gdal_open(outdbfile, GA_ReadOnly, 1);
5936  if (hdsOut == NULL) {
5937  if (pgraster != NULL) {
5938  rt_raster_destroy(raster);
5939  PG_FREE_IF_COPY(pgraster, 0);
5940  }
5941  elog(ERROR, "RASTER_addBandOutDB: Could not open out-db file with GDAL");
5942  PG_RETURN_NULL();
5943  }
5944 
5945  /* get offline raster's geotransform */
5946  if (GDALGetGeoTransform(hdsOut, ogt) != CE_None) {
5947  ogt[0] = 0;
5948  ogt[1] = 1;
5949  ogt[2] = 0;
5950  ogt[3] = 0;
5951  ogt[4] = 0;
5952  ogt[5] = -1;
5953  }
5954 
5955  /* raster doesn't exist, create it now */
5956  if (raster == NULL) {
5957  raster = rt_raster_new(GDALGetRasterXSize(hdsOut), GDALGetRasterYSize(hdsOut));
5958  if (rt_raster_is_empty(raster)) {
5959  elog(ERROR, "RASTER_addBandOutDB: Could not create new raster");
5960  PG_RETURN_NULL();
5961  }
5962  rt_raster_set_geotransform_matrix(raster, ogt);
5964 
5965  if (rt_util_gdal_sr_auth_info(hdsOut, &authname, &authcode) == ES_NONE) {
5966  if (
5967  authname != NULL &&
5968  strcmp(authname, "EPSG") == 0 &&
5969  authcode != NULL
5970  ) {
5971  rt_raster_set_srid(raster, atoi(authcode));
5972  }
5973  else
5974  elog(INFO, "Unknown SRS auth name and code from out-db file. Defaulting SRID of new raster to %d", SRID_UNKNOWN);
5975  }
5976  else
5977  elog(INFO, "Could not get SRS auth name and code from out-db file. Defaulting SRID of new raster to %d", SRID_UNKNOWN);
5978  }
5979 
5980  /* some raster info */
5981  width = rt_raster_get_width(raster);
5982  height = rt_raster_get_height(raster);
5983 
5984  /* are rasters aligned? */
5985  _rast = rt_raster_new(1, 1);
5987  rt_raster_set_srid(_rast, rt_raster_get_srid(raster));
5988  err = rt_raster_same_alignment(raster, _rast, &aligned, NULL);
5989  rt_raster_destroy(_rast);
5990 
5991  if (err != ES_NONE) {
5992  GDALClose(hdsOut);
5993  if (raster != NULL)
5994  rt_raster_destroy(raster);
5995  if (pgraster != NULL)
5996  PG_FREE_IF_COPY(pgraster, 0);
5997  elog(ERROR, "RASTER_addBandOutDB: Could not test alignment of out-db file");
5998  return ES_ERROR;
5999  }
6000  else if (!aligned)
6001  elog(WARNING, "The in-db representation of the out-db raster is not aligned. Band data may be incorrect");
6002 
6003  numbands = GDALGetRasterCount(hdsOut);
6004 
6005  /* build up srcnband */
6006  if (allbands) {
6007  numsrcnband = numbands;
6008  srcnband = palloc(sizeof(int) * numsrcnband);
6009  if (srcnband == NULL) {
6010  GDALClose(hdsOut);
6011  if (raster != NULL)
6012  rt_raster_destroy(raster);
6013  if (pgraster != NULL)
6014  PG_FREE_IF_COPY(pgraster, 0);
6015  elog(ERROR, "RASTER_addBandOutDB: Could not allocate memory for band indexes");
6016  PG_RETURN_NULL();
6017  }
6018 
6019  for (i = 0, j = 1; i < numsrcnband; i++, j++)
6020  srcnband[i] = j;
6021  }
6022 
6023  /* check band properties and add band */
6024  for (i = 0, j = dstnband - 1; i < numsrcnband; i++, j++) {
6025  /* valid index? */
6026  if (srcnband[i] < 1 || srcnband[i] > numbands) {
6027  elog(NOTICE, "Out-db file does not have a band at index %d. Returning original raster", srcnband[i]);
6028  GDALClose(hdsOut);
6029  if (raster != NULL)
6030  rt_raster_destroy(raster);
6031  if (pgraster != NULL)
6032  PG_RETURN_POINTER(pgraster);
6033  else
6034  PG_RETURN_NULL();
6035  }
6036 
6037  /* get outdb band */
6038  hbandOut = NULL;
6039  hbandOut = GDALGetRasterBand(hdsOut, srcnband[i]);
6040  if (NULL == hbandOut) {
6041  GDALClose(hdsOut);
6042  if (raster != NULL)
6043  rt_raster_destroy(raster);
6044  if (pgraster != NULL)
6045  PG_FREE_IF_COPY(pgraster, 0);
6046  elog(ERROR, "RASTER_addBandOutDB: Could not get band %d from GDAL dataset", srcnband[i]);
6047  PG_RETURN_NULL();
6048  }
6049 
6050  /* supported pixel type */
6051  gdpixtype = GDALGetRasterDataType(hbandOut);
6052  pt = rt_util_gdal_datatype_to_pixtype(gdpixtype);
6053  if (pt == PT_END) {
6054  elog(NOTICE, "Pixel type %s of band %d from GDAL dataset is not supported. Returning original raster", GDALGetDataTypeName(gdpixtype), srcnband[i]);
6055  GDALClose(hdsOut);
6056  if (raster != NULL)
6057  rt_raster_destroy(raster);
6058  if (pgraster != NULL)
6059  PG_RETURN_POINTER(pgraster);
6060  else
6061  PG_RETURN_NULL();
6062  }
6063 
6064  /* use out-db band's nodata value if nodataval not already set */
6065  if (!hasnodata)
6066  nodataval = GDALGetRasterNoDataValue(hbandOut, &hasnodata);
6067 
6068  /* add band */
6069  band = rt_band_new_offline(
6070  width, height,
6071  pt,
6072  hasnodata, nodataval,
6073  srcnband[i] - 1, outdbfile
6074  );
6075  if (band == NULL) {
6076  GDALClose(hdsOut);
6077  if (raster != NULL)
6078  rt_raster_destroy(raster);
6079  if (pgraster != NULL)
6080  PG_FREE_IF_COPY(pgraster, 0);
6081  elog(ERROR, "RASTER_addBandOutDB: Could not create new out-db band");
6082  PG_RETURN_NULL();
6083  }
6084 
6085  if (rt_raster_add_band(raster, band, j) < 0) {
6086  GDALClose(hdsOut);
6087  if (raster != NULL)
6088  rt_raster_destroy(raster);
6089  if (pgraster != NULL)
6090  PG_FREE_IF_COPY(pgraster, 0);
6091  elog(ERROR, "RASTER_addBandOutDB: Could not add new out-db band to raster");
6092  PG_RETURN_NULL();
6093  }
6094  }
6095 
6096  pgrtn = rt_raster_serialize(raster);
6097  rt_raster_destroy(raster);
6098  if (pgraster != NULL)
6099  PG_FREE_IF_COPY(pgraster, 0);
6100  if (!pgrtn)
6101  PG_RETURN_NULL();
6102 
6103  SET_VARSIZE(pgrtn, pgrtn->size);
6104  PG_RETURN_POINTER(pgrtn);
6105 }
6106 
6111 Datum RASTER_tile(PG_FUNCTION_ARGS)
6112 {
6113  FuncCallContext *funcctx;
6114  int call_cntr;
6115  int max_calls;
6116  int i = 0;
6117  int j = 0;
6118 
6119  struct tile_arg_t {
6120 
6121  struct {
6122  rt_raster raster;
6123  double gt[6];
6124  int srid;
6125  int width;
6126  int height;
6127  } raster;
6128 
6129  struct {
6130  int width;
6131  int height;
6132 
6133  int nx;
6134  int ny;
6135  } tile;
6136 
6137  int numbands;
6138  int *nbands;
6139 
6140  struct {
6141  int pad;
6142  double hasnodata;
6143  double nodataval;
6144  } pad;
6145  };
6146  struct tile_arg_t *arg1 = NULL;
6147  struct tile_arg_t *arg2 = NULL;
6148 
6149  if (SRF_IS_FIRSTCALL()) {
6150  MemoryContext oldcontext;
6151  rt_pgraster *pgraster = NULL;
6152  int numbands;
6153 
6154  ArrayType *array;
6155  Oid etype;
6156  Datum *e;
6157  bool *nulls;
6158 
6159  int16 typlen;
6160  bool typbyval;
6161  char typalign;
6162 
6163  POSTGIS_RT_DEBUG(2, "RASTER_tile: first call");
6164 
6165  /* create a function context for cross-call persistence */
6166  funcctx = SRF_FIRSTCALL_INIT();
6167 
6168  /* switch to memory context appropriate for multiple function calls */
6169  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
6170 
6171  /* Get input arguments */
6172  if (PG_ARGISNULL(0)) {
6173  MemoryContextSwitchTo(oldcontext);
6174  SRF_RETURN_DONE(funcctx);
6175  }
6176 
6177  /* allocate arg1 */
6178  arg1 = palloc(sizeof(struct tile_arg_t));
6179  if (arg1 == NULL) {
6180  MemoryContextSwitchTo(oldcontext);
6181  elog(ERROR, "RASTER_tile: Could not allocate memory for arguments");
6182  SRF_RETURN_DONE(funcctx);
6183  }
6184 
6185  pgraster = (rt_pgraster *) PG_DETOAST_DATUM_COPY(PG_GETARG_DATUM(0));
6186  arg1->raster.raster = rt_raster_deserialize(pgraster, FALSE);
6187  if (!arg1->raster.raster) {
6188  ereport(ERROR, (
6189  errcode(ERRCODE_OUT_OF_MEMORY),
6190  errmsg("Could not deserialize raster")
6191  ));
6192  pfree(arg1);
6193  PG_FREE_IF_COPY(pgraster, 0);
6194  MemoryContextSwitchTo(oldcontext);
6195  SRF_RETURN_DONE(funcctx);
6196  }
6197 
6198  /* raster has bands */
6199  numbands = rt_raster_get_num_bands(arg1->raster.raster);
6200  /*
6201  if (!numbands) {
6202  elog(NOTICE, "Raster provided has no bands");
6203  rt_raster_destroy(arg1->raster.raster);
6204  pfree(arg1);
6205  PG_FREE_IF_COPY(pgraster, 0);
6206  MemoryContextSwitchTo(oldcontext);
6207  SRF_RETURN_DONE(funcctx);
6208  }
6209  */
6210 
6211  /* width (1) */
6212  if (PG_ARGISNULL(1)) {
6213  elog(NOTICE, "Width cannot be NULL. Returning NULL");
6214  rt_raster_destroy(arg1->raster.raster);
6215  pfree(arg1);
6216  PG_FREE_IF_COPY(pgraster, 0);
6217  MemoryContextSwitchTo(oldcontext);
6218  SRF_RETURN_DONE(funcctx);
6219  }
6220  arg1->tile.width = PG_GETARG_INT32(1);
6221  if (arg1->tile.width < 1) {
6222  elog(NOTICE, "Width must be greater than zero. Returning NULL");
6223  rt_raster_destroy(arg1->raster.raster);
6224  pfree(arg1);
6225  PG_FREE_IF_COPY(pgraster, 0);
6226  MemoryContextSwitchTo(oldcontext);
6227  SRF_RETURN_DONE(funcctx);
6228  }
6229 
6230  /* height (2) */
6231  if (PG_ARGISNULL(2)) {
6232  elog(NOTICE, "Height cannot be NULL. Returning NULL");
6233  rt_raster_destroy(arg1->raster.raster);
6234  pfree(arg1);
6235  PG_FREE_IF_COPY(pgraster, 0);
6236  MemoryContextSwitchTo(oldcontext);
6237  SRF_RETURN_DONE(funcctx);
6238  }
6239  arg1->tile.height = PG_GETARG_INT32(2);
6240  if (arg1->tile.height < 1) {
6241  elog(NOTICE, "Height must be greater than zero. Returning NULL");
6242  rt_raster_destroy(arg1->raster.raster);
6243  pfree(arg1);
6244  PG_FREE_IF_COPY(pgraster, 0);
6245  MemoryContextSwitchTo(oldcontext);
6246  SRF_RETURN_DONE(funcctx);
6247  }
6248 
6249  /* nband, array (3) */
6250  if (numbands && !PG_ARGISNULL(3)) {
6251  array = PG_GETARG_ARRAYTYPE_P(3);
6252  etype = ARR_ELEMTYPE(array);
6253  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
6254 
6255  switch (etype) {
6256  case INT2OID:
6257  case INT4OID:
6258  break;
6259  default:
6260  rt_raster_destroy(arg1->raster.raster);
6261  pfree(arg1);
6262  PG_FREE_IF_COPY(pgraster, 0);
6263  MemoryContextSwitchTo(oldcontext);
6264  elog(ERROR, "RASTER_tile: Invalid data type for band indexes");
6265  SRF_RETURN_DON