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