PostGIS  3.4.0dev-r@@SVN_REVISION@@
rtpg_statistics.c
Go to the documentation of this file.
1 /*
2  *
3  * WKTRaster - Raster Types for PostGIS
4  * http://trac.osgeo.org/postgis/wiki/WKTRaster
5  *
6  * Copyright (C) 2013 Bborie Park <dustymugs@gmail.com>
7  * Copyright (C) 2011-2013 Regents of the University of California
8  * <bkpark@ucdavis.edu>
9  * Copyright (C) 2010-2011 Jorge Arevalo <jorge.arevalo@deimos-space.com>
10  * Copyright (C) 2010-2011 David Zwarg <dzwarg@azavea.com>
11  * Copyright (C) 2009-2011 Pierre Racine <pierre.racine@sbf.ulaval.ca>
12  * Copyright (C) 2009-2011 Mateusz Loskot <mateusz@loskot.net>
13  * Copyright (C) 2008-2009 Sandro Santilli <strk@kbt.io>
14  *
15  * This program is free software; you can redistribute it and/or
16  * modify it under the terms of the GNU General Public License
17  * as published by the Free Software Foundation; either version 2
18  * of the License, or (at your option) any later version.
19  *
20  * This program is distributed in the hope that it will be useful,
21  * but WITHOUT ANY WARRANTY; without even the implied warranty of
22  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23  * GNU General Public License for more details.
24  *
25  * You should have received a copy of the GNU General Public License
26  * along with this program; if not, write to the Free Software Foundation,
27  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
28  *
29  */
30 
31 #include <postgres.h>
32 #include <fmgr.h>
33 #include <utils/builtins.h> /* for text_to_cstring() */
34 #include "utils/lsyscache.h" /* for get_typlenbyvalalign */
35 #include "utils/array.h" /* for ArrayType */
36 #include "catalog/pg_type.h" /* for INT2OID, INT4OID, FLOAT4OID, FLOAT8OID and TEXTOID */
37 #include <executor/spi.h>
38 #include <funcapi.h> /* for SRF */
39 
40 #include "../../postgis_config.h"
41 
42 
43 #include "access/htup_details.h" /* for heap_form_tuple() */
44 
45 
46 #include "rtpostgis.h"
47 
48 /* Get summary stats */
49 Datum RASTER_summaryStats(PG_FUNCTION_ARGS);
50 Datum RASTER_summaryStatsCoverage(PG_FUNCTION_ARGS);
51 
52 Datum RASTER_summaryStats_transfn(PG_FUNCTION_ARGS);
53 Datum RASTER_summaryStats_finalfn(PG_FUNCTION_ARGS);
54 
55 /* get histogram */
56 Datum RASTER_histogram(PG_FUNCTION_ARGS);
57 
58 /* get quantiles */
59 Datum RASTER_quantile(PG_FUNCTION_ARGS);
60 
61 /* get counts of values */
62 Datum RASTER_valueCount(PG_FUNCTION_ARGS);
63 Datum RASTER_valueCountCoverage(PG_FUNCTION_ARGS);
64 
65 #define VALUES_LENGTH 6
66 
70 PG_FUNCTION_INFO_V1(RASTER_summaryStats);
71 Datum RASTER_summaryStats(PG_FUNCTION_ARGS)
72 {
73  rt_pgraster *pgraster = NULL;
74  rt_raster raster = NULL;
75  rt_band band = NULL;
76  int32_t bandindex = 1;
77  bool exclude_nodata_value = TRUE;
78  int num_bands = 0;
79  double sample = 0;
80  rt_bandstats stats = NULL;
81 
82  TupleDesc tupdesc;
83  Datum values[VALUES_LENGTH];
84  bool nulls[VALUES_LENGTH];
85  HeapTuple tuple;
86  Datum result;
87 
88  /* pgraster is null, return null */
89  if (PG_ARGISNULL(0))
90  PG_RETURN_NULL();
91  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
92 
93  raster = rt_raster_deserialize(pgraster, FALSE);
94  if (!raster) {
95  PG_FREE_IF_COPY(pgraster, 0);
96  elog(ERROR, "RASTER_summaryStats: Cannot deserialize raster");
97  PG_RETURN_NULL();
98  }
99 
100  /* band index is 1-based */
101  if (!PG_ARGISNULL(1))
102  bandindex = PG_GETARG_INT32(1);
103  num_bands = rt_raster_get_num_bands(raster);
104  if (bandindex < 1 || bandindex > num_bands) {
105  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
106  rt_raster_destroy(raster);
107  PG_FREE_IF_COPY(pgraster, 0);
108  PG_RETURN_NULL();
109  }
110 
111  /* exclude_nodata_value flag */
112  if (!PG_ARGISNULL(2))
113  exclude_nodata_value = PG_GETARG_BOOL(2);
114 
115  /* sample % */
116  if (!PG_ARGISNULL(3)) {
117  sample = PG_GETARG_FLOAT8(3);
118  if (sample < 0 || sample > 1) {
119  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
120  rt_raster_destroy(raster);
121  PG_FREE_IF_COPY(pgraster, 0);
122  PG_RETURN_NULL();
123  }
124  else if (FLT_EQ(sample, 0.0))
125  sample = 1;
126  }
127  else
128  sample = 1;
129 
130  /* get band */
131  band = rt_raster_get_band(raster, bandindex - 1);
132  if (!band) {
133  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
134  rt_raster_destroy(raster);
135  PG_FREE_IF_COPY(pgraster, 0);
136  PG_RETURN_NULL();
137  }
138 
139  /* we don't need the raw values, hence the zero parameter */
140  stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 0, NULL, NULL, NULL);
141  rt_band_destroy(band);
142  rt_raster_destroy(raster);
143  PG_FREE_IF_COPY(pgraster, 0);
144  if (NULL == stats) {
145  elog(NOTICE, "Cannot compute summary statistics for band at index %d. Returning NULL", bandindex);
146  PG_RETURN_NULL();
147  }
148 
149  /* Build a tuple descriptor for our result type */
150  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
151  ereport(ERROR, (
152  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
153  errmsg(
154  "function returning record called in context "
155  "that cannot accept type record"
156  )
157  ));
158  }
159 
160  BlessTupleDesc(tupdesc);
161 
162  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
163 
164  values[0] = Int64GetDatum(stats->count);
165  if (stats->count > 0) {
166  values[1] = Float8GetDatum(stats->sum);
167  values[2] = Float8GetDatum(stats->mean);
168  values[3] = Float8GetDatum(stats->stddev);
169  values[4] = Float8GetDatum(stats->min);
170  values[5] = Float8GetDatum(stats->max);
171  }
172  else {
173  nulls[1] = TRUE;
174  nulls[2] = TRUE;
175  nulls[3] = TRUE;
176  nulls[4] = TRUE;
177  nulls[5] = TRUE;
178  }
179 
180  /* build a tuple */
181  tuple = heap_form_tuple(tupdesc, values, nulls);
182 
183  /* make the tuple into a datum */
184  result = HeapTupleGetDatum(tuple);
185 
186  /* clean up */
187  pfree(stats);
188 
189  PG_RETURN_DATUM(result);
190 }
191 
195 PG_FUNCTION_INFO_V1(RASTER_summaryStatsCoverage);
196 Datum RASTER_summaryStatsCoverage(PG_FUNCTION_ARGS)
197 {
198  text *tablenametext = NULL;
199  char *tablename = NULL;
200  text *colnametext = NULL;
201  char *colname = NULL;
202  int32_t bandindex = 1;
203  bool exclude_nodata_value = TRUE;
204  double sample = 0;
205 
206  int len = 0;
207  char *sql = NULL;
208  int spi_result;
209  Portal portal;
210  TupleDesc tupdesc;
211  SPITupleTable *tuptable = NULL;
212  HeapTuple tuple;
213  Datum datum;
214  bool isNull = FALSE;
215 
216  rt_pgraster *pgraster = NULL;
217  rt_raster raster = NULL;
218  rt_band band = NULL;
219  int num_bands = 0;
220  uint64_t cK = 0;
221  double cM = 0;
222  double cQ = 0;
223  rt_bandstats stats = NULL;
224  rt_bandstats rtn = NULL;
225 
226  Datum values[VALUES_LENGTH];
227  bool nulls[VALUES_LENGTH];
228  Datum result;
229 
230  /* tablename is null, return null */
231  if (PG_ARGISNULL(0)) {
232  elog(NOTICE, "Table name must be provided");
233  PG_RETURN_NULL();
234  }
235  tablenametext = PG_GETARG_TEXT_P(0);
236  tablename = text_to_cstring(tablenametext);
237  if (!strlen(tablename)) {
238  elog(NOTICE, "Table name must be provided");
239  PG_RETURN_NULL();
240  }
241 
242  /* column name is null, return null */
243  if (PG_ARGISNULL(1)) {
244  elog(NOTICE, "Column name must be provided");
245  PG_RETURN_NULL();
246  }
247  colnametext = PG_GETARG_TEXT_P(1);
248  colname = text_to_cstring(colnametext);
249  if (!strlen(colname)) {
250  elog(NOTICE, "Column name must be provided");
251  PG_RETURN_NULL();
252  }
253 
254  /* band index is 1-based */
255  if (!PG_ARGISNULL(2))
256  bandindex = PG_GETARG_INT32(2);
257 
258  /* exclude_nodata_value flag */
259  if (!PG_ARGISNULL(3))
260  exclude_nodata_value = PG_GETARG_BOOL(3);
261 
262  /* sample % */
263  if (!PG_ARGISNULL(4)) {
264  sample = PG_GETARG_FLOAT8(4);
265  if (sample < 0 || sample > 1) {
266  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
267  rt_raster_destroy(raster);
268  PG_RETURN_NULL();
269  }
270  else if (FLT_EQ(sample, 0.0))
271  sample = 1;
272  }
273  else
274  sample = 1;
275 
276  /* iterate through rasters of coverage */
277  /* connect to database */
278  spi_result = SPI_connect();
279  if (spi_result != SPI_OK_CONNECT) {
280  pfree(sql);
281  elog(ERROR, "RASTER_summaryStatsCoverage: Cannot connect to database using SPI");
282  PG_RETURN_NULL();
283  }
284 
285  /* create sql */
286  len = sizeof(char) * (strlen("SELECT \"\" FROM \"\" WHERE \"\" IS NOT NULL") + (strlen(colname) * 2) + strlen(tablename) + 1);
287  sql = (char *) palloc(len);
288  if (NULL == sql) {
289  if (SPI_tuptable) SPI_freetuptable(tuptable);
290  SPI_finish();
291  elog(ERROR, "RASTER_summaryStatsCoverage: Cannot allocate memory for sql");
292  PG_RETURN_NULL();
293  }
294 
295  /* get cursor */
296  snprintf(sql, len, "SELECT \"%s\" FROM \"%s\" WHERE \"%s\" IS NOT NULL", colname, tablename, colname);
297  portal = SPI_cursor_open_with_args(
298  "coverage",
299  sql,
300  0, NULL,
301  NULL, NULL,
302  TRUE, 0
303  );
304  pfree(sql);
305 
306  /* process resultset */
307  SPI_cursor_fetch(portal, TRUE, 1);
308  while (SPI_processed == 1 && SPI_tuptable != NULL) {
309  tupdesc = SPI_tuptable->tupdesc;
310  tuple = SPI_tuptable->vals[0];
311 
312  datum = SPI_getbinval(tuple, tupdesc, 1, &isNull);
313  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
314  SPI_freetuptable(SPI_tuptable);
315  SPI_cursor_close(portal);
316  SPI_finish();
317 
318  if (NULL != rtn) pfree(rtn);
319  elog(ERROR, "RASTER_summaryStatsCoverage: Cannot get raster of coverage");
320  PG_RETURN_NULL();
321  }
322  else if (isNull) {
323  SPI_cursor_fetch(portal, TRUE, 1);
324  continue;
325  }
326 
327  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(datum);
328 
329  raster = rt_raster_deserialize(pgraster, FALSE);
330  if (!raster) {
331  SPI_freetuptable(SPI_tuptable);
332  SPI_cursor_close(portal);
333  SPI_finish();
334 
335  if (NULL != rtn) pfree(rtn);
336  elog(ERROR, "RASTER_summaryStatsCoverage: Cannot deserialize raster");
337  PG_RETURN_NULL();
338  }
339 
340  /* inspect number of bands */
341  num_bands = rt_raster_get_num_bands(raster);
342  if (bandindex < 1 || bandindex > num_bands) {
343  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
344 
345  rt_raster_destroy(raster);
346 
347  SPI_freetuptable(SPI_tuptable);
348  SPI_cursor_close(portal);
349  SPI_finish();
350 
351  if (NULL != rtn) pfree(rtn);
352  PG_RETURN_NULL();
353  }
354 
355  /* get band */
356  band = rt_raster_get_band(raster, bandindex - 1);
357  if (!band) {
358  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
359 
360  rt_raster_destroy(raster);
361 
362  SPI_freetuptable(SPI_tuptable);
363  SPI_cursor_close(portal);
364  SPI_finish();
365 
366  if (NULL != rtn) pfree(rtn);
367  PG_RETURN_NULL();
368  }
369 
370  /* we don't need the raw values, hence the zero parameter */
371  stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 0, &cK, &cM, &cQ);
372 
373  rt_band_destroy(band);
374  rt_raster_destroy(raster);
375 
376  if (NULL == stats) {
377  elog(NOTICE, "Cannot compute summary statistics for band at index %d. Returning NULL", bandindex);
378 
379  SPI_freetuptable(SPI_tuptable);
380  SPI_cursor_close(portal);
381  SPI_finish();
382 
383  if (NULL != rtn) pfree(rtn);
384  PG_RETURN_NULL();
385  }
386 
387  /* initialize rtn */
388  if (stats->count > 0) {
389  if (NULL == rtn) {
390  rtn = (rt_bandstats) SPI_palloc(sizeof(struct rt_bandstats_t));
391  if (NULL == rtn) {
392  SPI_freetuptable(SPI_tuptable);
393  SPI_cursor_close(portal);
394  SPI_finish();
395 
396  elog(ERROR, "RASTER_summaryStatsCoverage: Cannot allocate memory for summary stats of coverage");
397  PG_RETURN_NULL();
398  }
399 
400  rtn->sample = stats->sample;
401  rtn->count = stats->count;
402  rtn->min = stats->min;
403  rtn->max = stats->max;
404  rtn->sum = stats->sum;
405  rtn->mean = stats->mean;
406  rtn->stddev = -1;
407 
408  rtn->values = NULL;
409  rtn->sorted = 0;
410  }
411  else {
412  rtn->count += stats->count;
413  rtn->sum += stats->sum;
414 
415  if (stats->min < rtn->min)
416  rtn->min = stats->min;
417  if (stats->max > rtn->max)
418  rtn->max = stats->max;
419  }
420  }
421 
422  pfree(stats);
423 
424  /* next record */
425  SPI_cursor_fetch(portal, TRUE, 1);
426  }
427 
428  if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
429  SPI_cursor_close(portal);
430  SPI_finish();
431 
432  if (NULL == rtn) {
433  elog(ERROR, "RASTER_summaryStatsCoverage: Cannot compute coverage summary stats");
434  PG_RETURN_NULL();
435  }
436 
437  /* coverage mean and deviation */
438  rtn->mean = rtn->sum / rtn->count;
439  /* sample deviation */
440  if (rtn->sample > 0 && rtn->sample < 1)
441  rtn->stddev = sqrt(cQ / (rtn->count - 1));
442  /* standard deviation */
443  else
444  rtn->stddev = sqrt(cQ / rtn->count);
445 
446  /* Build a tuple descriptor for our result type */
447  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
448  ereport(ERROR, (
449  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
450  errmsg(
451  "function returning record called in context "
452  "that cannot accept type record"
453  )
454  ));
455  }
456 
457  BlessTupleDesc(tupdesc);
458 
459  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
460 
461  values[0] = Int64GetDatum(rtn->count);
462  if (rtn->count > 0) {
463  values[1] = Float8GetDatum(rtn->sum);
464  values[2] = Float8GetDatum(rtn->mean);
465  values[3] = Float8GetDatum(rtn->stddev);
466  values[4] = Float8GetDatum(rtn->min);
467  values[5] = Float8GetDatum(rtn->max);
468  }
469  else {
470  nulls[1] = TRUE;
471  nulls[2] = TRUE;
472  nulls[3] = TRUE;
473  nulls[4] = TRUE;
474  nulls[5] = TRUE;
475  }
476 
477  /* build a tuple */
478  tuple = heap_form_tuple(tupdesc, values, nulls);
479 
480  /* make the tuple into a datum */
481  result = HeapTupleGetDatum(tuple);
482 
483  /* clean up */
484  pfree(rtn);
485 
486  PG_RETURN_DATUM(result);
487 }
488 
489 /* ---------------------------------------------------------------- */
490 /* Aggregate ST_SummaryStats */
491 /* ---------------------------------------------------------------- */
492 
493 typedef struct rtpg_summarystats_arg_t *rtpg_summarystats_arg;
494 struct rtpg_summarystats_arg_t {
495  rt_bandstats stats;
496 
497  /* coefficients for one-pass standard deviation */
498  uint64_t cK;
499  double cM;
500  double cQ;
501 
502  int32_t band_index; /* one-based */
503  bool exclude_nodata_value;
504  double sample; /* value between 0 and 1 */
505 };
506 
507 static void
508 rtpg_summarystats_arg_destroy(rtpg_summarystats_arg arg) {
509  if (arg->stats != NULL)
510  pfree(arg->stats);
511 
512  pfree(arg);
513 }
514 
515 static rtpg_summarystats_arg
516 rtpg_summarystats_arg_init() {
517  rtpg_summarystats_arg arg = NULL;
518 
519  arg = palloc(sizeof(struct rtpg_summarystats_arg_t));
520  if (arg == NULL) {
521  elog(
522  ERROR,
523  "rtpg_summarystats_arg_init: Cannot allocate memory for function arguments"
524  );
525  return NULL;
526  }
527 
528  arg->stats = (rt_bandstats) palloc(sizeof(struct rt_bandstats_t));
529  if (arg->stats == NULL) {
530  rtpg_summarystats_arg_destroy(arg);
531  elog(
532  ERROR,
533  "rtpg_summarystats_arg_init: Cannot allocate memory for stats function argument"
534  );
535  return NULL;
536  }
537 
538  arg->stats->sample = 0;
539  arg->stats->count = 0;
540  arg->stats->min = 0;
541  arg->stats->max = 0;
542  arg->stats->sum = 0;
543  arg->stats->mean = 0;
544  arg->stats->stddev = -1;
545  arg->stats->values = NULL;
546  arg->stats->sorted = 0;
547 
548  arg->cK = 0;
549  arg->cM = 0;
550  arg->cQ = 0;
551 
552  arg->band_index = 1;
553  arg->exclude_nodata_value = TRUE;
554  arg->sample = 1;
555 
556  return arg;
557 }
558 
559 PG_FUNCTION_INFO_V1(RASTER_summaryStats_transfn);
560 Datum RASTER_summaryStats_transfn(PG_FUNCTION_ARGS)
561 {
562  MemoryContext aggcontext;
563  MemoryContext oldcontext;
564  rtpg_summarystats_arg state = NULL;
565  bool skiparg = FALSE;
566 
567  int i = 0;
568 
569  rt_pgraster *pgraster = NULL;
570  rt_raster raster = NULL;
571  rt_band band = NULL;
572  int num_bands = 0;
573  rt_bandstats stats = NULL;
574 
575  POSTGIS_RT_DEBUG(3, "Starting...");
576 
577  /* cannot be called directly as this is exclusive aggregate function */
578  if (!AggCheckCallContext(fcinfo, &aggcontext)) {
579  elog(
580  ERROR,
581  "RASTER_summaryStats_transfn: Cannot be called in a non-aggregate context"
582  );
583  PG_RETURN_NULL();
584  }
585 
586  /* switch to aggcontext */
587  oldcontext = MemoryContextSwitchTo(aggcontext);
588 
589  if (PG_ARGISNULL(0)) {
590  POSTGIS_RT_DEBUG(3, "Creating state variable");
591 
592  state = rtpg_summarystats_arg_init();
593  if (state == NULL) {
594  MemoryContextSwitchTo(oldcontext);
595  elog(
596  ERROR,
597  "RASTER_summaryStats_transfn: Cannot allocate memory for state variable"
598  );
599  PG_RETURN_NULL();
600  }
601 
602  skiparg = FALSE;
603  }
604  else {
605  POSTGIS_RT_DEBUG(3, "State variable already exists");
606  state = (rtpg_summarystats_arg) PG_GETARG_POINTER(0);
607  skiparg = TRUE;
608  }
609 
610  /* raster arg is NOT NULL */
611  if (!PG_ARGISNULL(1)) {
612  /* deserialize raster */
613  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
614 
615  /* Get raster object */
616  raster = rt_raster_deserialize(pgraster, FALSE);
617  if (raster == NULL) {
618 
619  rtpg_summarystats_arg_destroy(state);
620  PG_FREE_IF_COPY(pgraster, 1);
621 
622  MemoryContextSwitchTo(oldcontext);
623  elog(ERROR, "RASTER_summaryStats_transfn: Cannot deserialize raster");
624  PG_RETURN_NULL();
625  }
626  }
627 
628  do {
629  Oid calltype;
630  int nargs = 0;
631 
632  if (skiparg)
633  break;
634 
635  /* 4 or 5 total possible args */
636  nargs = PG_NARGS();
637  POSTGIS_RT_DEBUGF(4, "nargs = %d", nargs);
638 
639  for (i = 2; i < nargs; i++) {
640  if (PG_ARGISNULL(i))
641  continue;
642 
643  calltype = get_fn_expr_argtype(fcinfo->flinfo, i);
644 
645  /* band index */
646  if (
647  (calltype == INT2OID || calltype == INT4OID) &&
648  i == 2
649  ) {
650  if (calltype == INT2OID)
651  state->band_index = PG_GETARG_INT16(i);
652  else
653  state->band_index = PG_GETARG_INT32(i);
654 
655  /* basic check, > 0 */
656  if (state->band_index < 1) {
657 
658  rtpg_summarystats_arg_destroy(state);
659  if (raster != NULL) {
660  rt_raster_destroy(raster);
661  PG_FREE_IF_COPY(pgraster, 1);
662  }
663 
664  MemoryContextSwitchTo(oldcontext);
665  elog(
666  ERROR,
667  "RASTER_summaryStats_transfn: Invalid band index (must use 1-based). Returning NULL"
668  );
669  PG_RETURN_NULL();
670  }
671  }
672  /* exclude_nodata_value */
673  else if (
674  calltype == BOOLOID && (
675  i == 2 || i == 3
676  )
677  ) {
678  state->exclude_nodata_value = PG_GETARG_BOOL(i);
679  }
680  /* sample rate */
681  else if (
682  (calltype == FLOAT4OID || calltype == FLOAT8OID) &&
683  (i == 3 || i == 4)
684  ) {
685  if (calltype == FLOAT4OID)
686  state->sample = PG_GETARG_FLOAT4(i);
687  else
688  state->sample = PG_GETARG_FLOAT8(i);
689 
690  /* basic check, 0 <= sample <= 1 */
691  if (state->sample < 0. || state->sample > 1.) {
692 
693  rtpg_summarystats_arg_destroy(state);
694  if (raster != NULL) {
695  rt_raster_destroy(raster);
696  PG_FREE_IF_COPY(pgraster, 1);
697  }
698 
699  MemoryContextSwitchTo(oldcontext);
700  elog(
701  ERROR,
702  "Invalid sample percentage (must be between 0 and 1). Returning NULL"
703  );
704 
705  PG_RETURN_NULL();
706  }
707  else if (FLT_EQ(state->sample, 0.0))
708  state->sample = 1;
709  }
710  /* unknown arg */
711  else {
712  rtpg_summarystats_arg_destroy(state);
713  if (raster != NULL) {
714  rt_raster_destroy(raster);
715  PG_FREE_IF_COPY(pgraster, 1);
716  }
717 
718  MemoryContextSwitchTo(oldcontext);
719  elog(
720  ERROR,
721  "RASTER_summaryStats_transfn: Unknown function parameter at index %d",
722  i
723  );
724  PG_RETURN_NULL();
725  }
726  }
727  }
728  while (0);
729 
730  /* null raster, return */
731  if (PG_ARGISNULL(1)) {
732  POSTGIS_RT_DEBUG(4, "NULL raster so processing required");
733  MemoryContextSwitchTo(oldcontext);
734  PG_RETURN_POINTER(state);
735  }
736 
737  /* inspect number of bands */
738  num_bands = rt_raster_get_num_bands(raster);
739  if (state->band_index > num_bands) {
740  elog(
741  NOTICE,
742  "Raster does not have band at index %d. Skipping raster",
743  state->band_index
744  );
745 
746  rt_raster_destroy(raster);
747  PG_FREE_IF_COPY(pgraster, 1);
748 
749  MemoryContextSwitchTo(oldcontext);
750  PG_RETURN_POINTER(state);
751  }
752 
753  /* get band */
754  band = rt_raster_get_band(raster, state->band_index - 1);
755  if (!band) {
756  elog(
757  NOTICE, "Cannot find band at index %d. Skipping raster",
758  state->band_index
759  );
760 
761  rt_raster_destroy(raster);
762  PG_FREE_IF_COPY(pgraster, 1);
763 
764  MemoryContextSwitchTo(oldcontext);
765  PG_RETURN_POINTER(state);
766  }
767 
768  /* we don't need the raw values, hence the zero parameter */
769  stats = rt_band_get_summary_stats(
770  band, (int) state->exclude_nodata_value,
771  state->sample, 0,
772  &(state->cK), &(state->cM), &(state->cQ)
773  );
774 
775  rt_band_destroy(band);
776  rt_raster_destroy(raster);
777  PG_FREE_IF_COPY(pgraster, 1);
778 
779  if (NULL == stats) {
780  elog(
781  NOTICE,
782  "Cannot compute summary statistics for band at index %d. Returning NULL",
783  state->band_index
784  );
785 
786  rtpg_summarystats_arg_destroy(state);
787 
788  MemoryContextSwitchTo(oldcontext);
789  PG_RETURN_NULL();
790  }
791 
792  if (stats->count > 0) {
793  if (state->stats->count < 1) {
794  state->stats->sample = stats->sample;
795  state->stats->count = stats->count;
796  state->stats->min = stats->min;
797  state->stats->max = stats->max;
798  state->stats->sum = stats->sum;
799  state->stats->mean = stats->mean;
800  state->stats->stddev = -1;
801  }
802  else {
803  state->stats->count += stats->count;
804  state->stats->sum += stats->sum;
805 
806  if (stats->min < state->stats->min)
807  state->stats->min = stats->min;
808  if (stats->max > state->stats->max)
809  state->stats->max = stats->max;
810  }
811  }
812 
813  pfree(stats);
814 
815  /* switch back to local context */
816  MemoryContextSwitchTo(oldcontext);
817 
818  POSTGIS_RT_DEBUG(3, "Finished");
819 
820  PG_RETURN_POINTER(state);
821 }
822 
823 PG_FUNCTION_INFO_V1(RASTER_summaryStats_finalfn);
824 Datum RASTER_summaryStats_finalfn(PG_FUNCTION_ARGS)
825 {
826  rtpg_summarystats_arg state = NULL;
827 
828  TupleDesc tupdesc;
829  HeapTuple tuple;
830  Datum values[VALUES_LENGTH];
831  bool nulls[VALUES_LENGTH];
832  Datum result;
833 
834  POSTGIS_RT_DEBUG(3, "Starting...");
835 
836  /* cannot be called directly as this is exclusive aggregate function */
837  if (!AggCheckCallContext(fcinfo, NULL)) {
838  elog(ERROR, "RASTER_summaryStats_finalfn: Cannot be called in a non-aggregate context");
839  PG_RETURN_NULL();
840  }
841 
842  /* NULL, return null */
843  if (PG_ARGISNULL(0))
844  PG_RETURN_NULL();
845 
846  state = (rtpg_summarystats_arg) PG_GETARG_POINTER(0);
847 
848  if (NULL == state) {
849  elog(ERROR, "RASTER_summaryStats_finalfn: Cannot compute coverage summary stats");
850  PG_RETURN_NULL();
851  }
852 
853  /* coverage mean and deviation */
854  if (state->stats->count > 0) {
855  state->stats->mean = state->stats->sum / state->stats->count;
856  /* sample deviation */
857  if (state->stats->sample > 0 && state->stats->sample < 1)
858  state->stats->stddev = sqrt(state->cQ / (state->stats->count - 1));
859  /* standard deviation */
860  else
861  state->stats->stddev = sqrt(state->cQ / state->stats->count);
862  }
863 
864  /* Build a tuple descriptor for our result type */
865  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
866  rtpg_summarystats_arg_destroy(state);
867  ereport(ERROR, (
868  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
869  errmsg(
870  "function returning record called in context "
871  "that cannot accept type record"
872  )
873  ));
874  }
875 
876  BlessTupleDesc(tupdesc);
877 
878  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
879 
880  values[0] = Int64GetDatum(state->stats->count);
881  if (state->stats->count > 0) {
882  values[1] = Float8GetDatum(state->stats->sum);
883  values[2] = Float8GetDatum(state->stats->mean);
884  values[3] = Float8GetDatum(state->stats->stddev);
885  values[4] = Float8GetDatum(state->stats->min);
886  values[5] = Float8GetDatum(state->stats->max);
887  }
888  else {
889  nulls[1] = TRUE;
890  nulls[2] = TRUE;
891  nulls[3] = TRUE;
892  nulls[4] = TRUE;
893  nulls[5] = TRUE;
894  }
895 
896  /* build a tuple */
897  tuple = heap_form_tuple(tupdesc, values, nulls);
898 
899  /* make the tuple into a datum */
900  result = HeapTupleGetDatum(tuple);
901 
902  /* clean up */
903  /* For Windowing functions, it is important to leave */
904  /* the state intact, knowing that the aggcontext will be */
905  /* freed by PgSQL when the statement is complete. */
906  /* https://trac.osgeo.org/postgis/ticket/4770 */
907  // rtpg_summarystats_arg_destroy(state);
908 
909  PG_RETURN_DATUM(result);
910 }
911 
912 #undef VALUES_LENGTH
913 #define VALUES_LENGTH 4
914 
918 PG_FUNCTION_INFO_V1(RASTER_histogram);
919 Datum RASTER_histogram(PG_FUNCTION_ARGS)
920 {
921  FuncCallContext *funcctx;
922  TupleDesc tupdesc;
923 
924  int i;
925  rt_histogram hist;
926  rt_histogram hist2;
927  int call_cntr;
928  int max_calls;
929 
930  /* first call of function */
931  if (SRF_IS_FIRSTCALL()) {
932  MemoryContext oldcontext;
933 
934  rt_pgraster *pgraster = NULL;
935  rt_raster raster = NULL;
936  rt_band band = NULL;
937  int32_t bandindex = 1;
938  int num_bands = 0;
939  bool exclude_nodata_value = TRUE;
940  double sample = 0;
941  uint32_t bin_count = 0;
942  double *bin_width = NULL;
943  uint32_t bin_width_count = 0;
944  double width = 0;
945  bool right = FALSE;
946  double min = 0;
947  double max = 0;
948  rt_bandstats stats = NULL;
949  uint32_t count;
950 
951  int j;
952  int n;
953 
954  ArrayType *array;
955  Oid etype;
956  Datum *e;
957  bool *nulls;
958  int16 typlen;
959  bool typbyval;
960  char typalign;
961 
962  POSTGIS_RT_DEBUG(3, "RASTER_histogram: Starting");
963 
964  /* create a function context for cross-call persistence */
965  funcctx = SRF_FIRSTCALL_INIT();
966 
967  /* switch to memory context appropriate for multiple function calls */
968  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
969 
970  /* pgraster is null, return nothing */
971  if (PG_ARGISNULL(0)) {
972  MemoryContextSwitchTo(oldcontext);
973  SRF_RETURN_DONE(funcctx);
974  }
975  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
976 
977  raster = rt_raster_deserialize(pgraster, FALSE);
978  if (!raster) {
979  PG_FREE_IF_COPY(pgraster, 0);
980  MemoryContextSwitchTo(oldcontext);
981  elog(ERROR, "RASTER_histogram: Cannot deserialize raster");
982  SRF_RETURN_DONE(funcctx);
983  }
984 
985  /* band index is 1-based */
986  if (!PG_ARGISNULL(1))
987  bandindex = PG_GETARG_INT32(1);
988  num_bands = rt_raster_get_num_bands(raster);
989  if (bandindex < 1 || bandindex > num_bands) {
990  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
991  rt_raster_destroy(raster);
992  PG_FREE_IF_COPY(pgraster, 0);
993  MemoryContextSwitchTo(oldcontext);
994  SRF_RETURN_DONE(funcctx);
995  }
996 
997  /* exclude_nodata_value flag */
998  if (!PG_ARGISNULL(2))
999  exclude_nodata_value = PG_GETARG_BOOL(2);
1000 
1001  /* sample % */
1002  if (!PG_ARGISNULL(3)) {
1003  sample = PG_GETARG_FLOAT8(3);
1004  if (sample < 0 || sample > 1) {
1005  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
1006  rt_raster_destroy(raster);
1007  PG_FREE_IF_COPY(pgraster, 0);
1008  MemoryContextSwitchTo(oldcontext);
1009  SRF_RETURN_DONE(funcctx);
1010  }
1011  else if (FLT_EQ(sample, 0.0))
1012  sample = 1;
1013  }
1014  else
1015  sample = 1;
1016 
1017  /* bin_count */
1018  if (!PG_ARGISNULL(4)) {
1019  bin_count = PG_GETARG_INT32(4);
1020  if (bin_count < 1) bin_count = 0;
1021  }
1022 
1023  /* bin_width */
1024  if (!PG_ARGISNULL(5)) {
1025  array = PG_GETARG_ARRAYTYPE_P(5);
1026  etype = ARR_ELEMTYPE(array);
1027  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1028 
1029  switch (etype) {
1030  case FLOAT4OID:
1031  case FLOAT8OID:
1032  break;
1033  default:
1034  rt_raster_destroy(raster);
1035  PG_FREE_IF_COPY(pgraster, 0);
1036  MemoryContextSwitchTo(oldcontext);
1037  elog(ERROR, "RASTER_histogram: Invalid data type for width");
1038  SRF_RETURN_DONE(funcctx);
1039  break;
1040  }
1041 
1042  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1043  &nulls, &n);
1044 
1045  bin_width = palloc(sizeof(double) * n);
1046  for (i = 0, j = 0; i < n; i++) {
1047  if (nulls[i]) continue;
1048 
1049  switch (etype) {
1050  case FLOAT4OID:
1051  width = (double) DatumGetFloat4(e[i]);
1052  break;
1053  case FLOAT8OID:
1054  width = (double) DatumGetFloat8(e[i]);
1055  break;
1056  }
1057 
1058  if (width < 0 || FLT_EQ(width, 0.0)) {
1059  elog(NOTICE, "Invalid value for width (must be greater than 0). Returning NULL");
1060  pfree(bin_width);
1061  rt_raster_destroy(raster);
1062  PG_FREE_IF_COPY(pgraster, 0);
1063  MemoryContextSwitchTo(oldcontext);
1064  SRF_RETURN_DONE(funcctx);
1065  }
1066 
1067  bin_width[j] = width;
1068  POSTGIS_RT_DEBUGF(5, "bin_width[%d] = %f", j, bin_width[j]);
1069  j++;
1070  }
1071  bin_width_count = j;
1072 
1073  if (j < 1) {
1074  pfree(bin_width);
1075  bin_width = NULL;
1076  }
1077  }
1078 
1079  /* right */
1080  if (!PG_ARGISNULL(6))
1081  right = PG_GETARG_BOOL(6);
1082 
1083  /* min */
1084  if (!PG_ARGISNULL(7)) min = PG_GETARG_FLOAT8(7);
1085 
1086  /* max */
1087  if (!PG_ARGISNULL(8)) max = PG_GETARG_FLOAT8(8);
1088 
1089  /* get band */
1090  band = rt_raster_get_band(raster, bandindex - 1);
1091  if (!band) {
1092  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
1093  rt_raster_destroy(raster);
1094  PG_FREE_IF_COPY(pgraster, 0);
1095  MemoryContextSwitchTo(oldcontext);
1096  SRF_RETURN_DONE(funcctx);
1097  }
1098 
1099  /* get stats */
1100  stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 1, NULL, NULL, NULL);
1101  rt_band_destroy(band);
1102  rt_raster_destroy(raster);
1103  PG_FREE_IF_COPY(pgraster, 0);
1104  if (NULL == stats || NULL == stats->values) {
1105  elog(NOTICE, "Cannot compute summary statistics for band at index %d", bandindex);
1106  MemoryContextSwitchTo(oldcontext);
1107  SRF_RETURN_DONE(funcctx);
1108  }
1109  else if (stats->count < 1) {
1110  elog(NOTICE, "Cannot compute histogram for band at index %d as the band has no values", bandindex);
1111  MemoryContextSwitchTo(oldcontext);
1112  SRF_RETURN_DONE(funcctx);
1113  }
1114 
1115  /* get histogram */
1116  hist = rt_band_get_histogram(stats, (uint32_t)bin_count, bin_width, bin_width_count, right, min, max, &count);
1117  if (bin_width_count) pfree(bin_width);
1118  pfree(stats);
1119  if (NULL == hist || !count) {
1120  elog(NOTICE, "Cannot compute histogram for band at index %d", bandindex);
1121  MemoryContextSwitchTo(oldcontext);
1122  SRF_RETURN_DONE(funcctx);
1123  }
1124 
1125  POSTGIS_RT_DEBUGF(3, "%d bins returned", count);
1126 
1127  /* Store needed information */
1128  funcctx->user_fctx = hist;
1129 
1130  /* total number of tuples to be returned */
1131  funcctx->max_calls = count;
1132 
1133  /* Build a tuple descriptor for our result type */
1134  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
1135  ereport(ERROR, (
1136  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1137  errmsg(
1138  "function returning record called in context "
1139  "that cannot accept type record"
1140  )
1141  ));
1142  }
1143 
1144  BlessTupleDesc(tupdesc);
1145  funcctx->tuple_desc = tupdesc;
1146 
1147  MemoryContextSwitchTo(oldcontext);
1148  }
1149 
1150  /* stuff done on every call of the function */
1151  funcctx = SRF_PERCALL_SETUP();
1152 
1153  call_cntr = funcctx->call_cntr;
1154  max_calls = funcctx->max_calls;
1155  tupdesc = funcctx->tuple_desc;
1156  hist2 = funcctx->user_fctx;
1157 
1158  /* do when there is more left to send */
1159  if (call_cntr < max_calls) {
1160  Datum values[VALUES_LENGTH];
1161  bool nulls[VALUES_LENGTH];
1162  HeapTuple tuple;
1163  Datum result;
1164 
1165  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
1166 
1167  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
1168 
1169  values[0] = Float8GetDatum(hist2[call_cntr].min);
1170  values[1] = Float8GetDatum(hist2[call_cntr].max);
1171  values[2] = Int64GetDatum(hist2[call_cntr].count);
1172  values[3] = Float8GetDatum(hist2[call_cntr].percent);
1173 
1174  /* build a tuple */
1175  tuple = heap_form_tuple(tupdesc, values, nulls);
1176 
1177  /* make the tuple into a datum */
1178  result = HeapTupleGetDatum(tuple);
1179 
1180  SRF_RETURN_NEXT(funcctx, result);
1181  }
1182  /* do when there is no more left */
1183  else {
1184  pfree(hist2);
1185  SRF_RETURN_DONE(funcctx);
1186  }
1187 }
1188 
1189 #undef VALUES_LENGTH
1190 #define VALUES_LENGTH 2
1191 
1195 PG_FUNCTION_INFO_V1(RASTER_quantile);
1196 Datum RASTER_quantile(PG_FUNCTION_ARGS)
1197 {
1198  FuncCallContext *funcctx;
1199  TupleDesc tupdesc;
1200 
1201  int i;
1202  rt_quantile quant;
1203  rt_quantile quant2;
1204  int call_cntr;
1205  int max_calls;
1206 
1207  /* first call of function */
1208  if (SRF_IS_FIRSTCALL()) {
1209  MemoryContext oldcontext;
1210 
1211  rt_pgraster *pgraster = NULL;
1212  rt_raster raster = NULL;
1213  rt_band band = NULL;
1214  int32_t bandindex = 0;
1215  int num_bands = 0;
1216  bool exclude_nodata_value = TRUE;
1217  double sample = 0;
1218  double *quantiles = NULL;
1219  uint32_t quantiles_count = 0;
1220  double quantile = 0;
1221  rt_bandstats stats = NULL;
1222  uint32_t count;
1223 
1224  int j;
1225  int n;
1226 
1227  ArrayType *array;
1228  Oid etype;
1229  Datum *e;
1230  bool *nulls;
1231  int16 typlen;
1232  bool typbyval;
1233  char typalign;
1234 
1235  /* create a function context for cross-call persistence */
1236  funcctx = SRF_FIRSTCALL_INIT();
1237 
1238  /* switch to memory context appropriate for multiple function calls */
1239  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
1240 
1241  /* pgraster is null, return nothing */
1242  if (PG_ARGISNULL(0)) {
1243  MemoryContextSwitchTo(oldcontext);
1244  SRF_RETURN_DONE(funcctx);
1245  }
1246  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1247 
1248  raster = rt_raster_deserialize(pgraster, FALSE);
1249  if (!raster) {
1250  PG_FREE_IF_COPY(pgraster, 0);
1251  MemoryContextSwitchTo(oldcontext);
1252  elog(ERROR, "RASTER_quantile: Cannot deserialize raster");
1253  SRF_RETURN_DONE(funcctx);
1254  }
1255 
1256  /* band index is 1-based */
1257  bandindex = PG_GETARG_INT32(1);
1258  num_bands = rt_raster_get_num_bands(raster);
1259  if (bandindex < 1 || bandindex > num_bands) {
1260  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
1261  rt_raster_destroy(raster);
1262  PG_FREE_IF_COPY(pgraster, 0);
1263  MemoryContextSwitchTo(oldcontext);
1264  SRF_RETURN_DONE(funcctx);
1265  }
1266 
1267  /* exclude_nodata_value flag */
1268  if (!PG_ARGISNULL(2))
1269  exclude_nodata_value = PG_GETARG_BOOL(2);
1270 
1271  /* sample % */
1272  if (!PG_ARGISNULL(3)) {
1273  sample = PG_GETARG_FLOAT8(3);
1274  if (sample < 0 || sample > 1) {
1275  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
1276  rt_raster_destroy(raster);
1277  PG_FREE_IF_COPY(pgraster, 0);
1278  MemoryContextSwitchTo(oldcontext);
1279  SRF_RETURN_DONE(funcctx);
1280  }
1281  else if (FLT_EQ(sample, 0.0))
1282  sample = 1;
1283  }
1284  else
1285  sample = 1;
1286 
1287  /* quantiles */
1288  if (!PG_ARGISNULL(4)) {
1289  array = PG_GETARG_ARRAYTYPE_P(4);
1290  etype = ARR_ELEMTYPE(array);
1291  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1292 
1293  switch (etype) {
1294  case FLOAT4OID:
1295  case FLOAT8OID:
1296  break;
1297  default:
1298  rt_raster_destroy(raster);
1299  PG_FREE_IF_COPY(pgraster, 0);
1300  MemoryContextSwitchTo(oldcontext);
1301  elog(ERROR, "RASTER_quantile: Invalid data type for quantiles");
1302  SRF_RETURN_DONE(funcctx);
1303  break;
1304  }
1305 
1306  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1307  &nulls, &n);
1308 
1309  quantiles = palloc(sizeof(double) * n);
1310  for (i = 0, j = 0; i < n; i++) {
1311  if (nulls[i]) continue;
1312 
1313  switch (etype) {
1314  case FLOAT4OID:
1315  quantile = (double) DatumGetFloat4(e[i]);
1316  break;
1317  case FLOAT8OID:
1318  quantile = (double) DatumGetFloat8(e[i]);
1319  break;
1320  }
1321 
1322  if (quantile < 0 || quantile > 1) {
1323  elog(NOTICE, "Invalid value for quantile (must be between 0 and 1). Returning NULL");
1324  pfree(quantiles);
1325  rt_raster_destroy(raster);
1326  PG_FREE_IF_COPY(pgraster, 0);
1327  MemoryContextSwitchTo(oldcontext);
1328  SRF_RETURN_DONE(funcctx);
1329  }
1330 
1331  quantiles[j] = quantile;
1332  POSTGIS_RT_DEBUGF(5, "quantiles[%d] = %f", j, quantiles[j]);
1333  j++;
1334  }
1335  quantiles_count = j;
1336 
1337  if (j < 1) {
1338  pfree(quantiles);
1339  quantiles = NULL;
1340  }
1341  }
1342 
1343  /* get band */
1344  band = rt_raster_get_band(raster, bandindex - 1);
1345  if (!band) {
1346  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
1347  rt_raster_destroy(raster);
1348  PG_FREE_IF_COPY(pgraster, 0);
1349  MemoryContextSwitchTo(oldcontext);
1350  SRF_RETURN_DONE(funcctx);
1351  }
1352 
1353  /* get stats */
1354  stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 1, NULL, NULL, NULL);
1355  rt_band_destroy(band);
1356  rt_raster_destroy(raster);
1357  PG_FREE_IF_COPY(pgraster, 0);
1358  if (NULL == stats || NULL == stats->values) {
1359  elog(NOTICE, "Cannot retrieve summary statistics for band at index %d", bandindex);
1360  MemoryContextSwitchTo(oldcontext);
1361  SRF_RETURN_DONE(funcctx);
1362  }
1363  else if (stats->count < 1) {
1364  elog(NOTICE, "Cannot compute quantiles for band at index %d as the band has no values", bandindex);
1365  MemoryContextSwitchTo(oldcontext);
1366  SRF_RETURN_DONE(funcctx);
1367  }
1368 
1369  /* get quantiles */
1370  quant = rt_band_get_quantiles(stats, quantiles, quantiles_count, &count);
1371  if (quantiles_count) pfree(quantiles);
1372  pfree(stats);
1373  if (NULL == quant || !count) {
1374  elog(NOTICE, "Cannot compute quantiles for band at index %d", bandindex);
1375  MemoryContextSwitchTo(oldcontext);
1376  SRF_RETURN_DONE(funcctx);
1377  }
1378 
1379  POSTGIS_RT_DEBUGF(3, "%d quantiles returned", count);
1380 
1381  /* Store needed information */
1382  funcctx->user_fctx = quant;
1383 
1384  /* total number of tuples to be returned */
1385  funcctx->max_calls = count;
1386 
1387  /* Build a tuple descriptor for our result type */
1388  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
1389  ereport(ERROR, (
1390  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1391  errmsg(
1392  "function returning record called in context "
1393  "that cannot accept type record"
1394  )
1395  ));
1396  }
1397 
1398  BlessTupleDesc(tupdesc);
1399  funcctx->tuple_desc = tupdesc;
1400 
1401  MemoryContextSwitchTo(oldcontext);
1402  }
1403 
1404  /* stuff done on every call of the function */
1405  funcctx = SRF_PERCALL_SETUP();
1406 
1407  call_cntr = funcctx->call_cntr;
1408  max_calls = funcctx->max_calls;
1409  tupdesc = funcctx->tuple_desc;
1410  quant2 = funcctx->user_fctx;
1411 
1412  /* do when there is more left to send */
1413  if (call_cntr < max_calls) {
1414  Datum values[VALUES_LENGTH];
1415  bool nulls[VALUES_LENGTH];
1416  HeapTuple tuple;
1417  Datum result;
1418 
1419  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
1420 
1421  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
1422 
1423  values[0] = Float8GetDatum(quant2[call_cntr].quantile);
1424  values[1] = Float8GetDatum(quant2[call_cntr].value);
1425 
1426  /* build a tuple */
1427  tuple = heap_form_tuple(tupdesc, values, nulls);
1428 
1429  /* make the tuple into a datum */
1430  result = HeapTupleGetDatum(tuple);
1431 
1432  SRF_RETURN_NEXT(funcctx, result);
1433  }
1434  /* do when there is no more left */
1435  else {
1436  pfree(quant2);
1437  SRF_RETURN_DONE(funcctx);
1438  }
1439 }
1440 
1441 #undef VALUES_LENGTH
1442 #define VALUES_LENGTH 3
1443 
1444 /* get counts of values */
1445 PG_FUNCTION_INFO_V1(RASTER_valueCount);
1446 Datum RASTER_valueCount(PG_FUNCTION_ARGS) {
1447  FuncCallContext *funcctx;
1448  TupleDesc tupdesc;
1449 
1450  int i;
1451  rt_valuecount vcnts;
1452  rt_valuecount vcnts2;
1453  int call_cntr;
1454  int max_calls;
1455 
1456  /* first call of function */
1457  if (SRF_IS_FIRSTCALL()) {
1458  MemoryContext oldcontext;
1459 
1460  rt_pgraster *pgraster = NULL;
1461  rt_raster raster = NULL;
1462  rt_band band = NULL;
1463  int32_t bandindex = 0;
1464  int num_bands = 0;
1465  bool exclude_nodata_value = TRUE;
1466  double *search_values = NULL;
1467  uint32_t search_values_count = 0;
1468  double roundto = 0;
1469  uint32_t count;
1470 
1471  int j;
1472  int n;
1473 
1474  ArrayType *array;
1475  Oid etype;
1476  Datum *e;
1477  bool *nulls;
1478  int16 typlen;
1479  bool typbyval;
1480  char typalign;
1481 
1482  /* create a function context for cross-call persistence */
1483  funcctx = SRF_FIRSTCALL_INIT();
1484 
1485  /* switch to memory context appropriate for multiple function calls */
1486  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
1487 
1488  /* pgraster is null, return nothing */
1489  if (PG_ARGISNULL(0)) {
1490  MemoryContextSwitchTo(oldcontext);
1491  SRF_RETURN_DONE(funcctx);
1492  }
1493  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1494 
1495  raster = rt_raster_deserialize(pgraster, FALSE);
1496  if (!raster) {
1497  PG_FREE_IF_COPY(pgraster, 0);
1498  MemoryContextSwitchTo(oldcontext);
1499  elog(ERROR, "RASTER_valueCount: Cannot deserialize raster");
1500  SRF_RETURN_DONE(funcctx);
1501  }
1502 
1503  /* band index is 1-based */
1504  bandindex = PG_GETARG_INT32(1);
1505  num_bands = rt_raster_get_num_bands(raster);
1506  if (bandindex < 1 || bandindex > num_bands) {
1507  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
1508  rt_raster_destroy(raster);
1509  PG_FREE_IF_COPY(pgraster, 0);
1510  MemoryContextSwitchTo(oldcontext);
1511  SRF_RETURN_DONE(funcctx);
1512  }
1513 
1514  /* exclude_nodata_value flag */
1515  if (!PG_ARGISNULL(2))
1516  exclude_nodata_value = PG_GETARG_BOOL(2);
1517 
1518  /* search values */
1519  if (!PG_ARGISNULL(3)) {
1520  array = PG_GETARG_ARRAYTYPE_P(3);
1521  etype = ARR_ELEMTYPE(array);
1522  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1523 
1524  switch (etype) {
1525  case FLOAT4OID:
1526  case FLOAT8OID:
1527  break;
1528  default:
1529  rt_raster_destroy(raster);
1530  PG_FREE_IF_COPY(pgraster, 0);
1531  MemoryContextSwitchTo(oldcontext);
1532  elog(ERROR, "RASTER_valueCount: Invalid data type for values");
1533  SRF_RETURN_DONE(funcctx);
1534  break;
1535  }
1536 
1537  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1538  &nulls, &n);
1539 
1540  search_values = palloc(sizeof(double) * n);
1541  for (i = 0, j = 0; i < n; i++) {
1542  if (nulls[i]) continue;
1543 
1544  switch (etype) {
1545  case FLOAT4OID:
1546  search_values[j] = (double) DatumGetFloat4(e[i]);
1547  break;
1548  case FLOAT8OID:
1549  search_values[j] = (double) DatumGetFloat8(e[i]);
1550  break;
1551  }
1552 
1553  POSTGIS_RT_DEBUGF(5, "search_values[%d] = %f", j, search_values[j]);
1554  j++;
1555  }
1556  search_values_count = j;
1557 
1558  if (j < 1) {
1559  pfree(search_values);
1560  search_values = NULL;
1561  }
1562  }
1563 
1564  /* roundto */
1565  if (!PG_ARGISNULL(4)) {
1566  roundto = PG_GETARG_FLOAT8(4);
1567  if (roundto < 0.) roundto = 0;
1568  }
1569 
1570  /* get band */
1571  band = rt_raster_get_band(raster, bandindex - 1);
1572  if (!band) {
1573  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
1574  rt_raster_destroy(raster);
1575  PG_FREE_IF_COPY(pgraster, 0);
1576  MemoryContextSwitchTo(oldcontext);
1577  SRF_RETURN_DONE(funcctx);
1578  }
1579 
1580  /* get counts of values */
1581  vcnts = rt_band_get_value_count(band, (int) exclude_nodata_value, search_values, search_values_count, roundto, NULL, &count);
1582  rt_band_destroy(band);
1583  rt_raster_destroy(raster);
1584  PG_FREE_IF_COPY(pgraster, 0);
1585  if (NULL == vcnts || !count) {
1586  elog(NOTICE, "Cannot count the values for band at index %d", bandindex);
1587  MemoryContextSwitchTo(oldcontext);
1588  SRF_RETURN_DONE(funcctx);
1589  }
1590 
1591  POSTGIS_RT_DEBUGF(3, "%d value counts returned", count);
1592 
1593  /* Store needed information */
1594  funcctx->user_fctx = vcnts;
1595 
1596  /* total number of tuples to be returned */
1597  funcctx->max_calls = count;
1598 
1599  /* Build a tuple descriptor for our result type */
1600  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
1601  ereport(ERROR, (
1602  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1603  errmsg(
1604  "function returning record called in context "
1605  "that cannot accept type record"
1606  )
1607  ));
1608  }
1609 
1610  BlessTupleDesc(tupdesc);
1611  funcctx->tuple_desc = tupdesc;
1612 
1613  MemoryContextSwitchTo(oldcontext);
1614  }
1615 
1616  /* stuff done on every call of the function */
1617  funcctx = SRF_PERCALL_SETUP();
1618 
1619  call_cntr = funcctx->call_cntr;
1620  max_calls = funcctx->max_calls;
1621  tupdesc = funcctx->tuple_desc;
1622  vcnts2 = funcctx->user_fctx;
1623 
1624  /* do when there is more left to send */
1625  if (call_cntr < max_calls) {
1626  Datum values[VALUES_LENGTH];
1627  bool nulls[VALUES_LENGTH];
1628  HeapTuple tuple;
1629  Datum result;
1630 
1631  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
1632 
1633  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
1634 
1635  values[0] = Float8GetDatum(vcnts2[call_cntr].value);
1636  values[1] = UInt32GetDatum(vcnts2[call_cntr].count);
1637  values[2] = Float8GetDatum(vcnts2[call_cntr].percent);
1638 
1639  /* build a tuple */
1640  tuple = heap_form_tuple(tupdesc, values, nulls);
1641 
1642  /* make the tuple into a datum */
1643  result = HeapTupleGetDatum(tuple);
1644 
1645  SRF_RETURN_NEXT(funcctx, result);
1646  }
1647  /* do when there is no more left */
1648  else {
1649  pfree(vcnts2);
1650  SRF_RETURN_DONE(funcctx);
1651  }
1652 }
1653 
1654 /* get counts of values for a coverage */
1655 PG_FUNCTION_INFO_V1(RASTER_valueCountCoverage);
1656 Datum RASTER_valueCountCoverage(PG_FUNCTION_ARGS) {
1657  FuncCallContext *funcctx;
1658  TupleDesc tupdesc;
1659 
1660  uint32_t i;
1661  uint64_t covcount = 0;
1662  uint64_t covtotal = 0;
1663  rt_valuecount covvcnts = NULL;
1664  rt_valuecount covvcnts2;
1665  int call_cntr;
1666  int max_calls;
1667 
1668  POSTGIS_RT_DEBUG(3, "RASTER_valueCountCoverage: Starting");
1669 
1670  /* first call of function */
1671  if (SRF_IS_FIRSTCALL()) {
1672  MemoryContext oldcontext;
1673 
1674  text *tablenametext = NULL;
1675  char *tablename = NULL;
1676  text *colnametext = NULL;
1677  char *colname = NULL;
1678  int32_t bandindex = 1;
1679  bool exclude_nodata_value = TRUE;
1680  double *search_values = NULL;
1681  uint32_t search_values_count = 0;
1682  double roundto = 0;
1683 
1684  int len = 0;
1685  char *sql = NULL;
1686  int spi_result;
1687  Portal portal;
1688  HeapTuple tuple;
1689  Datum datum;
1690  bool isNull = FALSE;
1691  rt_pgraster *pgraster = NULL;
1692  rt_raster raster = NULL;
1693  rt_band band = NULL;
1694  int num_bands = 0;
1695  uint32_t count;
1696  uint32_t total;
1697  rt_valuecount vcnts = NULL;
1698  int exists = 0;
1699 
1700  uint32_t j;
1701  int n;
1702 
1703  ArrayType *array;
1704  Oid etype;
1705  Datum *e;
1706  bool *nulls;
1707  int16 typlen;
1708  bool typbyval;
1709  char typalign;
1710 
1711  /* create a function context for cross-call persistence */
1712  funcctx = SRF_FIRSTCALL_INIT();
1713 
1714  /* switch to memory context appropriate for multiple function calls */
1715  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
1716 
1717  /* tablename is null, return null */
1718  if (PG_ARGISNULL(0)) {
1719  elog(NOTICE, "Table name must be provided");
1720  MemoryContextSwitchTo(oldcontext);
1721  SRF_RETURN_DONE(funcctx);
1722  }
1723  tablenametext = PG_GETARG_TEXT_P(0);
1724  tablename = text_to_cstring(tablenametext);
1725  if (!strlen(tablename)) {
1726  elog(NOTICE, "Table name must be provided");
1727  MemoryContextSwitchTo(oldcontext);
1728  SRF_RETURN_DONE(funcctx);
1729  }
1730  POSTGIS_RT_DEBUGF(3, "tablename = %s", tablename);
1731 
1732  /* column name is null, return null */
1733  if (PG_ARGISNULL(1)) {
1734  elog(NOTICE, "Column name must be provided");
1735  MemoryContextSwitchTo(oldcontext);
1736  SRF_RETURN_DONE(funcctx);
1737  }
1738  colnametext = PG_GETARG_TEXT_P(1);
1739  colname = text_to_cstring(colnametext);
1740  if (!strlen(colname)) {
1741  elog(NOTICE, "Column name must be provided");
1742  MemoryContextSwitchTo(oldcontext);
1743  SRF_RETURN_DONE(funcctx);
1744  }
1745  POSTGIS_RT_DEBUGF(3, "colname = %s", colname);
1746 
1747  /* band index is 1-based */
1748  if (!PG_ARGISNULL(2))
1749  bandindex = PG_GETARG_INT32(2);
1750 
1751  /* exclude_nodata_value flag */
1752  if (!PG_ARGISNULL(3))
1753  exclude_nodata_value = PG_GETARG_BOOL(3);
1754 
1755  /* search values */
1756  if (!PG_ARGISNULL(4)) {
1757  array = PG_GETARG_ARRAYTYPE_P(4);
1758  etype = ARR_ELEMTYPE(array);
1759  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1760 
1761  switch (etype) {
1762  case FLOAT4OID:
1763  case FLOAT8OID:
1764  break;
1765  default:
1766  MemoryContextSwitchTo(oldcontext);
1767  elog(ERROR, "RASTER_valueCountCoverage: Invalid data type for values");
1768  SRF_RETURN_DONE(funcctx);
1769  break;
1770  }
1771 
1772  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1773  &nulls, &n);
1774 
1775  search_values = palloc(sizeof(double) * n);
1776  for (i = 0, j = 0; i < (uint32_t) n; i++) {
1777  if (nulls[i]) continue;
1778 
1779  switch (etype) {
1780  case FLOAT4OID:
1781  search_values[j] = (double) DatumGetFloat4(e[i]);
1782  break;
1783  case FLOAT8OID:
1784  search_values[j] = (double) DatumGetFloat8(e[i]);
1785  break;
1786  }
1787 
1788  POSTGIS_RT_DEBUGF(5, "search_values[%d] = %f", j, search_values[j]);
1789  j++;
1790  }
1791  search_values_count = j;
1792 
1793  if (j < 1) {
1794  pfree(search_values);
1795  search_values = NULL;
1796  }
1797  }
1798 
1799  /* roundto */
1800  if (!PG_ARGISNULL(5)) {
1801  roundto = PG_GETARG_FLOAT8(5);
1802  if (roundto < 0.) roundto = 0;
1803  }
1804 
1805  /* iterate through rasters of coverage */
1806  /* connect to database */
1807  spi_result = SPI_connect();
1808  if (spi_result != SPI_OK_CONNECT) {
1809 
1810  if (search_values_count) pfree(search_values);
1811 
1812  MemoryContextSwitchTo(oldcontext);
1813  elog(ERROR, "RASTER_valueCountCoverage: Cannot connect to database using SPI");
1814  SRF_RETURN_DONE(funcctx);
1815  }
1816 
1817  /* create sql */
1818  len = sizeof(char) * (strlen("SELECT \"\" FROM \"\" WHERE \"\" IS NOT NULL") + (strlen(colname) * 2) + strlen(tablename) + 1);
1819  sql = (char *) palloc(len);
1820  if (NULL == sql) {
1821 
1822  if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
1823  SPI_finish();
1824 
1825  if (search_values_count) pfree(search_values);
1826 
1827  MemoryContextSwitchTo(oldcontext);
1828  elog(ERROR, "RASTER_valueCountCoverage: Cannot allocate memory for sql");
1829  SRF_RETURN_DONE(funcctx);
1830  }
1831 
1832  /* get cursor */
1833  snprintf(sql, len, "SELECT \"%s\" FROM \"%s\" WHERE \"%s\" IS NOT NULL", colname, tablename, colname);
1834  POSTGIS_RT_DEBUGF(3, "RASTER_valueCountCoverage: %s", sql);
1835  portal = SPI_cursor_open_with_args(
1836  "coverage",
1837  sql,
1838  0, NULL,
1839  NULL, NULL,
1840  TRUE, 0
1841  );
1842  pfree(sql);
1843 
1844  /* process resultset */
1845  SPI_cursor_fetch(portal, TRUE, 1);
1846  while (SPI_processed == 1 && SPI_tuptable != NULL) {
1847  tupdesc = SPI_tuptable->tupdesc;
1848  tuple = SPI_tuptable->vals[0];
1849 
1850  datum = SPI_getbinval(tuple, tupdesc, 1, &isNull);
1851  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1852 
1853  SPI_freetuptable(SPI_tuptable);
1854  SPI_cursor_close(portal);
1855  SPI_finish();
1856 
1857  if (NULL != covvcnts) pfree(covvcnts);
1858  if (search_values_count) pfree(search_values);
1859 
1860  MemoryContextSwitchTo(oldcontext);
1861  elog(ERROR, "RASTER_valueCountCoverage: Cannot get raster of coverage");
1862  SRF_RETURN_DONE(funcctx);
1863  }
1864  else if (isNull) {
1865  SPI_cursor_fetch(portal, TRUE, 1);
1866  continue;
1867  }
1868 
1869  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(datum);
1870 
1871  raster = rt_raster_deserialize(pgraster, FALSE);
1872  if (!raster) {
1873 
1874  SPI_freetuptable(SPI_tuptable);
1875  SPI_cursor_close(portal);
1876  SPI_finish();
1877 
1878  if (NULL != covvcnts) pfree(covvcnts);
1879  if (search_values_count) pfree(search_values);
1880 
1881  MemoryContextSwitchTo(oldcontext);
1882  elog(ERROR, "RASTER_valueCountCoverage: Cannot deserialize raster");
1883  SRF_RETURN_DONE(funcctx);
1884  }
1885 
1886  /* inspect number of bands*/
1887  num_bands = rt_raster_get_num_bands(raster);
1888  if (bandindex < 1 || bandindex > num_bands) {
1889  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
1890 
1891  rt_raster_destroy(raster);
1892 
1893  SPI_freetuptable(SPI_tuptable);
1894  SPI_cursor_close(portal);
1895  SPI_finish();
1896 
1897  if (NULL != covvcnts) pfree(covvcnts);
1898  if (search_values_count) pfree(search_values);
1899 
1900  MemoryContextSwitchTo(oldcontext);
1901  SRF_RETURN_DONE(funcctx);
1902  }
1903 
1904  /* get band */
1905  band = rt_raster_get_band(raster, bandindex - 1);
1906  if (!band) {
1907  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
1908 
1909  rt_raster_destroy(raster);
1910 
1911  SPI_freetuptable(SPI_tuptable);
1912  SPI_cursor_close(portal);
1913  SPI_finish();
1914 
1915  if (NULL != covvcnts) pfree(covvcnts);
1916  if (search_values_count) pfree(search_values);
1917 
1918  MemoryContextSwitchTo(oldcontext);
1919  SRF_RETURN_DONE(funcctx);
1920  }
1921 
1922  /* get counts of values */
1923  vcnts = rt_band_get_value_count(band, (int) exclude_nodata_value, search_values, search_values_count, roundto, &total, &count);
1924  rt_band_destroy(band);
1925  rt_raster_destroy(raster);
1926  if (NULL == vcnts || !count) {
1927  elog(NOTICE, "Cannot count the values for band at index %d", bandindex);
1928 
1929  SPI_freetuptable(SPI_tuptable);
1930  SPI_cursor_close(portal);
1931  SPI_finish();
1932 
1933  if (NULL != covvcnts) free(covvcnts);
1934  if (search_values_count) pfree(search_values);
1935 
1936  MemoryContextSwitchTo(oldcontext);
1937  SRF_RETURN_DONE(funcctx);
1938  }
1939 
1940  POSTGIS_RT_DEBUGF(3, "%d value counts returned", count);
1941 
1942  if (NULL == covvcnts) {
1943  covvcnts = (rt_valuecount) SPI_palloc(sizeof(struct rt_valuecount_t) * count);
1944  if (NULL == covvcnts) {
1945 
1946  SPI_freetuptable(SPI_tuptable);
1947  SPI_cursor_close(portal);
1948  SPI_finish();
1949 
1950  if (search_values_count) pfree(search_values);
1951 
1952  MemoryContextSwitchTo(oldcontext);
1953  elog(ERROR, "RASTER_valueCountCoverage: Cannot allocate memory for value counts of coverage");
1954  SRF_RETURN_DONE(funcctx);
1955  }
1956 
1957  for (i = 0; i < count; i++) {
1958  covvcnts[i].value = vcnts[i].value;
1959  covvcnts[i].count = vcnts[i].count;
1960  covvcnts[i].percent = -1;
1961  }
1962 
1963  covcount = count;
1964  }
1965  else {
1966  for (i = 0; i < count; i++) {
1967  exists = 0;
1968 
1969  for (j = 0; j < covcount; j++) {
1970  if (FLT_EQ(vcnts[i].value, covvcnts[j].value)) {
1971  exists = 1;
1972  break;
1973  }
1974  }
1975 
1976  if (exists) {
1977  covvcnts[j].count += vcnts[i].count;
1978  }
1979  else {
1980  covcount++;
1981  covvcnts = SPI_repalloc(covvcnts, sizeof(struct rt_valuecount_t) * covcount);
1982  if (!covvcnts) {
1983  SPI_freetuptable(SPI_tuptable);
1984  SPI_cursor_close(portal);
1985  SPI_finish();
1986 
1987  if (search_values_count) pfree(search_values);
1988 
1989  MemoryContextSwitchTo(oldcontext);
1990  elog(ERROR, "RASTER_valueCountCoverage: Cannot change allocated memory for value counts of coverage");
1991  SRF_RETURN_DONE(funcctx);
1992  }
1993 
1994  covvcnts[covcount - 1].value = vcnts[i].value;
1995  covvcnts[covcount - 1].count = vcnts[i].count;
1996  covvcnts[covcount - 1].percent = -1;
1997  }
1998  }
1999  }
2000 
2001  covtotal += total;
2002 
2003  pfree(vcnts);
2004 
2005  /* next record */
2006  SPI_cursor_fetch(portal, TRUE, 1);
2007  }
2008 
2009  if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
2010  SPI_cursor_close(portal);
2011  SPI_finish();
2012 
2013  if (search_values_count) pfree(search_values);
2014 
2015  /* compute percentages */
2016  for (i = 0; i < covcount; i++) {
2017  covvcnts[i].percent = (double) covvcnts[i].count / covtotal;
2018  }
2019 
2020  /* Store needed information */
2021  funcctx->user_fctx = covvcnts;
2022 
2023  /* total number of tuples to be returned */
2024  funcctx->max_calls = covcount;
2025 
2026  /* Build a tuple descriptor for our result type */
2027  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
2028  ereport(ERROR, (
2029  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
2030  errmsg(
2031  "function returning record called in context "
2032  "that cannot accept type record"
2033  )
2034  ));
2035  }
2036 
2037  BlessTupleDesc(tupdesc);
2038  funcctx->tuple_desc = tupdesc;
2039 
2040  MemoryContextSwitchTo(oldcontext);
2041  }
2042 
2043  /* stuff done on every call of the function */
2044  funcctx = SRF_PERCALL_SETUP();
2045 
2046  call_cntr = funcctx->call_cntr;
2047  max_calls = funcctx->max_calls;
2048  tupdesc = funcctx->tuple_desc;
2049  covvcnts2 = funcctx->user_fctx;
2050 
2051  /* do when there is more left to send */
2052  if (call_cntr < max_calls) {
2053  Datum values[VALUES_LENGTH];
2054  bool nulls[VALUES_LENGTH];
2055  HeapTuple tuple;
2056  Datum result;
2057 
2058  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
2059 
2060  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
2061 
2062  values[0] = Float8GetDatum(covvcnts2[call_cntr].value);
2063  values[1] = UInt32GetDatum(covvcnts2[call_cntr].count);
2064  values[2] = Float8GetDatum(covvcnts2[call_cntr].percent);
2065 
2066  /* build a tuple */
2067  tuple = heap_form_tuple(tupdesc, values, nulls);
2068 
2069  /* make the tuple into a datum */
2070  result = HeapTupleGetDatum(tuple);
2071 
2072  SRF_RETURN_NEXT(funcctx, result);
2073  }
2074  /* do when there is no more left */
2075  else {
2076  pfree(covvcnts2);
2077  SRF_RETURN_DONE(funcctx);
2078  }
2079 }
2080 
result
char result[OUT_DOUBLE_BUFFER_SIZE]
Definition: cu_print.c:262
TRUE
#define TRUE
Definition: dbfopen.c:73
FALSE
#define FALSE
Definition: dbfopen.c:72
rt_raster_destroy
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:86
FLT_EQ
#define FLT_EQ(x, y)
Definition: librtcore.h:2387
rt_bandstats
struct rt_bandstats_t * rt_bandstats
Definition: librtcore.h:152
rt_band_destroy
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition: rt_band.c:340
rt_raster_get_num_bands
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:376
rt_valuecount
struct rt_valuecount_t * rt_valuecount
Definition: librtcore.h:155
rt_band_get_histogram
rt_histogram rt_band_get_histogram(rt_bandstats stats, uint32_t bin_count, double *bin_widths, uint32_t bin_widths_count, int right, double min, double max, uint32_t *rtn_count)
Count the distribution of data.
Definition: rt_statistics.c:412
rt_band_get_value_count
rt_valuecount rt_band_get_value_count(rt_band band, int exclude_nodata_value, double *search_values, uint32_t search_values_count, double roundto, uint32_t *rtn_total, uint32_t *rtn_count)
Count the number of times provided value(s) occur in the band.
Definition: rt_statistics.c:1608
rt_band_get_quantiles
rt_quantile rt_band_get_quantiles(rt_bandstats stats, double *quantiles, int quantiles_count, uint32_t *rtn_count)
Compute the default set of or requested quantiles for a set of data the quantile formula used is same...
Definition: rt_statistics.c:672
rt_band_get_summary_stats
rt_bandstats rt_band_get_summary_stats(rt_band band, int exclude_nodata_value, double sample, int inc_vals, uint64_t *cK, double *cM, double *cQ)
Compute summary statistics for a band.
Definition: rt_statistics.c:110
rt_raster_deserialize
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition: rt_serialize.c:725
rt_raster_get_band
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:385
free
void free(void *)
genraster.value
int value
Definition: genraster.py:62
genraster.count
int count
Definition: genraster.py:57
ovdump.band
band
Definition: ovdump.py:58
rtrowdump.raster
raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
RASTER_summaryStats
Datum RASTER_summaryStats(PG_FUNCTION_ARGS)
Definition: rtpg_statistics.c:71
RASTER_summaryStatsCoverage
Datum RASTER_summaryStatsCoverage(PG_FUNCTION_ARGS)
Definition: rtpg_statistics.c:196
rtpg_summarystats_arg
struct rtpg_summarystats_arg_t * rtpg_summarystats_arg
Definition: rtpg_statistics.c:493
RASTER_histogram
Datum RASTER_histogram(PG_FUNCTION_ARGS)
Definition: rtpg_statistics.c:919
rtpg_summarystats_arg_destroy
static void rtpg_summarystats_arg_destroy(rtpg_summarystats_arg arg)
Definition: rtpg_statistics.c:508
RASTER_valueCountCoverage
Datum RASTER_valueCountCoverage(PG_FUNCTION_ARGS)
Definition: rtpg_statistics.c:1656
RASTER_valueCount
Datum RASTER_valueCount(PG_FUNCTION_ARGS)
Definition: rtpg_statistics.c:1446
rtpg_summarystats_arg_init
static rtpg_summarystats_arg rtpg_summarystats_arg_init()
Definition: rtpg_statistics.c:516
VALUES_LENGTH
#define VALUES_LENGTH
Definition: rtpg_statistics.c:1442
RASTER_summaryStats_transfn
Datum RASTER_summaryStats_transfn(PG_FUNCTION_ARGS)
Definition: rtpg_statistics.c:560
RASTER_quantile
Datum RASTER_quantile(PG_FUNCTION_ARGS)
Definition: rtpg_statistics.c:1196
RASTER_summaryStats_finalfn
Datum RASTER_summaryStats_finalfn(PG_FUNCTION_ARGS)
Definition: rtpg_statistics.c:824
PG_FUNCTION_INFO_V1
PG_FUNCTION_INFO_V1(RASTER_summaryStats)
Get summary stats of a band.
POSTGIS_RT_DEBUG
#define POSTGIS_RT_DEBUG(level, msg)
Definition: rtpostgis.h:65
POSTGIS_RT_DEBUGF
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition: rtpostgis.h:69
rt_band_t
Definition: librtcore.h:2466
rt_bandstats_t::mean
double mean
Definition: librtcore.h:2518
rt_bandstats_t::sorted
int sorted
Definition: librtcore.h:2522
rt_bandstats_t::sample
double sample
Definition: librtcore.h:2512
rt_bandstats_t::sum
double sum
Definition: librtcore.h:2517
rt_bandstats_t::max
double max
Definition: librtcore.h:2516
rt_bandstats_t::count
uint32_t count
Definition: librtcore.h:2513
rt_bandstats_t::min
double min
Definition: librtcore.h:2515
rt_bandstats_t::values
double * values
Definition: librtcore.h:2521
rt_bandstats_t::stddev
double stddev
Definition: librtcore.h:2519
rt_bandstats_t
Definition: librtcore.h:2511
rt_histogram_t
Definition: librtcore.h:2526
rt_quantile_t
Definition: librtcore.h:2538
rt_raster_serialized_t
Struct definitions.
Definition: librtcore.h:2403
rt_raster_t
Definition: librtcore.h:2437
rt_valuecount_t::count
uint32_t count
Definition: librtcore.h:2578
rt_valuecount_t::percent
double percent
Definition: librtcore.h:2579
rt_valuecount_t::value
double value
Definition: librtcore.h:2577
rt_valuecount_t
Definition: librtcore.h:2576
rtpg_summarystats_arg_t::cQ
double cQ
Definition: rtpg_statistics.c:500
rtpg_summarystats_arg_t::cM
double cM
Definition: rtpg_statistics.c:499
rtpg_summarystats_arg_t::band_index
int32_t band_index
Definition: rtpg_statistics.c:502
rtpg_summarystats_arg_t::cK
uint64_t cK
Definition: rtpg_statistics.c:498
rtpg_summarystats_arg_t::stats
rt_bandstats stats
Definition: rtpg_statistics.c:495
rtpg_summarystats_arg_t::sample
double sample
Definition: rtpg_statistics.c:504
rtpg_summarystats_arg_t::exclude_nodata_value
bool exclude_nodata_value
Definition: rtpg_statistics.c:503
rtpg_summarystats_arg_t
Definition: rtpg_statistics.c:494