PostGIS  2.4.9dev-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 Datum RASTER_histogramCoverage(PG_FUNCTION_ARGS);
58 
59 /* get quantiles */
60 Datum RASTER_quantile(PG_FUNCTION_ARGS);
61 Datum RASTER_quantileCoverage(PG_FUNCTION_ARGS);
62 
63 /* get counts of values */
64 Datum RASTER_valueCount(PG_FUNCTION_ARGS);
65 Datum RASTER_valueCountCoverage(PG_FUNCTION_ARGS);
66 
67 #define VALUES_LENGTH 6
68 
72 PG_FUNCTION_INFO_V1(RASTER_summaryStats);
73 Datum RASTER_summaryStats(PG_FUNCTION_ARGS)
74 {
75  rt_pgraster *pgraster = NULL;
76  rt_raster raster = NULL;
77  rt_band band = NULL;
78  int32_t bandindex = 1;
79  bool exclude_nodata_value = TRUE;
80  int num_bands = 0;
81  double sample = 0;
82  rt_bandstats stats = NULL;
83 
84  TupleDesc tupdesc;
85  Datum values[VALUES_LENGTH];
86  bool nulls[VALUES_LENGTH];
87  HeapTuple tuple;
88  Datum result;
89 
90  /* pgraster is null, return null */
91  if (PG_ARGISNULL(0))
92  PG_RETURN_NULL();
93  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
94 
95  raster = rt_raster_deserialize(pgraster, FALSE);
96  if (!raster) {
97  PG_FREE_IF_COPY(pgraster, 0);
98  elog(ERROR, "RASTER_summaryStats: Cannot deserialize raster");
99  PG_RETURN_NULL();
100  }
101 
102  /* band index is 1-based */
103  if (!PG_ARGISNULL(1))
104  bandindex = PG_GETARG_INT32(1);
105  num_bands = rt_raster_get_num_bands(raster);
106  if (bandindex < 1 || bandindex > num_bands) {
107  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
108  rt_raster_destroy(raster);
109  PG_FREE_IF_COPY(pgraster, 0);
110  PG_RETURN_NULL();
111  }
112 
113  /* exclude_nodata_value flag */
114  if (!PG_ARGISNULL(2))
115  exclude_nodata_value = PG_GETARG_BOOL(2);
116 
117  /* sample % */
118  if (!PG_ARGISNULL(3)) {
119  sample = PG_GETARG_FLOAT8(3);
120  if (sample < 0 || sample > 1) {
121  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
122  rt_raster_destroy(raster);
123  PG_FREE_IF_COPY(pgraster, 0);
124  PG_RETURN_NULL();
125  }
126  else if (FLT_EQ(sample, 0.0))
127  sample = 1;
128  }
129  else
130  sample = 1;
131 
132  /* get band */
133  band = rt_raster_get_band(raster, bandindex - 1);
134  if (!band) {
135  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
136  rt_raster_destroy(raster);
137  PG_FREE_IF_COPY(pgraster, 0);
138  PG_RETURN_NULL();
139  }
140 
141  /* we don't need the raw values, hence the zero parameter */
142  stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 0, NULL, NULL, NULL);
143  rt_band_destroy(band);
144  rt_raster_destroy(raster);
145  PG_FREE_IF_COPY(pgraster, 0);
146  if (NULL == stats) {
147  elog(NOTICE, "Cannot compute summary statistics for band at index %d. Returning NULL", bandindex);
148  PG_RETURN_NULL();
149  }
150 
151  /* Build a tuple descriptor for our result type */
152  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
153  ereport(ERROR, (
154  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
155  errmsg(
156  "function returning record called in context "
157  "that cannot accept type record"
158  )
159  ));
160  }
161 
162  BlessTupleDesc(tupdesc);
163 
164  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
165 
166  values[0] = Int64GetDatum(stats->count);
167  if (stats->count > 0) {
168  values[1] = Float8GetDatum(stats->sum);
169  values[2] = Float8GetDatum(stats->mean);
170  values[3] = Float8GetDatum(stats->stddev);
171  values[4] = Float8GetDatum(stats->min);
172  values[5] = Float8GetDatum(stats->max);
173  }
174  else {
175  nulls[1] = TRUE;
176  nulls[2] = TRUE;
177  nulls[3] = TRUE;
178  nulls[4] = TRUE;
179  nulls[5] = TRUE;
180  }
181 
182  /* build a tuple */
183  tuple = heap_form_tuple(tupdesc, values, nulls);
184 
185  /* make the tuple into a datum */
186  result = HeapTupleGetDatum(tuple);
187 
188  /* clean up */
189  pfree(stats);
190 
191  PG_RETURN_DATUM(result);
192 }
193 
197 PG_FUNCTION_INFO_V1(RASTER_summaryStatsCoverage);
198 Datum RASTER_summaryStatsCoverage(PG_FUNCTION_ARGS)
199 {
200  text *tablenametext = NULL;
201  char *tablename = NULL;
202  text *colnametext = NULL;
203  char *colname = NULL;
204  int32_t bandindex = 1;
205  bool exclude_nodata_value = TRUE;
206  double sample = 0;
207 
208  int len = 0;
209  char *sql = NULL;
210  int spi_result;
211  Portal portal;
212  TupleDesc tupdesc;
213  SPITupleTable *tuptable = NULL;
214  HeapTuple tuple;
215  Datum datum;
216  bool isNull = FALSE;
217 
218  rt_pgraster *pgraster = NULL;
219  rt_raster raster = NULL;
220  rt_band band = NULL;
221  int num_bands = 0;
222  uint64_t cK = 0;
223  double cM = 0;
224  double cQ = 0;
225  rt_bandstats stats = NULL;
226  rt_bandstats rtn = NULL;
227 
228  Datum values[VALUES_LENGTH];
229  bool nulls[VALUES_LENGTH];
230  Datum result;
231 
232  /* tablename is null, return null */
233  if (PG_ARGISNULL(0)) {
234  elog(NOTICE, "Table name must be provided");
235  PG_RETURN_NULL();
236  }
237  tablenametext = PG_GETARG_TEXT_P(0);
238  tablename = text_to_cstring(tablenametext);
239  if (!strlen(tablename)) {
240  elog(NOTICE, "Table name must be provided");
241  PG_RETURN_NULL();
242  }
243 
244  /* column name is null, return null */
245  if (PG_ARGISNULL(1)) {
246  elog(NOTICE, "Column name must be provided");
247  PG_RETURN_NULL();
248  }
249  colnametext = PG_GETARG_TEXT_P(1);
250  colname = text_to_cstring(colnametext);
251  if (!strlen(colname)) {
252  elog(NOTICE, "Column name must be provided");
253  PG_RETURN_NULL();
254  }
255 
256  /* band index is 1-based */
257  if (!PG_ARGISNULL(2))
258  bandindex = PG_GETARG_INT32(2);
259 
260  /* exclude_nodata_value flag */
261  if (!PG_ARGISNULL(3))
262  exclude_nodata_value = PG_GETARG_BOOL(3);
263 
264  /* sample % */
265  if (!PG_ARGISNULL(4)) {
266  sample = PG_GETARG_FLOAT8(4);
267  if (sample < 0 || sample > 1) {
268  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
269  rt_raster_destroy(raster);
270  PG_RETURN_NULL();
271  }
272  else if (FLT_EQ(sample, 0.0))
273  sample = 1;
274  }
275  else
276  sample = 1;
277 
278  /* iterate through rasters of coverage */
279  /* connect to database */
280  spi_result = SPI_connect();
281  if (spi_result != SPI_OK_CONNECT) {
282  pfree(sql);
283  elog(ERROR, "RASTER_summaryStatsCoverage: Cannot connect to database using SPI");
284  PG_RETURN_NULL();
285  }
286 
287  /* create sql */
288  len = sizeof(char) * (strlen("SELECT \"\" FROM \"\" WHERE \"\" IS NOT NULL") + (strlen(colname) * 2) + strlen(tablename) + 1);
289  sql = (char *) palloc(len);
290  if (NULL == sql) {
291  if (SPI_tuptable) SPI_freetuptable(tuptable);
292  SPI_finish();
293  elog(ERROR, "RASTER_summaryStatsCoverage: Cannot allocate memory for sql");
294  PG_RETURN_NULL();
295  }
296 
297  /* get cursor */
298  snprintf(sql, len, "SELECT \"%s\" FROM \"%s\" WHERE \"%s\" IS NOT NULL", colname, tablename, colname);
299  portal = SPI_cursor_open_with_args(
300  "coverage",
301  sql,
302  0, NULL,
303  NULL, NULL,
304  TRUE, 0
305  );
306  pfree(sql);
307 
308  /* process resultset */
309  SPI_cursor_fetch(portal, TRUE, 1);
310  while (SPI_processed == 1 && SPI_tuptable != NULL) {
311  tupdesc = SPI_tuptable->tupdesc;
312  tuptable = SPI_tuptable;
313  tuple = tuptable->vals[0];
314 
315  datum = SPI_getbinval(tuple, tupdesc, 1, &isNull);
316  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
317 
318  if (SPI_tuptable) SPI_freetuptable(tuptable);
319  SPI_cursor_close(portal);
320  SPI_finish();
321 
322  if (NULL != rtn) pfree(rtn);
323  elog(ERROR, "RASTER_summaryStatsCoverage: Cannot get raster of coverage");
324  PG_RETURN_NULL();
325  }
326  else if (isNull) {
327  SPI_cursor_fetch(portal, TRUE, 1);
328  continue;
329  }
330 
331  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(datum);
332 
333  raster = rt_raster_deserialize(pgraster, FALSE);
334  if (!raster) {
335 
336  if (SPI_tuptable) SPI_freetuptable(tuptable);
337  SPI_cursor_close(portal);
338  SPI_finish();
339 
340  if (NULL != rtn) pfree(rtn);
341  elog(ERROR, "RASTER_summaryStatsCoverage: Cannot deserialize raster");
342  PG_RETURN_NULL();
343  }
344 
345  /* inspect number of bands */
346  num_bands = rt_raster_get_num_bands(raster);
347  if (bandindex < 1 || bandindex > num_bands) {
348  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
349 
350  rt_raster_destroy(raster);
351 
352  if (SPI_tuptable) SPI_freetuptable(tuptable);
353  SPI_cursor_close(portal);
354  SPI_finish();
355 
356  if (NULL != rtn) pfree(rtn);
357  PG_RETURN_NULL();
358  }
359 
360  /* get band */
361  band = rt_raster_get_band(raster, bandindex - 1);
362  if (!band) {
363  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
364 
365  rt_raster_destroy(raster);
366 
367  if (SPI_tuptable) SPI_freetuptable(tuptable);
368  SPI_cursor_close(portal);
369  SPI_finish();
370 
371  if (NULL != rtn) pfree(rtn);
372  PG_RETURN_NULL();
373  }
374 
375  /* we don't need the raw values, hence the zero parameter */
376  stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 0, &cK, &cM, &cQ);
377 
378  rt_band_destroy(band);
379  rt_raster_destroy(raster);
380 
381  if (NULL == stats) {
382  elog(NOTICE, "Cannot compute summary statistics for band at index %d. Returning NULL", bandindex);
383 
384  if (SPI_tuptable) SPI_freetuptable(tuptable);
385  SPI_cursor_close(portal);
386  SPI_finish();
387 
388  if (NULL != rtn) pfree(rtn);
389  PG_RETURN_NULL();
390  }
391 
392  /* initialize rtn */
393  if (stats->count > 0) {
394  if (NULL == rtn) {
395  rtn = (rt_bandstats) SPI_palloc(sizeof(struct rt_bandstats_t));
396  if (NULL == rtn) {
397 
398  if (SPI_tuptable) SPI_freetuptable(tuptable);
399  SPI_cursor_close(portal);
400  SPI_finish();
401 
402  elog(ERROR, "RASTER_summaryStatsCoverage: Cannot allocate memory for summary stats of coverage");
403  PG_RETURN_NULL();
404  }
405 
406  rtn->sample = stats->sample;
407  rtn->count = stats->count;
408  rtn->min = stats->min;
409  rtn->max = stats->max;
410  rtn->sum = stats->sum;
411  rtn->mean = stats->mean;
412  rtn->stddev = -1;
413 
414  rtn->values = NULL;
415  rtn->sorted = 0;
416  }
417  else {
418  rtn->count += stats->count;
419  rtn->sum += stats->sum;
420 
421  if (stats->min < rtn->min)
422  rtn->min = stats->min;
423  if (stats->max > rtn->max)
424  rtn->max = stats->max;
425  }
426  }
427 
428  pfree(stats);
429 
430  /* next record */
431  SPI_cursor_fetch(portal, TRUE, 1);
432  }
433 
434  if (SPI_tuptable) SPI_freetuptable(tuptable);
435  SPI_cursor_close(portal);
436  SPI_finish();
437 
438  if (NULL == rtn) {
439  elog(ERROR, "RASTER_summaryStatsCoverage: Cannot compute coverage summary stats");
440  PG_RETURN_NULL();
441  }
442 
443  /* coverage mean and deviation */
444  rtn->mean = rtn->sum / rtn->count;
445  /* sample deviation */
446  if (rtn->sample > 0 && rtn->sample < 1)
447  rtn->stddev = sqrt(cQ / (rtn->count - 1));
448  /* standard deviation */
449  else
450  rtn->stddev = sqrt(cQ / rtn->count);
451 
452  /* Build a tuple descriptor for our result type */
453  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
454  ereport(ERROR, (
455  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
456  errmsg(
457  "function returning record called in context "
458  "that cannot accept type record"
459  )
460  ));
461  }
462 
463  BlessTupleDesc(tupdesc);
464 
465  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
466 
467  values[0] = Int64GetDatum(rtn->count);
468  if (rtn->count > 0) {
469  values[1] = Float8GetDatum(rtn->sum);
470  values[2] = Float8GetDatum(rtn->mean);
471  values[3] = Float8GetDatum(rtn->stddev);
472  values[4] = Float8GetDatum(rtn->min);
473  values[5] = Float8GetDatum(rtn->max);
474  }
475  else {
476  nulls[1] = TRUE;
477  nulls[2] = TRUE;
478  nulls[3] = TRUE;
479  nulls[4] = TRUE;
480  nulls[5] = TRUE;
481  }
482 
483  /* build a tuple */
484  tuple = heap_form_tuple(tupdesc, values, nulls);
485 
486  /* make the tuple into a datum */
487  result = HeapTupleGetDatum(tuple);
488 
489  /* clean up */
490  pfree(rtn);
491 
492  PG_RETURN_DATUM(result);
493 }
494 
495 /* ---------------------------------------------------------------- */
496 /* Aggregate ST_SummaryStats */
497 /* ---------------------------------------------------------------- */
498 
499 typedef struct rtpg_summarystats_arg_t *rtpg_summarystats_arg;
500 struct rtpg_summarystats_arg_t {
501  rt_bandstats stats;
502 
503  /* coefficients for one-pass standard deviation */
504  uint64_t cK;
505  double cM;
506  double cQ;
507 
508  int32_t band_index; /* one-based */
509  bool exclude_nodata_value;
510  double sample; /* value between 0 and 1 */
511 };
512 
513 static void
514 rtpg_summarystats_arg_destroy(rtpg_summarystats_arg arg) {
515  if (arg->stats != NULL)
516  pfree(arg->stats);
517 
518  pfree(arg);
519 }
520 
521 static rtpg_summarystats_arg
522 rtpg_summarystats_arg_init() {
523  rtpg_summarystats_arg arg = NULL;
524 
525  arg = palloc(sizeof(struct rtpg_summarystats_arg_t));
526  if (arg == NULL) {
527  elog(
528  ERROR,
529  "rtpg_summarystats_arg_init: Cannot allocate memory for function arguments"
530  );
531  return NULL;
532  }
533 
534  arg->stats = (rt_bandstats) palloc(sizeof(struct rt_bandstats_t));
535  if (arg->stats == NULL) {
536  rtpg_summarystats_arg_destroy(arg);
537  elog(
538  ERROR,
539  "rtpg_summarystats_arg_init: Cannot allocate memory for stats function argument"
540  );
541  return NULL;
542  }
543 
544  arg->stats->sample = 0;
545  arg->stats->count = 0;
546  arg->stats->min = 0;
547  arg->stats->max = 0;
548  arg->stats->sum = 0;
549  arg->stats->mean = 0;
550  arg->stats->stddev = -1;
551  arg->stats->values = NULL;
552  arg->stats->sorted = 0;
553 
554  arg->cK = 0;
555  arg->cM = 0;
556  arg->cQ = 0;
557 
558  arg->band_index = 1;
559  arg->exclude_nodata_value = TRUE;
560  arg->sample = 1;
561 
562  return arg;
563 }
564 
565 PG_FUNCTION_INFO_V1(RASTER_summaryStats_transfn);
566 Datum RASTER_summaryStats_transfn(PG_FUNCTION_ARGS)
567 {
568  MemoryContext aggcontext;
569  MemoryContext oldcontext;
570  rtpg_summarystats_arg state = NULL;
571  bool skiparg = FALSE;
572 
573  int i = 0;
574 
575  rt_pgraster *pgraster = NULL;
576  rt_raster raster = NULL;
577  rt_band band = NULL;
578  int num_bands = 0;
579  rt_bandstats stats = NULL;
580 
581  POSTGIS_RT_DEBUG(3, "Starting...");
582 
583  /* cannot be called directly as this is exclusive aggregate function */
584  if (!AggCheckCallContext(fcinfo, &aggcontext)) {
585  elog(
586  ERROR,
587  "RASTER_summaryStats_transfn: Cannot be called in a non-aggregate context"
588  );
589  PG_RETURN_NULL();
590  }
591 
592  /* switch to aggcontext */
593  oldcontext = MemoryContextSwitchTo(aggcontext);
594 
595  if (PG_ARGISNULL(0)) {
596  POSTGIS_RT_DEBUG(3, "Creating state variable");
597 
598  state = rtpg_summarystats_arg_init();
599  if (state == NULL) {
600  MemoryContextSwitchTo(oldcontext);
601  elog(
602  ERROR,
603  "RASTER_summaryStats_transfn: Cannot allocate memory for state variable"
604  );
605  PG_RETURN_NULL();
606  }
607 
608  skiparg = FALSE;
609  }
610  else {
611  POSTGIS_RT_DEBUG(3, "State variable already exists");
612  state = (rtpg_summarystats_arg) PG_GETARG_POINTER(0);
613  skiparg = TRUE;
614  }
615 
616  /* raster arg is NOT NULL */
617  if (!PG_ARGISNULL(1)) {
618  /* deserialize raster */
619  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
620 
621  /* Get raster object */
622  raster = rt_raster_deserialize(pgraster, FALSE);
623  if (raster == NULL) {
624 
625  rtpg_summarystats_arg_destroy(state);
626  PG_FREE_IF_COPY(pgraster, 1);
627 
628  MemoryContextSwitchTo(oldcontext);
629  elog(ERROR, "RASTER_summaryStats_transfn: Cannot deserialize raster");
630  PG_RETURN_NULL();
631  }
632  }
633 
634  do {
635  Oid calltype;
636  int nargs = 0;
637 
638  if (skiparg)
639  break;
640 
641  /* 4 or 5 total possible args */
642  nargs = PG_NARGS();
643  POSTGIS_RT_DEBUGF(4, "nargs = %d", nargs);
644 
645  for (i = 2; i < nargs; i++) {
646  if (PG_ARGISNULL(i))
647  continue;
648 
649  calltype = get_fn_expr_argtype(fcinfo->flinfo, i);
650 
651  /* band index */
652  if (
653  (calltype == INT2OID || calltype == INT4OID) &&
654  i == 2
655  ) {
656  if (calltype == INT2OID)
657  state->band_index = PG_GETARG_INT16(i);
658  else
659  state->band_index = PG_GETARG_INT32(i);
660 
661  /* basic check, > 0 */
662  if (state->band_index < 1) {
663 
664  rtpg_summarystats_arg_destroy(state);
665  if (raster != NULL) {
666  rt_raster_destroy(raster);
667  PG_FREE_IF_COPY(pgraster, 1);
668  }
669 
670  MemoryContextSwitchTo(oldcontext);
671  elog(
672  ERROR,
673  "RASTER_summaryStats_transfn: Invalid band index (must use 1-based). Returning NULL"
674  );
675  PG_RETURN_NULL();
676  }
677  }
678  /* exclude_nodata_value */
679  else if (
680  calltype == BOOLOID && (
681  i == 2 || i == 3
682  )
683  ) {
684  state->exclude_nodata_value = PG_GETARG_BOOL(i);
685  }
686  /* sample rate */
687  else if (
688  (calltype == FLOAT4OID || calltype == FLOAT8OID) &&
689  (i == 3 || i == 4)
690  ) {
691  if (calltype == FLOAT4OID)
692  state->sample = PG_GETARG_FLOAT4(i);
693  else
694  state->sample = PG_GETARG_FLOAT8(i);
695 
696  /* basic check, 0 <= sample <= 1 */
697  if (state->sample < 0. || state->sample > 1.) {
698 
699  rtpg_summarystats_arg_destroy(state);
700  if (raster != NULL) {
701  rt_raster_destroy(raster);
702  PG_FREE_IF_COPY(pgraster, 1);
703  }
704 
705  MemoryContextSwitchTo(oldcontext);
706  elog(
707  ERROR,
708  "Invalid sample percentage (must be between 0 and 1). Returning NULL"
709  );
710 
711  PG_RETURN_NULL();
712  }
713  else if (FLT_EQ(state->sample, 0.0))
714  state->sample = 1;
715  }
716  /* unknown arg */
717  else {
718  rtpg_summarystats_arg_destroy(state);
719  if (raster != NULL) {
720  rt_raster_destroy(raster);
721  PG_FREE_IF_COPY(pgraster, 1);
722  }
723 
724  MemoryContextSwitchTo(oldcontext);
725  elog(
726  ERROR,
727  "RASTER_summaryStats_transfn: Unknown function parameter at index %d",
728  i
729  );
730  PG_RETURN_NULL();
731  }
732  }
733  }
734  while (0);
735 
736  /* null raster, return */
737  if (PG_ARGISNULL(1)) {
738  POSTGIS_RT_DEBUG(4, "NULL raster so processing required");
739  MemoryContextSwitchTo(oldcontext);
740  PG_RETURN_POINTER(state);
741  }
742 
743  /* inspect number of bands */
744  num_bands = rt_raster_get_num_bands(raster);
745  if (state->band_index > num_bands) {
746  elog(
747  NOTICE,
748  "Raster does not have band at index %d. Skipping raster",
749  state->band_index
750  );
751 
752  rt_raster_destroy(raster);
753  PG_FREE_IF_COPY(pgraster, 1);
754 
755  MemoryContextSwitchTo(oldcontext);
756  PG_RETURN_POINTER(state);
757  }
758 
759  /* get band */
760  band = rt_raster_get_band(raster, state->band_index - 1);
761  if (!band) {
762  elog(
763  NOTICE, "Cannot find band at index %d. Skipping raster",
764  state->band_index
765  );
766 
767  rt_raster_destroy(raster);
768  PG_FREE_IF_COPY(pgraster, 1);
769 
770  MemoryContextSwitchTo(oldcontext);
771  PG_RETURN_POINTER(state);
772  }
773 
774  /* we don't need the raw values, hence the zero parameter */
775  stats = rt_band_get_summary_stats(
776  band, (int) state->exclude_nodata_value,
777  state->sample, 0,
778  &(state->cK), &(state->cM), &(state->cQ)
779  );
780 
781  rt_band_destroy(band);
782  rt_raster_destroy(raster);
783  PG_FREE_IF_COPY(pgraster, 1);
784 
785  if (NULL == stats) {
786  elog(
787  NOTICE,
788  "Cannot compute summary statistics for band at index %d. Returning NULL",
789  state->band_index
790  );
791 
792  rtpg_summarystats_arg_destroy(state);
793 
794  MemoryContextSwitchTo(oldcontext);
795  PG_RETURN_NULL();
796  }
797 
798  if (stats->count > 0) {
799  if (state->stats->count < 1) {
800  state->stats->sample = stats->sample;
801  state->stats->count = stats->count;
802  state->stats->min = stats->min;
803  state->stats->max = stats->max;
804  state->stats->sum = stats->sum;
805  state->stats->mean = stats->mean;
806  state->stats->stddev = -1;
807  }
808  else {
809  state->stats->count += stats->count;
810  state->stats->sum += stats->sum;
811 
812  if (stats->min < state->stats->min)
813  state->stats->min = stats->min;
814  if (stats->max > state->stats->max)
815  state->stats->max = stats->max;
816  }
817  }
818 
819  pfree(stats);
820 
821  /* switch back to local context */
822  MemoryContextSwitchTo(oldcontext);
823 
824  POSTGIS_RT_DEBUG(3, "Finished");
825 
826  PG_RETURN_POINTER(state);
827 }
828 
829 PG_FUNCTION_INFO_V1(RASTER_summaryStats_finalfn);
830 Datum RASTER_summaryStats_finalfn(PG_FUNCTION_ARGS)
831 {
832  rtpg_summarystats_arg state = NULL;
833 
834  TupleDesc tupdesc;
835  HeapTuple tuple;
836  Datum values[VALUES_LENGTH];
837  bool nulls[VALUES_LENGTH];
838  Datum result;
839 
840  POSTGIS_RT_DEBUG(3, "Starting...");
841 
842  /* cannot be called directly as this is exclusive aggregate function */
843  if (!AggCheckCallContext(fcinfo, NULL)) {
844  elog(ERROR, "RASTER_summaryStats_finalfn: Cannot be called in a non-aggregate context");
845  PG_RETURN_NULL();
846  }
847 
848  /* NULL, return null */
849  if (PG_ARGISNULL(0))
850  PG_RETURN_NULL();
851 
852  state = (rtpg_summarystats_arg) PG_GETARG_POINTER(0);
853 
854  if (NULL == state) {
855  elog(ERROR, "RASTER_summaryStats_finalfn: Cannot compute coverage summary stats");
856  PG_RETURN_NULL();
857  }
858 
859  /* coverage mean and deviation */
860  if (state->stats->count > 0) {
861  state->stats->mean = state->stats->sum / state->stats->count;
862  /* sample deviation */
863  if (state->stats->sample > 0 && state->stats->sample < 1)
864  state->stats->stddev = sqrt(state->cQ / (state->stats->count - 1));
865  /* standard deviation */
866  else
867  state->stats->stddev = sqrt(state->cQ / state->stats->count);
868  }
869 
870  /* Build a tuple descriptor for our result type */
871  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
872  rtpg_summarystats_arg_destroy(state);
873  ereport(ERROR, (
874  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
875  errmsg(
876  "function returning record called in context "
877  "that cannot accept type record"
878  )
879  ));
880  }
881 
882  BlessTupleDesc(tupdesc);
883 
884  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
885 
886  values[0] = Int64GetDatum(state->stats->count);
887  if (state->stats->count > 0) {
888  values[1] = Float8GetDatum(state->stats->sum);
889  values[2] = Float8GetDatum(state->stats->mean);
890  values[3] = Float8GetDatum(state->stats->stddev);
891  values[4] = Float8GetDatum(state->stats->min);
892  values[5] = Float8GetDatum(state->stats->max);
893  }
894  else {
895  nulls[1] = TRUE;
896  nulls[2] = TRUE;
897  nulls[3] = TRUE;
898  nulls[4] = TRUE;
899  nulls[5] = TRUE;
900  }
901 
902  /* build a tuple */
903  tuple = heap_form_tuple(tupdesc, values, nulls);
904 
905  /* make the tuple into a datum */
906  result = HeapTupleGetDatum(tuple);
907 
908  /* clean up */
909  rtpg_summarystats_arg_destroy(state);
910 
911  PG_RETURN_DATUM(result);
912 }
913 
914 #undef VALUES_LENGTH
915 #define VALUES_LENGTH 4
916 
920 PG_FUNCTION_INFO_V1(RASTER_histogram);
921 Datum RASTER_histogram(PG_FUNCTION_ARGS)
922 {
923  FuncCallContext *funcctx;
924  TupleDesc tupdesc;
925 
926  int i;
927  rt_histogram hist;
928  rt_histogram hist2;
929  int call_cntr;
930  int max_calls;
931 
932  /* first call of function */
933  if (SRF_IS_FIRSTCALL()) {
934  MemoryContext oldcontext;
935 
936  rt_pgraster *pgraster = NULL;
937  rt_raster raster = NULL;
938  rt_band band = NULL;
939  int32_t bandindex = 1;
940  int num_bands = 0;
941  bool exclude_nodata_value = TRUE;
942  double sample = 0;
943  uint32_t bin_count = 0;
944  double *bin_width = NULL;
945  uint32_t bin_width_count = 0;
946  double width = 0;
947  bool right = FALSE;
948  double min = 0;
949  double max = 0;
950  rt_bandstats stats = NULL;
951  uint32_t count;
952 
953  int j;
954  int n;
955 
956  ArrayType *array;
957  Oid etype;
958  Datum *e;
959  bool *nulls;
960  int16 typlen;
961  bool typbyval;
962  char typalign;
963 
964  POSTGIS_RT_DEBUG(3, "RASTER_histogram: Starting");
965 
966  /* create a function context for cross-call persistence */
967  funcctx = SRF_FIRSTCALL_INIT();
968 
969  /* switch to memory context appropriate for multiple function calls */
970  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
971 
972  /* pgraster is null, return nothing */
973  if (PG_ARGISNULL(0)) {
974  MemoryContextSwitchTo(oldcontext);
975  SRF_RETURN_DONE(funcctx);
976  }
977  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
978 
979  raster = rt_raster_deserialize(pgraster, FALSE);
980  if (!raster) {
981  PG_FREE_IF_COPY(pgraster, 0);
982  MemoryContextSwitchTo(oldcontext);
983  elog(ERROR, "RASTER_histogram: Cannot deserialize raster");
984  SRF_RETURN_DONE(funcctx);
985  }
986 
987  /* band index is 1-based */
988  if (!PG_ARGISNULL(1))
989  bandindex = PG_GETARG_INT32(1);
990  num_bands = rt_raster_get_num_bands(raster);
991  if (bandindex < 1 || bandindex > num_bands) {
992  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
993  rt_raster_destroy(raster);
994  PG_FREE_IF_COPY(pgraster, 0);
995  MemoryContextSwitchTo(oldcontext);
996  SRF_RETURN_DONE(funcctx);
997  }
998 
999  /* exclude_nodata_value flag */
1000  if (!PG_ARGISNULL(2))
1001  exclude_nodata_value = PG_GETARG_BOOL(2);
1002 
1003  /* sample % */
1004  if (!PG_ARGISNULL(3)) {
1005  sample = PG_GETARG_FLOAT8(3);
1006  if (sample < 0 || sample > 1) {
1007  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
1008  rt_raster_destroy(raster);
1009  PG_FREE_IF_COPY(pgraster, 0);
1010  MemoryContextSwitchTo(oldcontext);
1011  SRF_RETURN_DONE(funcctx);
1012  }
1013  else if (FLT_EQ(sample, 0.0))
1014  sample = 1;
1015  }
1016  else
1017  sample = 1;
1018 
1019  /* bin_count */
1020  if (!PG_ARGISNULL(4)) {
1021  bin_count = PG_GETARG_INT32(4);
1022  if (bin_count < 1) bin_count = 0;
1023  }
1024 
1025  /* bin_width */
1026  if (!PG_ARGISNULL(5)) {
1027  array = PG_GETARG_ARRAYTYPE_P(5);
1028  etype = ARR_ELEMTYPE(array);
1029  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1030 
1031  switch (etype) {
1032  case FLOAT4OID:
1033  case FLOAT8OID:
1034  break;
1035  default:
1036  rt_raster_destroy(raster);
1037  PG_FREE_IF_COPY(pgraster, 0);
1038  MemoryContextSwitchTo(oldcontext);
1039  elog(ERROR, "RASTER_histogram: Invalid data type for width");
1040  SRF_RETURN_DONE(funcctx);
1041  break;
1042  }
1043 
1044  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1045  &nulls, &n);
1046 
1047  bin_width = palloc(sizeof(double) * n);
1048  for (i = 0, j = 0; i < n; i++) {
1049  if (nulls[i]) continue;
1050 
1051  switch (etype) {
1052  case FLOAT4OID:
1053  width = (double) DatumGetFloat4(e[i]);
1054  break;
1055  case FLOAT8OID:
1056  width = (double) DatumGetFloat8(e[i]);
1057  break;
1058  }
1059 
1060  if (width < 0 || FLT_EQ(width, 0.0)) {
1061  elog(NOTICE, "Invalid value for width (must be greater than 0). Returning NULL");
1062  pfree(bin_width);
1063  rt_raster_destroy(raster);
1064  PG_FREE_IF_COPY(pgraster, 0);
1065  MemoryContextSwitchTo(oldcontext);
1066  SRF_RETURN_DONE(funcctx);
1067  }
1068 
1069  bin_width[j] = width;
1070  POSTGIS_RT_DEBUGF(5, "bin_width[%d] = %f", j, bin_width[j]);
1071  j++;
1072  }
1073  bin_width_count = j;
1074 
1075  if (j < 1) {
1076  pfree(bin_width);
1077  bin_width = NULL;
1078  }
1079  }
1080 
1081  /* right */
1082  if (!PG_ARGISNULL(6))
1083  right = PG_GETARG_BOOL(6);
1084 
1085  /* min */
1086  if (!PG_ARGISNULL(7)) min = PG_GETARG_FLOAT8(7);
1087 
1088  /* max */
1089  if (!PG_ARGISNULL(8)) max = PG_GETARG_FLOAT8(8);
1090 
1091  /* get band */
1092  band = rt_raster_get_band(raster, bandindex - 1);
1093  if (!band) {
1094  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
1095  rt_raster_destroy(raster);
1096  PG_FREE_IF_COPY(pgraster, 0);
1097  MemoryContextSwitchTo(oldcontext);
1098  SRF_RETURN_DONE(funcctx);
1099  }
1100 
1101  /* get stats */
1102  stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 1, NULL, NULL, NULL);
1103  rt_band_destroy(band);
1104  rt_raster_destroy(raster);
1105  PG_FREE_IF_COPY(pgraster, 0);
1106  if (NULL == stats || NULL == stats->values) {
1107  elog(NOTICE, "Cannot compute summary statistics for band at index %d", bandindex);
1108  MemoryContextSwitchTo(oldcontext);
1109  SRF_RETURN_DONE(funcctx);
1110  }
1111  else if (stats->count < 1) {
1112  elog(NOTICE, "Cannot compute histogram for band at index %d as the band has no values", bandindex);
1113  MemoryContextSwitchTo(oldcontext);
1114  SRF_RETURN_DONE(funcctx);
1115  }
1116 
1117  /* get histogram */
1118  hist = rt_band_get_histogram(stats, bin_count, bin_width, bin_width_count, right, min, max, &count);
1119  if (bin_width_count) pfree(bin_width);
1120  pfree(stats);
1121  if (NULL == hist || !count) {
1122  elog(NOTICE, "Cannot compute histogram for band at index %d", bandindex);
1123  MemoryContextSwitchTo(oldcontext);
1124  SRF_RETURN_DONE(funcctx);
1125  }
1126 
1127  POSTGIS_RT_DEBUGF(3, "%d bins returned", count);
1128 
1129  /* Store needed information */
1130  funcctx->user_fctx = hist;
1131 
1132  /* total number of tuples to be returned */
1133  funcctx->max_calls = count;
1134 
1135  /* Build a tuple descriptor for our result type */
1136  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
1137  ereport(ERROR, (
1138  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1139  errmsg(
1140  "function returning record called in context "
1141  "that cannot accept type record"
1142  )
1143  ));
1144  }
1145 
1146  BlessTupleDesc(tupdesc);
1147  funcctx->tuple_desc = tupdesc;
1148 
1149  MemoryContextSwitchTo(oldcontext);
1150  }
1151 
1152  /* stuff done on every call of the function */
1153  funcctx = SRF_PERCALL_SETUP();
1154 
1155  call_cntr = funcctx->call_cntr;
1156  max_calls = funcctx->max_calls;
1157  tupdesc = funcctx->tuple_desc;
1158  hist2 = funcctx->user_fctx;
1159 
1160  /* do when there is more left to send */
1161  if (call_cntr < max_calls) {
1162  Datum values[VALUES_LENGTH];
1163  bool nulls[VALUES_LENGTH];
1164  HeapTuple tuple;
1165  Datum result;
1166 
1167  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
1168 
1169  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
1170 
1171  values[0] = Float8GetDatum(hist2[call_cntr].min);
1172  values[1] = Float8GetDatum(hist2[call_cntr].max);
1173  values[2] = Int64GetDatum(hist2[call_cntr].count);
1174  values[3] = Float8GetDatum(hist2[call_cntr].percent);
1175 
1176  /* build a tuple */
1177  tuple = heap_form_tuple(tupdesc, values, nulls);
1178 
1179  /* make the tuple into a datum */
1180  result = HeapTupleGetDatum(tuple);
1181 
1182  SRF_RETURN_NEXT(funcctx, result);
1183  }
1184  /* do when there is no more left */
1185  else {
1186  pfree(hist2);
1187  SRF_RETURN_DONE(funcctx);
1188  }
1189 }
1190 
1194 PG_FUNCTION_INFO_V1(RASTER_histogramCoverage);
1195 Datum RASTER_histogramCoverage(PG_FUNCTION_ARGS)
1196 {
1197  FuncCallContext *funcctx;
1198  TupleDesc tupdesc;
1199 
1200  int i;
1201  rt_histogram covhist = NULL;
1202  rt_histogram covhist2;
1203  int call_cntr;
1204  int max_calls;
1205 
1206  POSTGIS_RT_DEBUG(3, "RASTER_histogramCoverage: Starting");
1207 
1208  /* first call of function */
1209  if (SRF_IS_FIRSTCALL()) {
1210  MemoryContext oldcontext;
1211 
1212  text *tablenametext = NULL;
1213  char *tablename = NULL;
1214  text *colnametext = NULL;
1215  char *colname = NULL;
1216  int32_t bandindex = 1;
1217  bool exclude_nodata_value = TRUE;
1218  double sample = 0;
1219  uint32_t bin_count = 0;
1220  double *bin_width = NULL;
1221  uint32_t bin_width_count = 0;
1222  double width = 0;
1223  bool right = FALSE;
1224  uint32_t count;
1225 
1226  int len = 0;
1227  char *sql = NULL;
1228  char *tmp = NULL;
1229  double min = 0;
1230  double max = 0;
1231  int spi_result;
1232  Portal portal;
1233  SPITupleTable *tuptable = NULL;
1234  HeapTuple tuple;
1235  Datum datum;
1236  bool isNull = FALSE;
1237 
1238  rt_pgraster *pgraster = NULL;
1239  rt_raster raster = NULL;
1240  rt_band band = NULL;
1241  int num_bands = 0;
1242  rt_bandstats stats = NULL;
1243  rt_histogram hist;
1244  uint64_t sum = 0;
1245 
1246  int j;
1247  int n;
1248 
1249  ArrayType *array;
1250  Oid etype;
1251  Datum *e;
1252  bool *nulls;
1253  int16 typlen;
1254  bool typbyval;
1255  char typalign;
1256 
1257  POSTGIS_RT_DEBUG(3, "RASTER_histogramCoverage: first call of function");
1258 
1259  /* create a function context for cross-call persistence */
1260  funcctx = SRF_FIRSTCALL_INIT();
1261 
1262  /* switch to memory context appropriate for multiple function calls */
1263  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
1264 
1265  /* tablename is null, return null */
1266  if (PG_ARGISNULL(0)) {
1267  elog(NOTICE, "Table name must be provided");
1268  MemoryContextSwitchTo(oldcontext);
1269  SRF_RETURN_DONE(funcctx);
1270  }
1271  tablenametext = PG_GETARG_TEXT_P(0);
1272  tablename = text_to_cstring(tablenametext);
1273  if (!strlen(tablename)) {
1274  elog(NOTICE, "Table name must be provided");
1275  MemoryContextSwitchTo(oldcontext);
1276  SRF_RETURN_DONE(funcctx);
1277  }
1278  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: tablename = %s", tablename);
1279 
1280  /* column name is null, return null */
1281  if (PG_ARGISNULL(1)) {
1282  elog(NOTICE, "Column name must be provided");
1283  MemoryContextSwitchTo(oldcontext);
1284  SRF_RETURN_DONE(funcctx);
1285  }
1286  colnametext = PG_GETARG_TEXT_P(1);
1287  colname = text_to_cstring(colnametext);
1288  if (!strlen(colname)) {
1289  elog(NOTICE, "Column name must be provided");
1290  MemoryContextSwitchTo(oldcontext);
1291  SRF_RETURN_DONE(funcctx);
1292  }
1293  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: colname = %s", colname);
1294 
1295  /* band index is 1-based */
1296  if (!PG_ARGISNULL(2))
1297  bandindex = PG_GETARG_INT32(2);
1298 
1299  /* exclude_nodata_value flag */
1300  if (!PG_ARGISNULL(3))
1301  exclude_nodata_value = PG_GETARG_BOOL(3);
1302 
1303  /* sample % */
1304  if (!PG_ARGISNULL(4)) {
1305  sample = PG_GETARG_FLOAT8(4);
1306  if (sample < 0 || sample > 1) {
1307  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
1308  MemoryContextSwitchTo(oldcontext);
1309  SRF_RETURN_DONE(funcctx);
1310  }
1311  else if (FLT_EQ(sample, 0.0))
1312  sample = 1;
1313  }
1314  else
1315  sample = 1;
1316 
1317  /* bin_count */
1318  if (!PG_ARGISNULL(5)) {
1319  bin_count = PG_GETARG_INT32(5);
1320  if (bin_count < 1) bin_count = 0;
1321  }
1322 
1323  /* bin_width */
1324  if (!PG_ARGISNULL(6)) {
1325  array = PG_GETARG_ARRAYTYPE_P(6);
1326  etype = ARR_ELEMTYPE(array);
1327  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1328 
1329  switch (etype) {
1330  case FLOAT4OID:
1331  case FLOAT8OID:
1332  break;
1333  default:
1334  MemoryContextSwitchTo(oldcontext);
1335  elog(ERROR, "RASTER_histogramCoverage: Invalid data type for width");
1336  SRF_RETURN_DONE(funcctx);
1337  break;
1338  }
1339 
1340  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1341  &nulls, &n);
1342 
1343  bin_width = palloc(sizeof(double) * n);
1344  for (i = 0, j = 0; i < n; i++) {
1345  if (nulls[i]) continue;
1346 
1347  switch (etype) {
1348  case FLOAT4OID:
1349  width = (double) DatumGetFloat4(e[i]);
1350  break;
1351  case FLOAT8OID:
1352  width = (double) DatumGetFloat8(e[i]);
1353  break;
1354  }
1355 
1356  if (width < 0 || FLT_EQ(width, 0.0)) {
1357  elog(NOTICE, "Invalid value for width (must be greater than 0). Returning NULL");
1358  pfree(bin_width);
1359  MemoryContextSwitchTo(oldcontext);
1360  SRF_RETURN_DONE(funcctx);
1361  }
1362 
1363  bin_width[j] = width;
1364  POSTGIS_RT_DEBUGF(5, "bin_width[%d] = %f", j, bin_width[j]);
1365  j++;
1366  }
1367  bin_width_count = j;
1368 
1369  if (j < 1) {
1370  pfree(bin_width);
1371  bin_width = NULL;
1372  }
1373  }
1374 
1375  /* right */
1376  if (!PG_ARGISNULL(7))
1377  right = PG_GETARG_BOOL(7);
1378 
1379  /* connect to database */
1380  spi_result = SPI_connect();
1381  if (spi_result != SPI_OK_CONNECT) {
1382 
1383  if (bin_width_count) pfree(bin_width);
1384 
1385  MemoryContextSwitchTo(oldcontext);
1386  elog(ERROR, "RASTER_histogramCoverage: Cannot connect to database using SPI");
1387  SRF_RETURN_DONE(funcctx);
1388  }
1389 
1390  /* coverage stats */
1391  len = sizeof(char) * (strlen("SELECT min, max FROM _st_summarystats('','',,::boolean,)") + strlen(tablename) + strlen(colname) + (MAX_INT_CHARLEN * 2) + MAX_DBL_CHARLEN + 1);
1392  sql = (char *) palloc(len);
1393  if (NULL == sql) {
1394 
1395  if (SPI_tuptable) SPI_freetuptable(tuptable);
1396  SPI_finish();
1397 
1398  if (bin_width_count) pfree(bin_width);
1399 
1400  MemoryContextSwitchTo(oldcontext);
1401  elog(ERROR, "RASTER_histogramCoverage: Cannot allocate memory for sql");
1402  SRF_RETURN_DONE(funcctx);
1403  }
1404 
1405  /* get stats */
1406  snprintf(sql, len, "SELECT min, max FROM _st_summarystats('%s','%s',%d,%d::boolean,%f)", tablename, colname, bandindex, (exclude_nodata_value ? 1 : 0), sample);
1407  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: %s", sql);
1408  spi_result = SPI_execute(sql, TRUE, 0);
1409  pfree(sql);
1410  if (spi_result != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
1411 
1412  if (SPI_tuptable) SPI_freetuptable(tuptable);
1413  SPI_finish();
1414 
1415  if (bin_width_count) pfree(bin_width);
1416 
1417  MemoryContextSwitchTo(oldcontext);
1418  elog(ERROR, "RASTER_histogramCoverage: Cannot get summary stats of coverage");
1419  SRF_RETURN_DONE(funcctx);
1420  }
1421 
1422  tupdesc = SPI_tuptable->tupdesc;
1423  tuptable = SPI_tuptable;
1424  tuple = tuptable->vals[0];
1425 
1426  tmp = SPI_getvalue(tuple, tupdesc, 1);
1427  if (NULL == tmp || !strlen(tmp)) {
1428 
1429  if (SPI_tuptable) SPI_freetuptable(tuptable);
1430  SPI_finish();
1431 
1432  if (bin_width_count) pfree(bin_width);
1433 
1434  MemoryContextSwitchTo(oldcontext);
1435  elog(ERROR, "RASTER_histogramCoverage: Cannot get summary stats of coverage");
1436  SRF_RETURN_DONE(funcctx);
1437  }
1438  min = strtod(tmp, NULL);
1439  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: min = %f", min);
1440  pfree(tmp);
1441 
1442  tmp = SPI_getvalue(tuple, tupdesc, 2);
1443  if (NULL == tmp || !strlen(tmp)) {
1444 
1445  if (SPI_tuptable) SPI_freetuptable(tuptable);
1446  SPI_finish();
1447 
1448  if (bin_width_count) pfree(bin_width);
1449 
1450  MemoryContextSwitchTo(oldcontext);
1451  elog(ERROR, "RASTER_histogramCoverage: Cannot get summary stats of coverage");
1452  SRF_RETURN_DONE(funcctx);
1453  }
1454  max = strtod(tmp, NULL);
1455  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: max = %f", max);
1456  pfree(tmp);
1457 
1458  /* iterate through rasters of coverage */
1459  /* create sql */
1460  len = sizeof(char) * (strlen("SELECT \"\" FROM \"\" WHERE \"\" IS NOT NULL") + (strlen(colname) * 2) + strlen(tablename) + 1);
1461  sql = (char *) palloc(len);
1462  if (NULL == sql) {
1463 
1464  if (SPI_tuptable) SPI_freetuptable(tuptable);
1465  SPI_finish();
1466 
1467  if (bin_width_count) pfree(bin_width);
1468 
1469  MemoryContextSwitchTo(oldcontext);
1470  elog(ERROR, "RASTER_histogramCoverage: Cannot allocate memory for sql");
1471  SRF_RETURN_DONE(funcctx);
1472  }
1473 
1474  /* get cursor */
1475  snprintf(sql, len, "SELECT \"%s\" FROM \"%s\" WHERE \"%s\" IS NOT NULL", colname, tablename, colname);
1476  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: %s", sql);
1477  portal = SPI_cursor_open_with_args(
1478  "coverage",
1479  sql,
1480  0, NULL,
1481  NULL, NULL,
1482  TRUE, 0
1483  );
1484  pfree(sql);
1485 
1486  /* process resultset */
1487  SPI_cursor_fetch(portal, TRUE, 1);
1488  while (SPI_processed == 1 && SPI_tuptable != NULL) {
1489  tupdesc = SPI_tuptable->tupdesc;
1490  tuptable = SPI_tuptable;
1491  tuple = tuptable->vals[0];
1492 
1493  datum = SPI_getbinval(tuple, tupdesc, 1, &isNull);
1494  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1495 
1496  if (SPI_tuptable) SPI_freetuptable(tuptable);
1497  SPI_cursor_close(portal);
1498  SPI_finish();
1499 
1500  if (NULL != covhist) pfree(covhist);
1501  if (bin_width_count) pfree(bin_width);
1502 
1503  MemoryContextSwitchTo(oldcontext);
1504  elog(ERROR, "RASTER_histogramCoverage: Cannot get raster of coverage");
1505  SRF_RETURN_DONE(funcctx);
1506  }
1507  else if (isNull) {
1508  SPI_cursor_fetch(portal, TRUE, 1);
1509  continue;
1510  }
1511 
1512  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(datum);
1513 
1514  raster = rt_raster_deserialize(pgraster, FALSE);
1515  if (!raster) {
1516 
1517  if (SPI_tuptable) SPI_freetuptable(tuptable);
1518  SPI_cursor_close(portal);
1519  SPI_finish();
1520 
1521  if (NULL != covhist) pfree(covhist);
1522  if (bin_width_count) pfree(bin_width);
1523 
1524  MemoryContextSwitchTo(oldcontext);
1525  elog(ERROR, "RASTER_histogramCoverage: Cannot deserialize raster");
1526  SRF_RETURN_DONE(funcctx);
1527  }
1528 
1529  /* inspect number of bands*/
1530  num_bands = rt_raster_get_num_bands(raster);
1531  if (bandindex < 1 || bandindex > num_bands) {
1532  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
1533 
1534  rt_raster_destroy(raster);
1535 
1536  if (SPI_tuptable) SPI_freetuptable(tuptable);
1537  SPI_cursor_close(portal);
1538  SPI_finish();
1539 
1540  if (NULL != covhist) pfree(covhist);
1541  if (bin_width_count) pfree(bin_width);
1542 
1543  MemoryContextSwitchTo(oldcontext);
1544  SRF_RETURN_DONE(funcctx);
1545  }
1546 
1547  /* get band */
1548  band = rt_raster_get_band(raster, bandindex - 1);
1549  if (!band) {
1550  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
1551 
1552  rt_raster_destroy(raster);
1553 
1554  if (SPI_tuptable) SPI_freetuptable(tuptable);
1555  SPI_cursor_close(portal);
1556  SPI_finish();
1557 
1558  if (NULL != covhist) pfree(covhist);
1559  if (bin_width_count) pfree(bin_width);
1560 
1561  MemoryContextSwitchTo(oldcontext);
1562  SRF_RETURN_DONE(funcctx);
1563  }
1564 
1565  /* we need the raw values, hence the non-zero parameter */
1566  stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 1, NULL, NULL, NULL);
1567 
1568  rt_band_destroy(band);
1569  rt_raster_destroy(raster);
1570 
1571  if (NULL == stats) {
1572  elog(NOTICE, "Cannot compute summary statistics for band at index %d. Returning NULL", bandindex);
1573 
1574  if (SPI_tuptable) SPI_freetuptable(tuptable);
1575  SPI_cursor_close(portal);
1576  SPI_finish();
1577 
1578  if (NULL != covhist) pfree(covhist);
1579  if (bin_width_count) pfree(bin_width);
1580 
1581  MemoryContextSwitchTo(oldcontext);
1582  SRF_RETURN_DONE(funcctx);
1583  }
1584 
1585  /* get histogram */
1586  if (stats->count > 0) {
1587  hist = rt_band_get_histogram(stats, bin_count, bin_width, bin_width_count, right, min, max, &count);
1588  pfree(stats);
1589  if (NULL == hist || !count) {
1590  elog(NOTICE, "Cannot compute histogram for band at index %d", bandindex);
1591 
1592  if (SPI_tuptable) SPI_freetuptable(tuptable);
1593  SPI_cursor_close(portal);
1594  SPI_finish();
1595 
1596  if (NULL != covhist) pfree(covhist);
1597  if (bin_width_count) pfree(bin_width);
1598 
1599  MemoryContextSwitchTo(oldcontext);
1600  SRF_RETURN_DONE(funcctx);
1601  }
1602 
1603  POSTGIS_RT_DEBUGF(3, "%d bins returned", count);
1604 
1605  /* coverage histogram */
1606  if (NULL == covhist) {
1607  covhist = (rt_histogram) SPI_palloc(sizeof(struct rt_histogram_t) * count);
1608  if (NULL == covhist) {
1609 
1610  pfree(hist);
1611  if (SPI_tuptable) SPI_freetuptable(tuptable);
1612  SPI_cursor_close(portal);
1613  SPI_finish();
1614 
1615  if (bin_width_count) pfree(bin_width);
1616 
1617  MemoryContextSwitchTo(oldcontext);
1618  elog(ERROR, "RASTER_histogramCoverage: Cannot allocate memory for histogram of coverage");
1619  SRF_RETURN_DONE(funcctx);
1620  }
1621 
1622  for (i = 0; i < count; i++) {
1623  sum += hist[i].count;
1624  covhist[i].count = hist[i].count;
1625  covhist[i].percent = 0;
1626  covhist[i].min = hist[i].min;
1627  covhist[i].max = hist[i].max;
1628  }
1629  }
1630  else {
1631  for (i = 0; i < count; i++) {
1632  sum += hist[i].count;
1633  covhist[i].count += hist[i].count;
1634  }
1635  }
1636 
1637  pfree(hist);
1638 
1639  /* assuming bin_count wasn't set, force consistency */
1640  if (bin_count <= 0) bin_count = count;
1641  }
1642 
1643  /* next record */
1644  SPI_cursor_fetch(portal, TRUE, 1);
1645  }
1646 
1647  if (SPI_tuptable) SPI_freetuptable(tuptable);
1648  SPI_cursor_close(portal);
1649  SPI_finish();
1650 
1651  if (bin_width_count) pfree(bin_width);
1652 
1653  /* finish percent of histogram */
1654  if (sum > 0) {
1655  for (i = 0; i < count; i++)
1656  covhist[i].percent = covhist[i].count / (double) sum;
1657  }
1658 
1659  /* Store needed information */
1660  funcctx->user_fctx = covhist;
1661 
1662  /* total number of tuples to be returned */
1663  funcctx->max_calls = count;
1664 
1665  /* Build a tuple descriptor for our result type */
1666  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
1667  ereport(ERROR, (
1668  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1669  errmsg(
1670  "function returning record called in context "
1671  "that cannot accept type record"
1672  )
1673  ));
1674  }
1675 
1676  BlessTupleDesc(tupdesc);
1677  funcctx->tuple_desc = tupdesc;
1678 
1679  MemoryContextSwitchTo(oldcontext);
1680  }
1681 
1682  /* stuff done on every call of the function */
1683  funcctx = SRF_PERCALL_SETUP();
1684 
1685  call_cntr = funcctx->call_cntr;
1686  max_calls = funcctx->max_calls;
1687  tupdesc = funcctx->tuple_desc;
1688  covhist2 = funcctx->user_fctx;
1689 
1690  /* do when there is more left to send */
1691  if (call_cntr < max_calls) {
1692  Datum values[VALUES_LENGTH];
1693  bool nulls[VALUES_LENGTH];
1694  HeapTuple tuple;
1695  Datum result;
1696 
1697  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
1698 
1699  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
1700 
1701  values[0] = Float8GetDatum(covhist2[call_cntr].min);
1702  values[1] = Float8GetDatum(covhist2[call_cntr].max);
1703  values[2] = Int64GetDatum(covhist2[call_cntr].count);
1704  values[3] = Float8GetDatum(covhist2[call_cntr].percent);
1705 
1706  /* build a tuple */
1707  tuple = heap_form_tuple(tupdesc, values, nulls);
1708 
1709  /* make the tuple into a datum */
1710  result = HeapTupleGetDatum(tuple);
1711 
1712  SRF_RETURN_NEXT(funcctx, result);
1713  }
1714  /* do when there is no more left */
1715  else {
1716  pfree(covhist2);
1717  SRF_RETURN_DONE(funcctx);
1718  }
1719 }
1720 
1721 #undef VALUES_LENGTH
1722 #define VALUES_LENGTH 2
1723 
1727 PG_FUNCTION_INFO_V1(RASTER_quantile);
1728 Datum RASTER_quantile(PG_FUNCTION_ARGS)
1729 {
1730  FuncCallContext *funcctx;
1731  TupleDesc tupdesc;
1732 
1733  int i;
1734  rt_quantile quant;
1735  rt_quantile quant2;
1736  int call_cntr;
1737  int max_calls;
1738 
1739  /* first call of function */
1740  if (SRF_IS_FIRSTCALL()) {
1741  MemoryContext oldcontext;
1742 
1743  rt_pgraster *pgraster = NULL;
1744  rt_raster raster = NULL;
1745  rt_band band = NULL;
1746  int32_t bandindex = 0;
1747  int num_bands = 0;
1748  bool exclude_nodata_value = TRUE;
1749  double sample = 0;
1750  double *quantiles = NULL;
1751  uint32_t quantiles_count = 0;
1752  double quantile = 0;
1753  rt_bandstats stats = NULL;
1754  uint32_t count;
1755 
1756  int j;
1757  int n;
1758 
1759  ArrayType *array;
1760  Oid etype;
1761  Datum *e;
1762  bool *nulls;
1763  int16 typlen;
1764  bool typbyval;
1765  char typalign;
1766 
1767  /* create a function context for cross-call persistence */
1768  funcctx = SRF_FIRSTCALL_INIT();
1769 
1770  /* switch to memory context appropriate for multiple function calls */
1771  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
1772 
1773  /* pgraster is null, return nothing */
1774  if (PG_ARGISNULL(0)) {
1775  MemoryContextSwitchTo(oldcontext);
1776  SRF_RETURN_DONE(funcctx);
1777  }
1778  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1779 
1780  raster = rt_raster_deserialize(pgraster, FALSE);
1781  if (!raster) {
1782  PG_FREE_IF_COPY(pgraster, 0);
1783  MemoryContextSwitchTo(oldcontext);
1784  elog(ERROR, "RASTER_quantile: Cannot deserialize raster");
1785  SRF_RETURN_DONE(funcctx);
1786  }
1787 
1788  /* band index is 1-based */
1789  bandindex = PG_GETARG_INT32(1);
1790  num_bands = rt_raster_get_num_bands(raster);
1791  if (bandindex < 1 || bandindex > num_bands) {
1792  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
1793  rt_raster_destroy(raster);
1794  PG_FREE_IF_COPY(pgraster, 0);
1795  MemoryContextSwitchTo(oldcontext);
1796  SRF_RETURN_DONE(funcctx);
1797  }
1798 
1799  /* exclude_nodata_value flag */
1800  if (!PG_ARGISNULL(2))
1801  exclude_nodata_value = PG_GETARG_BOOL(2);
1802 
1803  /* sample % */
1804  if (!PG_ARGISNULL(3)) {
1805  sample = PG_GETARG_FLOAT8(3);
1806  if (sample < 0 || sample > 1) {
1807  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
1808  rt_raster_destroy(raster);
1809  PG_FREE_IF_COPY(pgraster, 0);
1810  MemoryContextSwitchTo(oldcontext);
1811  SRF_RETURN_DONE(funcctx);
1812  }
1813  else if (FLT_EQ(sample, 0.0))
1814  sample = 1;
1815  }
1816  else
1817  sample = 1;
1818 
1819  /* quantiles */
1820  if (!PG_ARGISNULL(4)) {
1821  array = PG_GETARG_ARRAYTYPE_P(4);
1822  etype = ARR_ELEMTYPE(array);
1823  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1824 
1825  switch (etype) {
1826  case FLOAT4OID:
1827  case FLOAT8OID:
1828  break;
1829  default:
1830  rt_raster_destroy(raster);
1831  PG_FREE_IF_COPY(pgraster, 0);
1832  MemoryContextSwitchTo(oldcontext);
1833  elog(ERROR, "RASTER_quantile: Invalid data type for quantiles");
1834  SRF_RETURN_DONE(funcctx);
1835  break;
1836  }
1837 
1838  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1839  &nulls, &n);
1840 
1841  quantiles = palloc(sizeof(double) * n);
1842  for (i = 0, j = 0; i < n; i++) {
1843  if (nulls[i]) continue;
1844 
1845  switch (etype) {
1846  case FLOAT4OID:
1847  quantile = (double) DatumGetFloat4(e[i]);
1848  break;
1849  case FLOAT8OID:
1850  quantile = (double) DatumGetFloat8(e[i]);
1851  break;
1852  }
1853 
1854  if (quantile < 0 || quantile > 1) {
1855  elog(NOTICE, "Invalid value for quantile (must be between 0 and 1). Returning NULL");
1856  pfree(quantiles);
1857  rt_raster_destroy(raster);
1858  PG_FREE_IF_COPY(pgraster, 0);
1859  MemoryContextSwitchTo(oldcontext);
1860  SRF_RETURN_DONE(funcctx);
1861  }
1862 
1863  quantiles[j] = quantile;
1864  POSTGIS_RT_DEBUGF(5, "quantiles[%d] = %f", j, quantiles[j]);
1865  j++;
1866  }
1867  quantiles_count = j;
1868 
1869  if (j < 1) {
1870  pfree(quantiles);
1871  quantiles = NULL;
1872  }
1873  }
1874 
1875  /* get band */
1876  band = rt_raster_get_band(raster, bandindex - 1);
1877  if (!band) {
1878  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
1879  rt_raster_destroy(raster);
1880  PG_FREE_IF_COPY(pgraster, 0);
1881  MemoryContextSwitchTo(oldcontext);
1882  SRF_RETURN_DONE(funcctx);
1883  }
1884 
1885  /* get stats */
1886  stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 1, NULL, NULL, NULL);
1887  rt_band_destroy(band);
1888  rt_raster_destroy(raster);
1889  PG_FREE_IF_COPY(pgraster, 0);
1890  if (NULL == stats || NULL == stats->values) {
1891  elog(NOTICE, "Cannot retrieve summary statistics for band at index %d", bandindex);
1892  MemoryContextSwitchTo(oldcontext);
1893  SRF_RETURN_DONE(funcctx);
1894  }
1895  else if (stats->count < 1) {
1896  elog(NOTICE, "Cannot compute quantiles for band at index %d as the band has no values", bandindex);
1897  MemoryContextSwitchTo(oldcontext);
1898  SRF_RETURN_DONE(funcctx);
1899  }
1900 
1901  /* get quantiles */
1902  quant = rt_band_get_quantiles(stats, quantiles, quantiles_count, &count);
1903  if (quantiles_count) pfree(quantiles);
1904  pfree(stats);
1905  if (NULL == quant || !count) {
1906  elog(NOTICE, "Cannot compute quantiles for band at index %d", bandindex);
1907  MemoryContextSwitchTo(oldcontext);
1908  SRF_RETURN_DONE(funcctx);
1909  }
1910 
1911  POSTGIS_RT_DEBUGF(3, "%d quantiles returned", count);
1912 
1913  /* Store needed information */
1914  funcctx->user_fctx = quant;
1915 
1916  /* total number of tuples to be returned */
1917  funcctx->max_calls = count;
1918 
1919  /* Build a tuple descriptor for our result type */
1920  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
1921  ereport(ERROR, (
1922  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1923  errmsg(
1924  "function returning record called in context "
1925  "that cannot accept type record"
1926  )
1927  ));
1928  }
1929 
1930  BlessTupleDesc(tupdesc);
1931  funcctx->tuple_desc = tupdesc;
1932 
1933  MemoryContextSwitchTo(oldcontext);
1934  }
1935 
1936  /* stuff done on every call of the function */
1937  funcctx = SRF_PERCALL_SETUP();
1938 
1939  call_cntr = funcctx->call_cntr;
1940  max_calls = funcctx->max_calls;
1941  tupdesc = funcctx->tuple_desc;
1942  quant2 = funcctx->user_fctx;
1943 
1944  /* do when there is more left to send */
1945  if (call_cntr < max_calls) {
1946  Datum values[VALUES_LENGTH];
1947  bool nulls[VALUES_LENGTH];
1948  HeapTuple tuple;
1949  Datum result;
1950 
1951  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
1952 
1953  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
1954 
1955  values[0] = Float8GetDatum(quant2[call_cntr].quantile);
1956  values[1] = Float8GetDatum(quant2[call_cntr].value);
1957 
1958  /* build a tuple */
1959  tuple = heap_form_tuple(tupdesc, values, nulls);
1960 
1961  /* make the tuple into a datum */
1962  result = HeapTupleGetDatum(tuple);
1963 
1964  SRF_RETURN_NEXT(funcctx, result);
1965  }
1966  /* do when there is no more left */
1967  else {
1968  pfree(quant2);
1969  SRF_RETURN_DONE(funcctx);
1970  }
1971 }
1972 
1976 PG_FUNCTION_INFO_V1(RASTER_quantileCoverage);
1977 Datum RASTER_quantileCoverage(PG_FUNCTION_ARGS)
1978 {
1979  FuncCallContext *funcctx;
1980  TupleDesc tupdesc;
1981 
1982  int i;
1983  rt_quantile covquant = NULL;
1984  rt_quantile covquant2;
1985  int call_cntr;
1986  int max_calls;
1987 
1988  POSTGIS_RT_DEBUG(3, "RASTER_quantileCoverage: Starting");
1989 
1990  /* first call of function */
1991  if (SRF_IS_FIRSTCALL()) {
1992  MemoryContext oldcontext;
1993 
1994  text *tablenametext = NULL;
1995  char *tablename = NULL;
1996  text *colnametext = NULL;
1997  char *colname = NULL;
1998  int32_t bandindex = 1;
1999  bool exclude_nodata_value = TRUE;
2000  double sample = 0;
2001  double *quantiles = NULL;
2002  uint32_t quantiles_count = 0;
2003  double quantile = 0;
2004  uint32_t count;
2005 
2006  int len = 0;
2007  char *sql = NULL;
2008  char *tmp = NULL;
2009  uint64_t cov_count = 0;
2010  int spi_result;
2011  Portal portal;
2012  SPITupleTable *tuptable = NULL;
2013  HeapTuple tuple;
2014  Datum datum;
2015  bool isNull = FALSE;
2016 
2017  rt_pgraster *pgraster = NULL;
2018  rt_raster raster = NULL;
2019  rt_band band = NULL;
2020  int num_bands = 0;
2021  struct quantile_llist *qlls = NULL;
2022  uint32_t qlls_count;
2023 
2024  int j;
2025  int n;
2026 
2027  ArrayType *array;
2028  Oid etype;
2029  Datum *e;
2030  bool *nulls;
2031  int16 typlen;
2032  bool typbyval;
2033  char typalign;
2034 
2035  POSTGIS_RT_DEBUG(3, "RASTER_quantileCoverage: first call of function");
2036 
2037  /* create a function context for cross-call persistence */
2038  funcctx = SRF_FIRSTCALL_INIT();
2039 
2040  /* switch to memory context appropriate for multiple function calls */
2041  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
2042 
2043  /* tablename is null, return null */
2044  if (PG_ARGISNULL(0)) {
2045  elog(NOTICE, "Table name must be provided");
2046  MemoryContextSwitchTo(oldcontext);
2047  SRF_RETURN_DONE(funcctx);
2048  }
2049  tablenametext = PG_GETARG_TEXT_P(0);
2050  tablename = text_to_cstring(tablenametext);
2051  if (!strlen(tablename)) {
2052  elog(NOTICE, "Table name must be provided");
2053  MemoryContextSwitchTo(oldcontext);
2054  SRF_RETURN_DONE(funcctx);
2055  }
2056  POSTGIS_RT_DEBUGF(3, "RASTER_quantileCoverage: tablename = %s", tablename);
2057 
2058  /* column name is null, return null */
2059  if (PG_ARGISNULL(1)) {
2060  elog(NOTICE, "Column name must be provided");
2061  MemoryContextSwitchTo(oldcontext);
2062  SRF_RETURN_DONE(funcctx);
2063  }
2064  colnametext = PG_GETARG_TEXT_P(1);
2065  colname = text_to_cstring(colnametext);
2066  if (!strlen(colname)) {
2067  elog(NOTICE, "Column name must be provided");
2068  MemoryContextSwitchTo(oldcontext);
2069  SRF_RETURN_DONE(funcctx);
2070  }
2071  POSTGIS_RT_DEBUGF(3, "RASTER_quantileCoverage: colname = %s", colname);
2072 
2073  /* band index is 1-based */
2074  if (!PG_ARGISNULL(2))
2075  bandindex = PG_GETARG_INT32(2);
2076 
2077  /* exclude_nodata_value flag */
2078  if (!PG_ARGISNULL(3))
2079  exclude_nodata_value = PG_GETARG_BOOL(3);
2080 
2081  /* sample % */
2082  if (!PG_ARGISNULL(4)) {
2083  sample = PG_GETARG_FLOAT8(4);
2084  if (sample < 0 || sample > 1) {
2085  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
2086  MemoryContextSwitchTo(oldcontext);
2087  SRF_RETURN_DONE(funcctx);
2088  }
2089  else if (FLT_EQ(sample, 0.0))
2090  sample = 1;
2091  }
2092  else
2093  sample = 1;
2094 
2095  /* quantiles */
2096  if (!PG_ARGISNULL(5)) {
2097  array = PG_GETARG_ARRAYTYPE_P(5);
2098  etype = ARR_ELEMTYPE(array);
2099  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
2100 
2101  switch (etype) {
2102  case FLOAT4OID:
2103  case FLOAT8OID:
2104  break;
2105  default:
2106  MemoryContextSwitchTo(oldcontext);
2107  elog(ERROR, "RASTER_quantileCoverage: Invalid data type for quantiles");
2108  SRF_RETURN_DONE(funcctx);
2109  break;
2110  }
2111 
2112  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
2113  &nulls, &n);
2114 
2115  quantiles = palloc(sizeof(double) * n);
2116  for (i = 0, j = 0; i < n; i++) {
2117  if (nulls[i]) continue;
2118 
2119  switch (etype) {
2120  case FLOAT4OID:
2121  quantile = (double) DatumGetFloat4(e[i]);
2122  break;
2123  case FLOAT8OID:
2124  quantile = (double) DatumGetFloat8(e[i]);
2125  break;
2126  }
2127 
2128  if (quantile < 0 || quantile > 1) {
2129  elog(NOTICE, "Invalid value for quantile (must be between 0 and 1). Returning NULL");
2130  pfree(quantiles);
2131  MemoryContextSwitchTo(oldcontext);
2132  SRF_RETURN_DONE(funcctx);
2133  }
2134 
2135  quantiles[j] = quantile;
2136  POSTGIS_RT_DEBUGF(5, "quantiles[%d] = %f", j, quantiles[j]);
2137  j++;
2138  }
2139  quantiles_count = j;
2140 
2141  if (j < 1) {
2142  pfree(quantiles);
2143  quantiles = NULL;
2144  }
2145  }
2146 
2147  /* coverage stats */
2148  /* connect to database */
2149  spi_result = SPI_connect();
2150  if (spi_result != SPI_OK_CONNECT) {
2151  MemoryContextSwitchTo(oldcontext);
2152  elog(ERROR, "RASTER_quantileCoverage: Cannot connect to database using SPI");
2153  SRF_RETURN_DONE(funcctx);
2154  }
2155 
2156  len = sizeof(char) * (strlen("SELECT count FROM _st_summarystats('','',,::boolean,)") + strlen(tablename) + strlen(colname) + (MAX_INT_CHARLEN * 2) + MAX_DBL_CHARLEN + 1);
2157  sql = (char *) palloc(len);
2158  if (NULL == sql) {
2159 
2160  if (SPI_tuptable) SPI_freetuptable(tuptable);
2161  SPI_finish();
2162 
2163  MemoryContextSwitchTo(oldcontext);
2164  elog(ERROR, "RASTER_quantileCoverage: Cannot allocate memory for sql");
2165  SRF_RETURN_DONE(funcctx);
2166  }
2167 
2168  /* get stats */
2169  snprintf(sql, len, "SELECT count FROM _st_summarystats('%s','%s',%d,%d::boolean,%f)", tablename, colname, bandindex, (exclude_nodata_value ? 1 : 0), sample);
2170  POSTGIS_RT_DEBUGF(3, "stats sql: %s", sql);
2171  spi_result = SPI_execute(sql, TRUE, 0);
2172  pfree(sql);
2173  if (spi_result != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
2174 
2175  if (SPI_tuptable) SPI_freetuptable(tuptable);
2176  SPI_finish();
2177 
2178  MemoryContextSwitchTo(oldcontext);
2179  elog(ERROR, "RASTER_quantileCoverage: Cannot get summary stats of coverage");
2180  SRF_RETURN_DONE(funcctx);
2181  }
2182 
2183  tupdesc = SPI_tuptable->tupdesc;
2184  tuptable = SPI_tuptable;
2185  tuple = tuptable->vals[0];
2186 
2187  tmp = SPI_getvalue(tuple, tupdesc, 1);
2188  if (NULL == tmp || !strlen(tmp)) {
2189 
2190  if (SPI_tuptable) SPI_freetuptable(tuptable);
2191  SPI_finish();
2192 
2193  MemoryContextSwitchTo(oldcontext);
2194  elog(ERROR, "RASTER_quantileCoverage: Cannot get summary stats of coverage");
2195  SRF_RETURN_DONE(funcctx);
2196  }
2197  cov_count = strtol(tmp, NULL, 10);
2198  POSTGIS_RT_DEBUGF(3, "covcount = %d", (int) cov_count);
2199  pfree(tmp);
2200 
2201  /* iterate through rasters of coverage */
2202  /* create sql */
2203  len = sizeof(char) * (strlen("SELECT \"\" FROM \"\" WHERE \"\" IS NOT NULL") + (strlen(colname) * 2) + strlen(tablename) + 1);
2204  sql = (char *) palloc(len);
2205  if (NULL == sql) {
2206 
2207  if (SPI_tuptable) SPI_freetuptable(tuptable);
2208  SPI_finish();
2209 
2210  MemoryContextSwitchTo(oldcontext);
2211  elog(ERROR, "RASTER_quantileCoverage: Cannot allocate memory for sql");
2212  SRF_RETURN_DONE(funcctx);
2213  }
2214 
2215  /* get cursor */
2216  snprintf(sql, len, "SELECT \"%s\" FROM \"%s\" WHERE \"%s\" IS NOT NULL", colname, tablename, colname);
2217  POSTGIS_RT_DEBUGF(3, "coverage sql: %s", sql);
2218  portal = SPI_cursor_open_with_args(
2219  "coverage",
2220  sql,
2221  0, NULL,
2222  NULL, NULL,
2223  TRUE, 0
2224  );
2225  pfree(sql);
2226 
2227  /* process resultset */
2228  SPI_cursor_fetch(portal, TRUE, 1);
2229  while (SPI_processed == 1 && SPI_tuptable != NULL) {
2230  if (NULL != covquant) pfree(covquant);
2231 
2232  tupdesc = SPI_tuptable->tupdesc;
2233  tuptable = SPI_tuptable;
2234  tuple = tuptable->vals[0];
2235 
2236  datum = SPI_getbinval(tuple, tupdesc, 1, &isNull);
2237  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
2238 
2239  if (SPI_tuptable) SPI_freetuptable(tuptable);
2240  SPI_cursor_close(portal);
2241  SPI_finish();
2242 
2243  MemoryContextSwitchTo(oldcontext);
2244  elog(ERROR, "RASTER_quantileCoverage: Cannot get raster of coverage");
2245  SRF_RETURN_DONE(funcctx);
2246  }
2247  else if (isNull) {
2248  SPI_cursor_fetch(portal, TRUE, 1);
2249  continue;
2250  }
2251 
2252  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(datum);
2253 
2254  raster = rt_raster_deserialize(pgraster, FALSE);
2255  if (!raster) {
2256 
2257  if (SPI_tuptable) SPI_freetuptable(tuptable);
2258  SPI_cursor_close(portal);
2259  SPI_finish();
2260 
2261  MemoryContextSwitchTo(oldcontext);
2262  elog(ERROR, "RASTER_quantileCoverage: Cannot deserialize raster");
2263  SRF_RETURN_DONE(funcctx);
2264  }
2265 
2266  /* inspect number of bands*/
2267  num_bands = rt_raster_get_num_bands(raster);
2268  if (bandindex < 1 || bandindex > num_bands) {
2269  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
2270 
2271  rt_raster_destroy(raster);
2272 
2273  if (SPI_tuptable) SPI_freetuptable(tuptable);
2274  SPI_cursor_close(portal);
2275  SPI_finish();
2276 
2277  MemoryContextSwitchTo(oldcontext);
2278  SRF_RETURN_DONE(funcctx);
2279  }
2280 
2281  /* get band */
2282  band = rt_raster_get_band(raster, bandindex - 1);
2283  if (!band) {
2284  elog(NOTICE, "Cannot find raster band of index %d. Returning NULL", bandindex);
2285 
2286  rt_raster_destroy(raster);
2287 
2288  if (SPI_tuptable) SPI_freetuptable(tuptable);
2289  SPI_cursor_close(portal);
2290  SPI_finish();
2291 
2292  MemoryContextSwitchTo(oldcontext);
2293  SRF_RETURN_DONE(funcctx);
2294  }
2295 
2296  covquant = rt_band_get_quantiles_stream(
2297  band,
2298  exclude_nodata_value, sample, cov_count,
2299  &qlls, &qlls_count,
2300  quantiles, quantiles_count,
2301  &count
2302  );
2303 
2304  rt_band_destroy(band);
2305  rt_raster_destroy(raster);
2306 
2307  if (NULL == covquant || !count) {
2308  elog(NOTICE, "Cannot compute quantiles for band at index %d", bandindex);
2309 
2310  if (SPI_tuptable) SPI_freetuptable(tuptable);
2311  SPI_cursor_close(portal);
2312  SPI_finish();
2313 
2314  MemoryContextSwitchTo(oldcontext);
2315  SRF_RETURN_DONE(funcctx);
2316  }
2317 
2318  /* next record */
2319  SPI_cursor_fetch(portal, TRUE, 1);
2320  }
2321 
2322  covquant2 = SPI_palloc(sizeof(struct rt_quantile_t) * count);
2323  for (i = 0; i < count; i++) {
2324  covquant2[i].quantile = covquant[i].quantile;
2325  covquant2[i].has_value = covquant[i].has_value;
2326  if (covquant2[i].has_value)
2327  covquant2[i].value = covquant[i].value;
2328  }
2329 
2330  if (NULL != covquant) pfree(covquant);
2331  quantile_llist_destroy(&qlls, qlls_count);
2332 
2333  if (SPI_tuptable) SPI_freetuptable(tuptable);
2334  SPI_cursor_close(portal);
2335  SPI_finish();
2336 
2337  if (quantiles_count) pfree(quantiles);
2338 
2339  POSTGIS_RT_DEBUGF(3, "%d quantiles returned", count);
2340 
2341  /* Store needed information */
2342  funcctx->user_fctx = covquant2;
2343 
2344  /* total number of tuples to be returned */
2345  funcctx->max_calls = count;
2346 
2347  /* Build a tuple descriptor for our result type */
2348  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
2349  ereport(ERROR, (
2350  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
2351  errmsg(
2352  "function returning record called in context "
2353  "that cannot accept type record"
2354  )
2355  ));
2356  }
2357 
2358  BlessTupleDesc(tupdesc);
2359  funcctx->tuple_desc = tupdesc;
2360 
2361  MemoryContextSwitchTo(oldcontext);
2362  }
2363 
2364  /* stuff done on every call of the function */
2365  funcctx = SRF_PERCALL_SETUP();
2366 
2367  call_cntr = funcctx->call_cntr;
2368  max_calls = funcctx->max_calls;
2369  tupdesc = funcctx->tuple_desc;
2370  covquant2 = funcctx->user_fctx;
2371 
2372  /* do when there is more left to send */
2373  if (call_cntr < max_calls) {
2374  Datum values[VALUES_LENGTH];
2375  bool nulls[VALUES_LENGTH];
2376  HeapTuple tuple;
2377  Datum result;
2378 
2379  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
2380 
2381  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
2382 
2383  values[0] = Float8GetDatum(covquant2[call_cntr].quantile);
2384  if (covquant2[call_cntr].has_value)
2385  values[1] = Float8GetDatum(covquant2[call_cntr].value);
2386  else
2387  nulls[1] = TRUE;
2388 
2389  /* build a tuple */
2390  tuple = heap_form_tuple(tupdesc, values, nulls);
2391 
2392  /* make the tuple into a datum */
2393  result = HeapTupleGetDatum(tuple);
2394 
2395  SRF_RETURN_NEXT(funcctx, result);
2396  }
2397  /* do when there is no more left */
2398  else {
2399  POSTGIS_RT_DEBUG(3, "done");
2400  pfree(covquant2);
2401  SRF_RETURN_DONE(funcctx);
2402  }
2403 }
2404 
2405 #undef VALUES_LENGTH
2406 #define VALUES_LENGTH 3
2407 
2408 /* get counts of values */
2409 PG_FUNCTION_INFO_V1(RASTER_valueCount);
2410 Datum RASTER_valueCount(PG_FUNCTION_ARGS) {
2411  FuncCallContext *funcctx;
2412  TupleDesc tupdesc;
2413 
2414  int i;
2415  rt_valuecount vcnts;
2416  rt_valuecount vcnts2;
2417  int call_cntr;
2418  int max_calls;
2419 
2420  /* first call of function */
2421  if (SRF_IS_FIRSTCALL()) {
2422  MemoryContext oldcontext;
2423 
2424  rt_pgraster *pgraster = NULL;
2425  rt_raster raster = NULL;
2426  rt_band band = NULL;
2427  int32_t bandindex = 0;
2428  int num_bands = 0;
2429  bool exclude_nodata_value = TRUE;
2430  double *search_values = NULL;
2431  uint32_t search_values_count = 0;
2432  double roundto = 0;
2433  uint32_t count;
2434 
2435  int j;
2436  int n;
2437 
2438  ArrayType *array;
2439  Oid etype;
2440  Datum *e;
2441  bool *nulls;
2442  int16 typlen;
2443  bool typbyval;
2444  char typalign;
2445 
2446  /* create a function context for cross-call persistence */
2447  funcctx = SRF_FIRSTCALL_INIT();
2448 
2449  /* switch to memory context appropriate for multiple function calls */
2450  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
2451 
2452  /* pgraster is null, return nothing */
2453  if (PG_ARGISNULL(0)) {
2454  MemoryContextSwitchTo(oldcontext);
2455  SRF_RETURN_DONE(funcctx);
2456  }
2457  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2458 
2459  raster = rt_raster_deserialize(pgraster, FALSE);
2460  if (!raster) {
2461  PG_FREE_IF_COPY(pgraster, 0);
2462  MemoryContextSwitchTo(oldcontext);
2463  elog(ERROR, "RASTER_valueCount: Cannot deserialize raster");
2464  SRF_RETURN_DONE(funcctx);
2465  }
2466 
2467  /* band index is 1-based */
2468  bandindex = PG_GETARG_INT32(1);
2469  num_bands = rt_raster_get_num_bands(raster);
2470  if (bandindex < 1 || bandindex > num_bands) {
2471  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
2472  rt_raster_destroy(raster);
2473  PG_FREE_IF_COPY(pgraster, 0);
2474  MemoryContextSwitchTo(oldcontext);
2475  SRF_RETURN_DONE(funcctx);
2476  }
2477 
2478  /* exclude_nodata_value flag */
2479  if (!PG_ARGISNULL(2))
2480  exclude_nodata_value = PG_GETARG_BOOL(2);
2481 
2482  /* search values */
2483  if (!PG_ARGISNULL(3)) {
2484  array = PG_GETARG_ARRAYTYPE_P(3);
2485  etype = ARR_ELEMTYPE(array);
2486  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
2487 
2488  switch (etype) {
2489  case FLOAT4OID:
2490  case FLOAT8OID:
2491  break;
2492  default:
2493  rt_raster_destroy(raster);
2494  PG_FREE_IF_COPY(pgraster, 0);
2495  MemoryContextSwitchTo(oldcontext);
2496  elog(ERROR, "RASTER_valueCount: Invalid data type for values");
2497  SRF_RETURN_DONE(funcctx);
2498  break;
2499  }
2500 
2501  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
2502  &nulls, &n);
2503 
2504  search_values = palloc(sizeof(double) * n);
2505  for (i = 0, j = 0; i < n; i++) {
2506  if (nulls[i]) continue;
2507 
2508  switch (etype) {
2509  case FLOAT4OID:
2510  search_values[j] = (double) DatumGetFloat4(e[i]);
2511  break;
2512  case FLOAT8OID:
2513  search_values[j] = (double) DatumGetFloat8(e[i]);
2514  break;
2515  }
2516 
2517  POSTGIS_RT_DEBUGF(5, "search_values[%d] = %f", j, search_values[j]);
2518  j++;
2519  }
2520  search_values_count = j;
2521 
2522  if (j < 1) {
2523  pfree(search_values);
2524  search_values = NULL;
2525  }
2526  }
2527 
2528  /* roundto */
2529  if (!PG_ARGISNULL(4)) {
2530  roundto = PG_GETARG_FLOAT8(4);
2531  if (roundto < 0.) roundto = 0;
2532  }
2533 
2534  /* get band */
2535  band = rt_raster_get_band(raster, bandindex - 1);
2536  if (!band) {
2537  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
2538  rt_raster_destroy(raster);
2539  PG_FREE_IF_COPY(pgraster, 0);
2540  MemoryContextSwitchTo(oldcontext);
2541  SRF_RETURN_DONE(funcctx);
2542  }
2543 
2544  /* get counts of values */
2545  vcnts = rt_band_get_value_count(band, (int) exclude_nodata_value, search_values, search_values_count, roundto, NULL, &count);
2546  rt_band_destroy(band);
2547  rt_raster_destroy(raster);
2548  PG_FREE_IF_COPY(pgraster, 0);
2549  if (NULL == vcnts || !count) {
2550  elog(NOTICE, "Cannot count the values for band at index %d", bandindex);
2551  MemoryContextSwitchTo(oldcontext);
2552  SRF_RETURN_DONE(funcctx);
2553  }
2554 
2555  POSTGIS_RT_DEBUGF(3, "%d value counts returned", count);
2556 
2557  /* Store needed information */
2558  funcctx->user_fctx = vcnts;
2559 
2560  /* total number of tuples to be returned */
2561  funcctx->max_calls = count;
2562 
2563  /* Build a tuple descriptor for our result type */
2564  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
2565  ereport(ERROR, (
2566  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
2567  errmsg(
2568  "function returning record called in context "
2569  "that cannot accept type record"
2570  )
2571  ));
2572  }
2573 
2574  BlessTupleDesc(tupdesc);
2575  funcctx->tuple_desc = tupdesc;
2576 
2577  MemoryContextSwitchTo(oldcontext);
2578  }
2579 
2580  /* stuff done on every call of the function */
2581  funcctx = SRF_PERCALL_SETUP();
2582 
2583  call_cntr = funcctx->call_cntr;
2584  max_calls = funcctx->max_calls;
2585  tupdesc = funcctx->tuple_desc;
2586  vcnts2 = funcctx->user_fctx;
2587 
2588  /* do when there is more left to send */
2589  if (call_cntr < max_calls) {
2590  Datum values[VALUES_LENGTH];
2591  bool nulls[VALUES_LENGTH];
2592  HeapTuple tuple;
2593  Datum result;
2594 
2595  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
2596 
2597  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
2598 
2599  values[0] = Float8GetDatum(vcnts2[call_cntr].value);
2600  values[1] = UInt32GetDatum(vcnts2[call_cntr].count);
2601  values[2] = Float8GetDatum(vcnts2[call_cntr].percent);
2602 
2603  /* build a tuple */
2604  tuple = heap_form_tuple(tupdesc, values, nulls);
2605 
2606  /* make the tuple into a datum */
2607  result = HeapTupleGetDatum(tuple);
2608 
2609  SRF_RETURN_NEXT(funcctx, result);
2610  }
2611  /* do when there is no more left */
2612  else {
2613  pfree(vcnts2);
2614  SRF_RETURN_DONE(funcctx);
2615  }
2616 }
2617 
2618 /* get counts of values for a coverage */
2619 PG_FUNCTION_INFO_V1(RASTER_valueCountCoverage);
2620 Datum RASTER_valueCountCoverage(PG_FUNCTION_ARGS) {
2621  FuncCallContext *funcctx;
2622  TupleDesc tupdesc;
2623 
2624  int i;
2625  uint64_t covcount = 0;
2626  uint64_t covtotal = 0;
2627  rt_valuecount covvcnts = NULL;
2628  rt_valuecount covvcnts2;
2629  int call_cntr;
2630  int max_calls;
2631 
2632  POSTGIS_RT_DEBUG(3, "RASTER_valueCountCoverage: Starting");
2633 
2634  /* first call of function */
2635  if (SRF_IS_FIRSTCALL()) {
2636  MemoryContext oldcontext;
2637 
2638  text *tablenametext = NULL;
2639  char *tablename = NULL;
2640  text *colnametext = NULL;
2641  char *colname = NULL;
2642  int32_t bandindex = 1;
2643  bool exclude_nodata_value = TRUE;
2644  double *search_values = NULL;
2645  uint32_t search_values_count = 0;
2646  double roundto = 0;
2647 
2648  int len = 0;
2649  char *sql = NULL;
2650  int spi_result;
2651  Portal portal;
2652  SPITupleTable *tuptable = NULL;
2653  HeapTuple tuple;
2654  Datum datum;
2655  bool isNull = FALSE;
2656  rt_pgraster *pgraster = NULL;
2657  rt_raster raster = NULL;
2658  rt_band band = NULL;
2659  int num_bands = 0;
2660  uint32_t count;
2661  uint32_t total;
2662  rt_valuecount vcnts = NULL;
2663  int exists = 0;
2664 
2665  int j;
2666  int n;
2667 
2668  ArrayType *array;
2669  Oid etype;
2670  Datum *e;
2671  bool *nulls;
2672  int16 typlen;
2673  bool typbyval;
2674  char typalign;
2675 
2676  /* create a function context for cross-call persistence */
2677  funcctx = SRF_FIRSTCALL_INIT();
2678 
2679  /* switch to memory context appropriate for multiple function calls */
2680  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
2681 
2682  /* tablename is null, return null */
2683  if (PG_ARGISNULL(0)) {
2684  elog(NOTICE, "Table name must be provided");
2685  MemoryContextSwitchTo(oldcontext);
2686  SRF_RETURN_DONE(funcctx);
2687  }
2688  tablenametext = PG_GETARG_TEXT_P(0);
2689  tablename = text_to_cstring(tablenametext);
2690  if (!strlen(tablename)) {
2691  elog(NOTICE, "Table name must be provided");
2692  MemoryContextSwitchTo(oldcontext);
2693  SRF_RETURN_DONE(funcctx);
2694  }
2695  POSTGIS_RT_DEBUGF(3, "tablename = %s", tablename);
2696 
2697  /* column name is null, return null */
2698  if (PG_ARGISNULL(1)) {
2699  elog(NOTICE, "Column name must be provided");
2700  MemoryContextSwitchTo(oldcontext);
2701  SRF_RETURN_DONE(funcctx);
2702  }
2703  colnametext = PG_GETARG_TEXT_P(1);
2704  colname = text_to_cstring(colnametext);
2705  if (!strlen(colname)) {
2706  elog(NOTICE, "Column name must be provided");
2707  MemoryContextSwitchTo(oldcontext);
2708  SRF_RETURN_DONE(funcctx);
2709  }
2710  POSTGIS_RT_DEBUGF(3, "colname = %s", colname);
2711 
2712  /* band index is 1-based */
2713  if (!PG_ARGISNULL(2))
2714  bandindex = PG_GETARG_INT32(2);
2715 
2716  /* exclude_nodata_value flag */
2717  if (!PG_ARGISNULL(3))
2718  exclude_nodata_value = PG_GETARG_BOOL(3);
2719 
2720  /* search values */
2721  if (!PG_ARGISNULL(4)) {
2722  array = PG_GETARG_ARRAYTYPE_P(4);
2723  etype = ARR_ELEMTYPE(array);
2724  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
2725 
2726  switch (etype) {
2727  case FLOAT4OID:
2728  case FLOAT8OID:
2729  break;
2730  default:
2731  MemoryContextSwitchTo(oldcontext);
2732  elog(ERROR, "RASTER_valueCountCoverage: Invalid data type for values");
2733  SRF_RETURN_DONE(funcctx);
2734  break;
2735  }
2736 
2737  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
2738  &nulls, &n);
2739 
2740  search_values = palloc(sizeof(double) * n);
2741  for (i = 0, j = 0; i < n; i++) {
2742  if (nulls[i]) continue;
2743 
2744  switch (etype) {
2745  case FLOAT4OID:
2746  search_values[j] = (double) DatumGetFloat4(e[i]);
2747  break;
2748  case FLOAT8OID:
2749  search_values[j] = (double) DatumGetFloat8(e[i]);
2750  break;
2751  }
2752 
2753  POSTGIS_RT_DEBUGF(5, "search_values[%d] = %f", j, search_values[j]);
2754  j++;
2755  }
2756  search_values_count = j;
2757 
2758  if (j < 1) {
2759  pfree(search_values);
2760  search_values = NULL;
2761  }
2762  }
2763 
2764  /* roundto */
2765  if (!PG_ARGISNULL(5)) {
2766  roundto = PG_GETARG_FLOAT8(5);
2767  if (roundto < 0.) roundto = 0;
2768  }
2769 
2770  /* iterate through rasters of coverage */
2771  /* connect to database */
2772  spi_result = SPI_connect();
2773  if (spi_result != SPI_OK_CONNECT) {
2774 
2775  if (search_values_count) pfree(search_values);
2776 
2777  MemoryContextSwitchTo(oldcontext);
2778  elog(ERROR, "RASTER_valueCountCoverage: Cannot connect to database using SPI");
2779  SRF_RETURN_DONE(funcctx);
2780  }
2781 
2782  /* create sql */
2783  len = sizeof(char) * (strlen("SELECT \"\" FROM \"\" WHERE \"\" IS NOT NULL") + (strlen(colname) * 2) + strlen(tablename) + 1);
2784  sql = (char *) palloc(len);
2785  if (NULL == sql) {
2786 
2787  if (SPI_tuptable) SPI_freetuptable(tuptable);
2788  SPI_finish();
2789 
2790  if (search_values_count) pfree(search_values);
2791 
2792  MemoryContextSwitchTo(oldcontext);
2793  elog(ERROR, "RASTER_valueCountCoverage: Cannot allocate memory for sql");
2794  SRF_RETURN_DONE(funcctx);
2795  }
2796 
2797  /* get cursor */
2798  snprintf(sql, len, "SELECT \"%s\" FROM \"%s\" WHERE \"%s\" IS NOT NULL", colname, tablename, colname);
2799  POSTGIS_RT_DEBUGF(3, "RASTER_valueCountCoverage: %s", sql);
2800  portal = SPI_cursor_open_with_args(
2801  "coverage",
2802  sql,
2803  0, NULL,
2804  NULL, NULL,
2805  TRUE, 0
2806  );
2807  pfree(sql);
2808 
2809  /* process resultset */
2810  SPI_cursor_fetch(portal, TRUE, 1);
2811  while (SPI_processed == 1 && SPI_tuptable != NULL) {
2812  tupdesc = SPI_tuptable->tupdesc;
2813  tuptable = SPI_tuptable;
2814  tuple = tuptable->vals[0];
2815 
2816  datum = SPI_getbinval(tuple, tupdesc, 1, &isNull);
2817  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
2818 
2819  if (SPI_tuptable) SPI_freetuptable(tuptable);
2820  SPI_cursor_close(portal);
2821  SPI_finish();
2822 
2823  if (NULL != covvcnts) pfree(covvcnts);
2824  if (search_values_count) pfree(search_values);
2825 
2826  MemoryContextSwitchTo(oldcontext);
2827  elog(ERROR, "RASTER_valueCountCoverage: Cannot get raster of coverage");
2828  SRF_RETURN_DONE(funcctx);
2829  }
2830  else if (isNull) {
2831  SPI_cursor_fetch(portal, TRUE, 1);
2832  continue;
2833  }
2834 
2835  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(datum);
2836 
2837  raster = rt_raster_deserialize(pgraster, FALSE);
2838  if (!raster) {
2839 
2840  if (SPI_tuptable) SPI_freetuptable(tuptable);
2841  SPI_cursor_close(portal);
2842  SPI_finish();
2843 
2844  if (NULL != covvcnts) pfree(covvcnts);
2845  if (search_values_count) pfree(search_values);
2846 
2847  MemoryContextSwitchTo(oldcontext);
2848  elog(ERROR, "RASTER_valueCountCoverage: Cannot deserialize raster");
2849  SRF_RETURN_DONE(funcctx);
2850  }
2851 
2852  /* inspect number of bands*/
2853  num_bands = rt_raster_get_num_bands(raster);
2854  if (bandindex < 1 || bandindex > num_bands) {
2855  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
2856 
2857  rt_raster_destroy(raster);
2858 
2859  if (SPI_tuptable) SPI_freetuptable(tuptable);
2860  SPI_cursor_close(portal);
2861  SPI_finish();
2862 
2863  if (NULL != covvcnts) pfree(covvcnts);
2864  if (search_values_count) pfree(search_values);
2865 
2866  MemoryContextSwitchTo(oldcontext);
2867  SRF_RETURN_DONE(funcctx);
2868  }
2869 
2870  /* get band */
2871  band = rt_raster_get_band(raster, bandindex - 1);
2872  if (!band) {
2873  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
2874 
2875  rt_raster_destroy(raster);
2876 
2877  if (SPI_tuptable) SPI_freetuptable(tuptable);
2878  SPI_cursor_close(portal);
2879  SPI_finish();
2880 
2881  if (NULL != covvcnts) pfree(covvcnts);
2882  if (search_values_count) pfree(search_values);
2883 
2884  MemoryContextSwitchTo(oldcontext);
2885  SRF_RETURN_DONE(funcctx);
2886  }
2887 
2888  /* get counts of values */
2889  vcnts = rt_band_get_value_count(band, (int) exclude_nodata_value, search_values, search_values_count, roundto, &total, &count);
2890  rt_band_destroy(band);
2891  rt_raster_destroy(raster);
2892  if (NULL == vcnts || !count) {
2893  elog(NOTICE, "Cannot count the values for band at index %d", bandindex);
2894 
2895  if (SPI_tuptable) SPI_freetuptable(tuptable);
2896  SPI_cursor_close(portal);
2897  SPI_finish();
2898 
2899  if (NULL != covvcnts) free(covvcnts);
2900  if (search_values_count) pfree(search_values);
2901 
2902  MemoryContextSwitchTo(oldcontext);
2903  SRF_RETURN_DONE(funcctx);
2904  }
2905 
2906  POSTGIS_RT_DEBUGF(3, "%d value counts returned", count);
2907 
2908  if (NULL == covvcnts) {
2909  covvcnts = (rt_valuecount) SPI_palloc(sizeof(struct rt_valuecount_t) * count);
2910  if (NULL == covvcnts) {
2911 
2912  if (SPI_tuptable) SPI_freetuptable(tuptable);
2913  SPI_cursor_close(portal);
2914  SPI_finish();
2915 
2916  if (search_values_count) pfree(search_values);
2917 
2918  MemoryContextSwitchTo(oldcontext);
2919  elog(ERROR, "RASTER_valueCountCoverage: Cannot allocate memory for value counts of coverage");
2920  SRF_RETURN_DONE(funcctx);
2921  }
2922 
2923  for (i = 0; i < count; i++) {
2924  covvcnts[i].value = vcnts[i].value;
2925  covvcnts[i].count = vcnts[i].count;
2926  covvcnts[i].percent = -1;
2927  }
2928 
2929  covcount = count;
2930  }
2931  else {
2932  for (i = 0; i < count; i++) {
2933  exists = 0;
2934 
2935  for (j = 0; j < covcount; j++) {
2936  if (FLT_EQ(vcnts[i].value, covvcnts[j].value)) {
2937  exists = 1;
2938  break;
2939  }
2940  }
2941 
2942  if (exists) {
2943  covvcnts[j].count += vcnts[i].count;
2944  }
2945  else {
2946  covcount++;
2947  covvcnts = SPI_repalloc(covvcnts, sizeof(struct rt_valuecount_t) * covcount);
2948  if (NULL == covvcnts) {
2949 
2950  if (SPI_tuptable) SPI_freetuptable(tuptable);
2951  SPI_cursor_close(portal);
2952  SPI_finish();
2953 
2954  if (search_values_count) pfree(search_values);
2955  if (NULL != covvcnts) free(covvcnts);
2956 
2957  MemoryContextSwitchTo(oldcontext);
2958  elog(ERROR, "RASTER_valueCountCoverage: Cannot change allocated memory for value counts of coverage");
2959  SRF_RETURN_DONE(funcctx);
2960  }
2961 
2962  covvcnts[covcount - 1].value = vcnts[i].value;
2963  covvcnts[covcount - 1].count = vcnts[i].count;
2964  covvcnts[covcount - 1].percent = -1;
2965  }
2966  }
2967  }
2968 
2969  covtotal += total;
2970 
2971  pfree(vcnts);
2972 
2973  /* next record */
2974  SPI_cursor_fetch(portal, TRUE, 1);
2975  }
2976 
2977  if (SPI_tuptable) SPI_freetuptable(tuptable);
2978  SPI_cursor_close(portal);
2979  SPI_finish();
2980 
2981  if (search_values_count) pfree(search_values);
2982 
2983  /* compute percentages */
2984  for (i = 0; i < covcount; i++) {
2985  covvcnts[i].percent = (double) covvcnts[i].count / covtotal;
2986  }
2987 
2988  /* Store needed information */
2989  funcctx->user_fctx = covvcnts;
2990 
2991  /* total number of tuples to be returned */
2992  funcctx->max_calls = covcount;
2993 
2994  /* Build a tuple descriptor for our result type */
2995  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
2996  ereport(ERROR, (
2997  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
2998  errmsg(
2999  "function returning record called in context "
3000  "that cannot accept type record"
3001  )
3002  ));
3003  }
3004 
3005  BlessTupleDesc(tupdesc);
3006  funcctx->tuple_desc = tupdesc;
3007 
3008  MemoryContextSwitchTo(oldcontext);
3009  }
3010 
3011  /* stuff done on every call of the function */
3012  funcctx = SRF_PERCALL_SETUP();
3013 
3014  call_cntr = funcctx->call_cntr;
3015  max_calls = funcctx->max_calls;
3016  tupdesc = funcctx->tuple_desc;
3017  covvcnts2 = funcctx->user_fctx;
3018 
3019  /* do when there is more left to send */
3020  if (call_cntr < max_calls) {
3021  Datum values[VALUES_LENGTH];
3022  bool nulls[VALUES_LENGTH];
3023  HeapTuple tuple;
3024  Datum result;
3025 
3026  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
3027 
3028  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
3029 
3030  values[0] = Float8GetDatum(covvcnts2[call_cntr].value);
3031  values[1] = UInt32GetDatum(covvcnts2[call_cntr].count);
3032  values[2] = Float8GetDatum(covvcnts2[call_cntr].percent);
3033 
3034  /* build a tuple */
3035  tuple = heap_form_tuple(tupdesc, values, nulls);
3036 
3037  /* make the tuple into a datum */
3038  result = HeapTupleGetDatum(tuple);
3039 
3040  SRF_RETURN_NEXT(funcctx, result);
3041  }
3042  /* do when there is no more left */
3043  else {
3044  pfree(covvcnts2);
3045  SRF_RETURN_DONE(funcctx);
3046  }
3047 }
3048 
RASTER_summaryStats
Datum RASTER_summaryStats(PG_FUNCTION_ARGS)
Definition: rtpg_statistics.c:73
rt_bandstats_t::mean
double mean
Definition: librtcore.h:2316
rtpg_summarystats_arg_init
static rtpg_summarystats_arg rtpg_summarystats_arg_init()
Definition: rtpg_statistics.c:522
VALUES_LENGTH
#define VALUES_LENGTH
Definition: rtpg_statistics.c:2406
rtpg_summarystats_arg_t::stats
rt_bandstats stats
Definition: rtpg_statistics.c:501
quantile_llist_destroy
int quantile_llist_destroy(struct quantile_llist **list, uint32_t list_count)
Definition: rt_statistics.c:878
rt_band_t
Definition: librtcore.h:2264
rt_raster_get_num_bands
int rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:372
rt_bandstats_t::count
uint32_t count
Definition: librtcore.h:2311
quantile_llist::quantile
double quantile
Definition: librtcore.h:2345
rtrowdump.raster
raster
Be careful!! Zeros function&#39;s input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
rt_valuecount_t::count
uint32_t count
Definition: librtcore.h:2376
rt_bandstats_t::sorted
int sorted
Definition: librtcore.h:2320
ovdump.band
band
Definition: ovdump.py:57
rt_quantile_t::quantile
double quantile
Definition: librtcore.h:2337
rt_bandstats_t::sum
double sum
Definition: librtcore.h:2315
FLT_EQ
#define FLT_EQ(x, y)
Definition: librtcore.h:2185
rt_histogram_t::max
double max
Definition: librtcore.h:2329
rtpg_summarystats_arg_t::cM
double cM
Definition: rtpg_statistics.c:505
RASTER_valueCount
Datum RASTER_valueCount(PG_FUNCTION_ARGS)
Definition: rtpg_statistics.c:2410
rt_band_destroy
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition: rt_band.c:242
rt_valuecount_t::value
double value
Definition: librtcore.h:2375
rtpg_summarystats_arg_t::sample
double sample
Definition: rtpg_statistics.c:510
rt_band_get_quantiles_stream
rt_quantile rt_band_get_quantiles_stream(rt_band band, int exclude_nodata_value, double sample, uint64_t cov_count, struct quantile_llist **qlls, uint32_t *qlls_count, double *quantiles, int quantiles_count, uint32_t *rtn_count)
Compute the default set of or requested quantiles for a coverage.
Definition: rt_statistics.c:1010
rtpg_summarystats_arg
struct rtpg_summarystats_arg_t * rtpg_summarystats_arg
Definition: rtpg_statistics.c:499
rtpg_summarystats_arg_t::cQ
double cQ
Definition: rtpg_statistics.c:506
rt_quantile_t::value
double value
Definition: librtcore.h:2338
POSTGIS_RT_DEBUGF
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition: rtpostgis.h:65
uint32_t
unsigned int uint32_t
Definition: uthash.h:78
MAX_INT_CHARLEN
#define MAX_INT_CHARLEN
Definition: rtpostgis.h:76
RASTER_summaryStats_transfn
Datum RASTER_summaryStats_transfn(PG_FUNCTION_ARGS)
Definition: rtpg_statistics.c:566
rt_bandstats_t::min
double min
Definition: librtcore.h:2313
rtpg_summarystats_arg_t
Definition: rtpg_statistics.c:500
rt_histogram_t::percent
double percent
Definition: librtcore.h:2326
quantile_llist
Definition: librtcore.h:2343
RASTER_histogramCoverage
Datum RASTER_histogramCoverage(PG_FUNCTION_ARGS)
Definition: rtpg_statistics.c:1195
rt_bandstats_t
Definition: librtcore.h:2309
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_raster_get_band
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:381
rt_quantile_t
Definition: librtcore.h:2336
genraster.count
int count
Definition: genraster.py:56
MAX_DBL_CHARLEN
#define MAX_DBL_CHARLEN
Definition: rtpostgis.h:75
PG_FUNCTION_INFO_V1
PG_FUNCTION_INFO_V1(RASTER_summaryStats)
Get summary stats of a band.
rt_bandstats_t::values
double * values
Definition: librtcore.h:2319
rtpg_summarystats_arg_t::cK
uint64_t cK
Definition: rtpg_statistics.c:504
rt_raster_destroy
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:82
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_histogram_t::min
double min
Definition: librtcore.h:2328
rt_histogram
struct rt_histogram_t * rt_histogram
Definition: librtcore.h:151
rt_bandstats_t::sample
double sample
Definition: librtcore.h:2310
rtpg_summarystats_arg_t::exclude_nodata_value
bool exclude_nodata_value
Definition: rtpg_statistics.c:509
rt_valuecount_t::percent
double percent
Definition: librtcore.h:2377
rt_valuecount
struct rt_valuecount_t * rt_valuecount
Definition: librtcore.h:153
rt_bandstats
struct rt_bandstats_t * rt_bandstats
Definition: librtcore.h:150
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
RASTER_summaryStats_finalfn
Datum RASTER_summaryStats_finalfn(PG_FUNCTION_ARGS)
Definition: rtpg_statistics.c:830
FALSE
#define FALSE
Definition: dbfopen.c:168
RASTER_valueCountCoverage
Datum RASTER_valueCountCoverage(PG_FUNCTION_ARGS)
Definition: rtpg_statistics.c:2620
rt_raster_serialized_t
Struct definitions.
Definition: librtcore.h:2201
rt_histogram_t
Definition: librtcore.h:2324
rtpg_summarystats_arg_destroy
static void rtpg_summarystats_arg_destroy(rtpg_summarystats_arg arg)
Definition: rtpg_statistics.c:514
RASTER_histogram
Datum RASTER_histogram(PG_FUNCTION_ARGS)
Definition: rtpg_statistics.c:921
rt_band_get_histogram
rt_histogram rt_band_get_histogram(rt_bandstats stats, int bin_count, double *bin_widths, int bin_widths_count, int right, double min, double max, uint32_t *rtn_count)
Count the distribution of data.
Definition: rt_statistics.c:412
rt_bandstats_t::max
double max
Definition: librtcore.h:2314
rt_quantile_t::has_value
uint32_t has_value
Definition: librtcore.h:2339
rtpg_summarystats_arg_t::band_index
int32_t band_index
Definition: rtpg_statistics.c:508
RASTER_quantile
Datum RASTER_quantile(PG_FUNCTION_ARGS)
Definition: rtpg_statistics.c:1728
free
void free(void *)
RASTER_quantileCoverage
Datum RASTER_quantileCoverage(PG_FUNCTION_ARGS)
Definition: rtpg_statistics.c:1977
rt_raster_t
Definition: librtcore.h:2235
genraster.value
int value
Definition: genraster.py:61
rt_bandstats_t::stddev
double stddev
Definition: librtcore.h:2317
POSTGIS_RT_DEBUG
#define POSTGIS_RT_DEBUG(level, msg)
Definition: rtpostgis.h:61
RASTER_summaryStatsCoverage
Datum RASTER_summaryStatsCoverage(PG_FUNCTION_ARGS)
Definition: rtpg_statistics.c:198
TRUE
#define TRUE
Definition: dbfopen.c:169
rt_raster_deserialize
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition: rt_serialize.c:717
rt_valuecount_t
Definition: librtcore.h:2374
rt_histogram_t::count
uint32_t count
Definition: librtcore.h:2325