PostGIS  2.5.7dev-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  rtpg_summarystats_arg_destroy(state);
906 
907  PG_RETURN_DATUM(result);
908 }
909 
910 #undef VALUES_LENGTH
911 #define VALUES_LENGTH 4
912 
916 PG_FUNCTION_INFO_V1(RASTER_histogram);
917 Datum RASTER_histogram(PG_FUNCTION_ARGS)
918 {
919  FuncCallContext *funcctx;
920  TupleDesc tupdesc;
921 
922  int i;
923  rt_histogram hist;
924  rt_histogram hist2;
925  int call_cntr;
926  int max_calls;
927 
928  /* first call of function */
929  if (SRF_IS_FIRSTCALL()) {
930  MemoryContext oldcontext;
931 
932  rt_pgraster *pgraster = NULL;
933  rt_raster raster = NULL;
934  rt_band band = NULL;
935  int32_t bandindex = 1;
936  int num_bands = 0;
937  bool exclude_nodata_value = TRUE;
938  double sample = 0;
939  uint32_t bin_count = 0;
940  double *bin_width = NULL;
941  uint32_t bin_width_count = 0;
942  double width = 0;
943  bool right = FALSE;
944  double min = 0;
945  double max = 0;
946  rt_bandstats stats = NULL;
947  uint32_t count;
948 
949  int j;
950  int n;
951 
952  ArrayType *array;
953  Oid etype;
954  Datum *e;
955  bool *nulls;
956  int16 typlen;
957  bool typbyval;
958  char typalign;
959 
960  POSTGIS_RT_DEBUG(3, "RASTER_histogram: Starting");
961 
962  /* create a function context for cross-call persistence */
963  funcctx = SRF_FIRSTCALL_INIT();
964 
965  /* switch to memory context appropriate for multiple function calls */
966  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
967 
968  /* pgraster is null, return nothing */
969  if (PG_ARGISNULL(0)) {
970  MemoryContextSwitchTo(oldcontext);
971  SRF_RETURN_DONE(funcctx);
972  }
973  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
974 
975  raster = rt_raster_deserialize(pgraster, FALSE);
976  if (!raster) {
977  PG_FREE_IF_COPY(pgraster, 0);
978  MemoryContextSwitchTo(oldcontext);
979  elog(ERROR, "RASTER_histogram: Cannot deserialize raster");
980  SRF_RETURN_DONE(funcctx);
981  }
982 
983  /* band index is 1-based */
984  if (!PG_ARGISNULL(1))
985  bandindex = PG_GETARG_INT32(1);
986  num_bands = rt_raster_get_num_bands(raster);
987  if (bandindex < 1 || bandindex > num_bands) {
988  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
989  rt_raster_destroy(raster);
990  PG_FREE_IF_COPY(pgraster, 0);
991  MemoryContextSwitchTo(oldcontext);
992  SRF_RETURN_DONE(funcctx);
993  }
994 
995  /* exclude_nodata_value flag */
996  if (!PG_ARGISNULL(2))
997  exclude_nodata_value = PG_GETARG_BOOL(2);
998 
999  /* sample % */
1000  if (!PG_ARGISNULL(3)) {
1001  sample = PG_GETARG_FLOAT8(3);
1002  if (sample < 0 || sample > 1) {
1003  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
1004  rt_raster_destroy(raster);
1005  PG_FREE_IF_COPY(pgraster, 0);
1006  MemoryContextSwitchTo(oldcontext);
1007  SRF_RETURN_DONE(funcctx);
1008  }
1009  else if (FLT_EQ(sample, 0.0))
1010  sample = 1;
1011  }
1012  else
1013  sample = 1;
1014 
1015  /* bin_count */
1016  if (!PG_ARGISNULL(4)) {
1017  bin_count = PG_GETARG_INT32(4);
1018  if (bin_count < 1) bin_count = 0;
1019  }
1020 
1021  /* bin_width */
1022  if (!PG_ARGISNULL(5)) {
1023  array = PG_GETARG_ARRAYTYPE_P(5);
1024  etype = ARR_ELEMTYPE(array);
1025  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1026 
1027  switch (etype) {
1028  case FLOAT4OID:
1029  case FLOAT8OID:
1030  break;
1031  default:
1032  rt_raster_destroy(raster);
1033  PG_FREE_IF_COPY(pgraster, 0);
1034  MemoryContextSwitchTo(oldcontext);
1035  elog(ERROR, "RASTER_histogram: Invalid data type for width");
1036  SRF_RETURN_DONE(funcctx);
1037  break;
1038  }
1039 
1040  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1041  &nulls, &n);
1042 
1043  bin_width = palloc(sizeof(double) * n);
1044  for (i = 0, j = 0; i < n; i++) {
1045  if (nulls[i]) continue;
1046 
1047  switch (etype) {
1048  case FLOAT4OID:
1049  width = (double) DatumGetFloat4(e[i]);
1050  break;
1051  case FLOAT8OID:
1052  width = (double) DatumGetFloat8(e[i]);
1053  break;
1054  }
1055 
1056  if (width < 0 || FLT_EQ(width, 0.0)) {
1057  elog(NOTICE, "Invalid value for width (must be greater than 0). Returning NULL");
1058  pfree(bin_width);
1059  rt_raster_destroy(raster);
1060  PG_FREE_IF_COPY(pgraster, 0);
1061  MemoryContextSwitchTo(oldcontext);
1062  SRF_RETURN_DONE(funcctx);
1063  }
1064 
1065  bin_width[j] = width;
1066  POSTGIS_RT_DEBUGF(5, "bin_width[%d] = %f", j, bin_width[j]);
1067  j++;
1068  }
1069  bin_width_count = j;
1070 
1071  if (j < 1) {
1072  pfree(bin_width);
1073  bin_width = NULL;
1074  }
1075  }
1076 
1077  /* right */
1078  if (!PG_ARGISNULL(6))
1079  right = PG_GETARG_BOOL(6);
1080 
1081  /* min */
1082  if (!PG_ARGISNULL(7)) min = PG_GETARG_FLOAT8(7);
1083 
1084  /* max */
1085  if (!PG_ARGISNULL(8)) max = PG_GETARG_FLOAT8(8);
1086 
1087  /* get band */
1088  band = rt_raster_get_band(raster, bandindex - 1);
1089  if (!band) {
1090  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
1091  rt_raster_destroy(raster);
1092  PG_FREE_IF_COPY(pgraster, 0);
1093  MemoryContextSwitchTo(oldcontext);
1094  SRF_RETURN_DONE(funcctx);
1095  }
1096 
1097  /* get stats */
1098  stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 1, NULL, NULL, NULL);
1099  rt_band_destroy(band);
1100  rt_raster_destroy(raster);
1101  PG_FREE_IF_COPY(pgraster, 0);
1102  if (NULL == stats || NULL == stats->values) {
1103  elog(NOTICE, "Cannot compute summary statistics for band at index %d", bandindex);
1104  MemoryContextSwitchTo(oldcontext);
1105  SRF_RETURN_DONE(funcctx);
1106  }
1107  else if (stats->count < 1) {
1108  elog(NOTICE, "Cannot compute histogram for band at index %d as the band has no values", bandindex);
1109  MemoryContextSwitchTo(oldcontext);
1110  SRF_RETURN_DONE(funcctx);
1111  }
1112 
1113  /* get histogram */
1114  hist = rt_band_get_histogram(stats, (uint32_t)bin_count, bin_width, bin_width_count, right, min, max, &count);
1115  if (bin_width_count) pfree(bin_width);
1116  pfree(stats);
1117  if (NULL == hist || !count) {
1118  elog(NOTICE, "Cannot compute histogram for band at index %d", bandindex);
1119  MemoryContextSwitchTo(oldcontext);
1120  SRF_RETURN_DONE(funcctx);
1121  }
1122 
1123  POSTGIS_RT_DEBUGF(3, "%d bins returned", count);
1124 
1125  /* Store needed information */
1126  funcctx->user_fctx = hist;
1127 
1128  /* total number of tuples to be returned */
1129  funcctx->max_calls = count;
1130 
1131  /* Build a tuple descriptor for our result type */
1132  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
1133  ereport(ERROR, (
1134  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1135  errmsg(
1136  "function returning record called in context "
1137  "that cannot accept type record"
1138  )
1139  ));
1140  }
1141 
1142  BlessTupleDesc(tupdesc);
1143  funcctx->tuple_desc = tupdesc;
1144 
1145  MemoryContextSwitchTo(oldcontext);
1146  }
1147 
1148  /* stuff done on every call of the function */
1149  funcctx = SRF_PERCALL_SETUP();
1150 
1151  call_cntr = funcctx->call_cntr;
1152  max_calls = funcctx->max_calls;
1153  tupdesc = funcctx->tuple_desc;
1154  hist2 = funcctx->user_fctx;
1155 
1156  /* do when there is more left to send */
1157  if (call_cntr < max_calls) {
1158  Datum values[VALUES_LENGTH];
1159  bool nulls[VALUES_LENGTH];
1160  HeapTuple tuple;
1161  Datum result;
1162 
1163  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
1164 
1165  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
1166 
1167  values[0] = Float8GetDatum(hist2[call_cntr].min);
1168  values[1] = Float8GetDatum(hist2[call_cntr].max);
1169  values[2] = Int64GetDatum(hist2[call_cntr].count);
1170  values[3] = Float8GetDatum(hist2[call_cntr].percent);
1171 
1172  /* build a tuple */
1173  tuple = heap_form_tuple(tupdesc, values, nulls);
1174 
1175  /* make the tuple into a datum */
1176  result = HeapTupleGetDatum(tuple);
1177 
1178  SRF_RETURN_NEXT(funcctx, result);
1179  }
1180  /* do when there is no more left */
1181  else {
1182  pfree(hist2);
1183  SRF_RETURN_DONE(funcctx);
1184  }
1185 }
1186 
1190 PG_FUNCTION_INFO_V1(RASTER_histogramCoverage);
1191 Datum RASTER_histogramCoverage(PG_FUNCTION_ARGS)
1192 {
1193  FuncCallContext *funcctx;
1194  TupleDesc tupdesc;
1195 
1196  uint32_t i;
1197  rt_histogram covhist = NULL;
1198  rt_histogram covhist2;
1199  int call_cntr;
1200  int max_calls;
1201 
1202  POSTGIS_RT_DEBUG(3, "RASTER_histogramCoverage: Starting");
1203 
1204  /* first call of function */
1205  if (SRF_IS_FIRSTCALL()) {
1206  MemoryContext oldcontext;
1207 
1208  text *tablenametext = NULL;
1209  char *tablename = NULL;
1210  text *colnametext = NULL;
1211  char *colname = NULL;
1212  int32_t bandindex = 1;
1213  bool exclude_nodata_value = TRUE;
1214  double sample = 0;
1215  uint32_t bin_count = 0;
1216  double *bin_width = NULL;
1217  uint32_t bin_width_count = 0;
1218  double width = 0;
1219  bool right = FALSE;
1220  uint32_t count;
1221 
1222  int len = 0;
1223  char *sql = NULL;
1224  char *tmp = NULL;
1225  double min = 0;
1226  double max = 0;
1227  int spi_result;
1228  Portal portal;
1229  HeapTuple tuple;
1230  Datum datum;
1231  bool isNull = FALSE;
1232 
1233  rt_pgraster *pgraster = NULL;
1234  rt_raster raster = NULL;
1235  rt_band band = NULL;
1236  int num_bands = 0;
1237  rt_bandstats stats = NULL;
1238  rt_histogram hist;
1239  uint64_t sum = 0;
1240 
1241  int j;
1242  int n;
1243 
1244  ArrayType *array;
1245  Oid etype;
1246  Datum *e;
1247  bool *nulls;
1248  int16 typlen;
1249  bool typbyval;
1250  char typalign;
1251 
1252  POSTGIS_RT_DEBUG(3, "RASTER_histogramCoverage: first call of function");
1253 
1254  /* create a function context for cross-call persistence */
1255  funcctx = SRF_FIRSTCALL_INIT();
1256 
1257  /* switch to memory context appropriate for multiple function calls */
1258  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
1259 
1260  /* tablename is null, return null */
1261  if (PG_ARGISNULL(0)) {
1262  elog(NOTICE, "Table name must be provided");
1263  MemoryContextSwitchTo(oldcontext);
1264  SRF_RETURN_DONE(funcctx);
1265  }
1266  tablenametext = PG_GETARG_TEXT_P(0);
1267  tablename = text_to_cstring(tablenametext);
1268  if (!strlen(tablename)) {
1269  elog(NOTICE, "Table name must be provided");
1270  MemoryContextSwitchTo(oldcontext);
1271  SRF_RETURN_DONE(funcctx);
1272  }
1273  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: tablename = %s", tablename);
1274 
1275  /* column name is null, return null */
1276  if (PG_ARGISNULL(1)) {
1277  elog(NOTICE, "Column name must be provided");
1278  MemoryContextSwitchTo(oldcontext);
1279  SRF_RETURN_DONE(funcctx);
1280  }
1281  colnametext = PG_GETARG_TEXT_P(1);
1282  colname = text_to_cstring(colnametext);
1283  if (!strlen(colname)) {
1284  elog(NOTICE, "Column name must be provided");
1285  MemoryContextSwitchTo(oldcontext);
1286  SRF_RETURN_DONE(funcctx);
1287  }
1288  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: colname = %s", colname);
1289 
1290  /* band index is 1-based */
1291  if (!PG_ARGISNULL(2))
1292  bandindex = PG_GETARG_INT32(2);
1293 
1294  /* exclude_nodata_value flag */
1295  if (!PG_ARGISNULL(3))
1296  exclude_nodata_value = PG_GETARG_BOOL(3);
1297 
1298  /* sample % */
1299  if (!PG_ARGISNULL(4)) {
1300  sample = PG_GETARG_FLOAT8(4);
1301  if (sample < 0 || sample > 1) {
1302  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
1303  MemoryContextSwitchTo(oldcontext);
1304  SRF_RETURN_DONE(funcctx);
1305  }
1306  else if (FLT_EQ(sample, 0.0))
1307  sample = 1;
1308  }
1309  else
1310  sample = 1;
1311 
1312  /* bin_count */
1313  if (!PG_ARGISNULL(5)) {
1314  bin_count = PG_GETARG_INT32(5);
1315  if (bin_count < 1) bin_count = 0;
1316  }
1317 
1318  /* bin_width */
1319  if (!PG_ARGISNULL(6)) {
1320  array = PG_GETARG_ARRAYTYPE_P(6);
1321  etype = ARR_ELEMTYPE(array);
1322  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1323 
1324  switch (etype) {
1325  case FLOAT4OID:
1326  case FLOAT8OID:
1327  break;
1328  default:
1329  MemoryContextSwitchTo(oldcontext);
1330  elog(ERROR, "RASTER_histogramCoverage: Invalid data type for width");
1331  SRF_RETURN_DONE(funcctx);
1332  break;
1333  }
1334 
1335  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1336  &nulls, &n);
1337 
1338  bin_width = palloc(sizeof(double) * n);
1339  for (i = 0, j = 0; i < (uint32_t) n; i++) {
1340  if (nulls[i]) continue;
1341 
1342  switch (etype) {
1343  case FLOAT4OID:
1344  width = (double) DatumGetFloat4(e[i]);
1345  break;
1346  case FLOAT8OID:
1347  width = (double) DatumGetFloat8(e[i]);
1348  break;
1349  }
1350 
1351  if (width < 0 || FLT_EQ(width, 0.0)) {
1352  elog(NOTICE, "Invalid value for width (must be greater than 0). Returning NULL");
1353  pfree(bin_width);
1354  MemoryContextSwitchTo(oldcontext);
1355  SRF_RETURN_DONE(funcctx);
1356  }
1357 
1358  bin_width[j] = width;
1359  POSTGIS_RT_DEBUGF(5, "bin_width[%d] = %f", j, bin_width[j]);
1360  j++;
1361  }
1362  bin_width_count = j;
1363 
1364  if (j < 1) {
1365  pfree(bin_width);
1366  bin_width = NULL;
1367  }
1368  }
1369 
1370  /* right */
1371  if (!PG_ARGISNULL(7))
1372  right = PG_GETARG_BOOL(7);
1373 
1374  /* connect to database */
1375  spi_result = SPI_connect();
1376  if (spi_result != SPI_OK_CONNECT) {
1377 
1378  if (bin_width_count) pfree(bin_width);
1379 
1380  MemoryContextSwitchTo(oldcontext);
1381  elog(ERROR, "RASTER_histogramCoverage: Cannot connect to database using SPI");
1382  SRF_RETURN_DONE(funcctx);
1383  }
1384 
1385  /* coverage stats */
1386  len = sizeof(char) * (strlen("SELECT min, max FROM _st_summarystats('','',,::boolean,)") + strlen(tablename) + strlen(colname) + (MAX_INT_CHARLEN * 2) + MAX_DBL_CHARLEN + 1);
1387  sql = (char *) palloc(len);
1388  if (NULL == sql) {
1389 
1390  if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
1391  SPI_finish();
1392 
1393  if (bin_width_count) pfree(bin_width);
1394 
1395  MemoryContextSwitchTo(oldcontext);
1396  elog(ERROR, "RASTER_histogramCoverage: Cannot allocate memory for sql");
1397  SRF_RETURN_DONE(funcctx);
1398  }
1399 
1400  /* get stats */
1401  snprintf(sql, len, "SELECT min, max FROM _st_summarystats('%s','%s',%d,%d::boolean,%f)", tablename, colname, bandindex, (exclude_nodata_value ? 1 : 0), sample);
1402  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: %s", sql);
1403  spi_result = SPI_execute(sql, TRUE, 0);
1404  pfree(sql);
1405  if (spi_result != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
1406 
1407  if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
1408  SPI_finish();
1409 
1410  if (bin_width_count) pfree(bin_width);
1411 
1412  MemoryContextSwitchTo(oldcontext);
1413  elog(ERROR, "RASTER_histogramCoverage: Cannot get summary stats of coverage");
1414  SRF_RETURN_DONE(funcctx);
1415  }
1416 
1417  tupdesc = SPI_tuptable->tupdesc;
1418  tuple = SPI_tuptable->vals[0];
1419 
1420  tmp = SPI_getvalue(tuple, tupdesc, 1);
1421  if (NULL == tmp || !strlen(tmp)) {
1422 
1423  SPI_freetuptable(SPI_tuptable);
1424  SPI_finish();
1425 
1426  if (bin_width_count) pfree(bin_width);
1427 
1428  MemoryContextSwitchTo(oldcontext);
1429  elog(ERROR, "RASTER_histogramCoverage: Cannot get summary stats of coverage");
1430  SRF_RETURN_DONE(funcctx);
1431  }
1432  min = strtod(tmp, NULL);
1433  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: min = %f", min);
1434  pfree(tmp);
1435 
1436  tmp = SPI_getvalue(tuple, tupdesc, 2);
1437  if (NULL == tmp || !strlen(tmp)) {
1438 
1439  if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
1440  SPI_finish();
1441 
1442  if (bin_width_count) pfree(bin_width);
1443 
1444  MemoryContextSwitchTo(oldcontext);
1445  elog(ERROR, "RASTER_histogramCoverage: Cannot get summary stats of coverage");
1446  SRF_RETURN_DONE(funcctx);
1447  }
1448  max = strtod(tmp, NULL);
1449  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: max = %f", max);
1450  pfree(tmp);
1451 
1452  /* iterate through rasters of coverage */
1453  /* create sql */
1454  len = sizeof(char) * (strlen("SELECT \"\" FROM \"\" WHERE \"\" IS NOT NULL") + (strlen(colname) * 2) + strlen(tablename) + 1);
1455  sql = (char *) palloc(len);
1456  if (NULL == sql) {
1457 
1458  if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
1459  SPI_finish();
1460 
1461  if (bin_width_count) pfree(bin_width);
1462 
1463  MemoryContextSwitchTo(oldcontext);
1464  elog(ERROR, "RASTER_histogramCoverage: Cannot allocate memory for sql");
1465  SRF_RETURN_DONE(funcctx);
1466  }
1467 
1468  /* get cursor */
1469  snprintf(sql, len, "SELECT \"%s\" FROM \"%s\" WHERE \"%s\" IS NOT NULL", colname, tablename, colname);
1470  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: %s", sql);
1471  portal = SPI_cursor_open_with_args(
1472  "coverage",
1473  sql,
1474  0, NULL,
1475  NULL, NULL,
1476  TRUE, 0
1477  );
1478  pfree(sql);
1479 
1480  /* process resultset */
1481  SPI_cursor_fetch(portal, TRUE, 1);
1482  while (SPI_processed == 1 && SPI_tuptable != NULL) {
1483  tupdesc = SPI_tuptable->tupdesc;
1484  tuple = SPI_tuptable->vals[0];
1485 
1486  datum = SPI_getbinval(tuple, tupdesc, 1, &isNull);
1487  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1488  SPI_freetuptable(SPI_tuptable);
1489  SPI_cursor_close(portal);
1490  SPI_finish();
1491 
1492  if (NULL != covhist) pfree(covhist);
1493  if (bin_width_count) pfree(bin_width);
1494 
1495  MemoryContextSwitchTo(oldcontext);
1496  elog(ERROR, "RASTER_histogramCoverage: Cannot get raster of coverage");
1497  SRF_RETURN_DONE(funcctx);
1498  }
1499  else if (isNull) {
1500  SPI_cursor_fetch(portal, TRUE, 1);
1501  continue;
1502  }
1503 
1504  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(datum);
1505 
1506  raster = rt_raster_deserialize(pgraster, FALSE);
1507  if (!raster) {
1508 
1509  if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
1510  SPI_cursor_close(portal);
1511  SPI_finish();
1512 
1513  if (NULL != covhist) pfree(covhist);
1514  if (bin_width_count) pfree(bin_width);
1515 
1516  MemoryContextSwitchTo(oldcontext);
1517  elog(ERROR, "RASTER_histogramCoverage: Cannot deserialize raster");
1518  SRF_RETURN_DONE(funcctx);
1519  }
1520 
1521  /* inspect number of bands*/
1522  num_bands = rt_raster_get_num_bands(raster);
1523  if (bandindex < 1 || bandindex > num_bands) {
1524  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
1525 
1526  rt_raster_destroy(raster);
1527 
1528  if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
1529  SPI_cursor_close(portal);
1530  SPI_finish();
1531 
1532  if (NULL != covhist) pfree(covhist);
1533  if (bin_width_count) pfree(bin_width);
1534 
1535  MemoryContextSwitchTo(oldcontext);
1536  SRF_RETURN_DONE(funcctx);
1537  }
1538 
1539  /* get band */
1540  band = rt_raster_get_band(raster, bandindex - 1);
1541  if (!band) {
1542  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
1543 
1544  rt_raster_destroy(raster);
1545 
1546  if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
1547  SPI_cursor_close(portal);
1548  SPI_finish();
1549 
1550  if (NULL != covhist) pfree(covhist);
1551  if (bin_width_count) pfree(bin_width);
1552 
1553  MemoryContextSwitchTo(oldcontext);
1554  SRF_RETURN_DONE(funcctx);
1555  }
1556 
1557  /* we need the raw values, hence the non-zero parameter */
1558  stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 1, NULL, NULL, NULL);
1559 
1560  rt_band_destroy(band);
1561  rt_raster_destroy(raster);
1562 
1563  if (NULL == stats) {
1564  elog(NOTICE, "Cannot compute summary statistics for band at index %d. Returning NULL", bandindex);
1565 
1566  if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
1567  SPI_cursor_close(portal);
1568  SPI_finish();
1569 
1570  if (NULL != covhist) pfree(covhist);
1571  if (bin_width_count) pfree(bin_width);
1572 
1573  MemoryContextSwitchTo(oldcontext);
1574  SRF_RETURN_DONE(funcctx);
1575  }
1576 
1577  /* get histogram */
1578  if (stats->count > 0) {
1579  hist = rt_band_get_histogram(stats, bin_count, bin_width, bin_width_count, right, min, max, &count);
1580  pfree(stats);
1581  if (NULL == hist || !count) {
1582  elog(NOTICE, "Cannot compute histogram for band at index %d", bandindex);
1583 
1584  if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
1585  SPI_cursor_close(portal);
1586  SPI_finish();
1587 
1588  if (NULL != covhist) pfree(covhist);
1589  if (bin_width_count) pfree(bin_width);
1590 
1591  MemoryContextSwitchTo(oldcontext);
1592  SRF_RETURN_DONE(funcctx);
1593  }
1594 
1595  POSTGIS_RT_DEBUGF(3, "%d bins returned", count);
1596 
1597  /* coverage histogram */
1598  if (NULL == covhist) {
1599  covhist = (rt_histogram) SPI_palloc(sizeof(struct rt_histogram_t) * count);
1600  if (NULL == covhist) {
1601 
1602  pfree(hist);
1603  if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
1604  SPI_cursor_close(portal);
1605  SPI_finish();
1606 
1607  if (bin_width_count) pfree(bin_width);
1608 
1609  MemoryContextSwitchTo(oldcontext);
1610  elog(ERROR, "RASTER_histogramCoverage: Cannot allocate memory for histogram of coverage");
1611  SRF_RETURN_DONE(funcctx);
1612  }
1613 
1614  for (i = 0; i < count; i++) {
1615  sum += hist[i].count;
1616  covhist[i].count = hist[i].count;
1617  covhist[i].percent = 0;
1618  covhist[i].min = hist[i].min;
1619  covhist[i].max = hist[i].max;
1620  }
1621  }
1622  else {
1623  for (i = 0; i < count; i++) {
1624  sum += hist[i].count;
1625  covhist[i].count += hist[i].count;
1626  }
1627  }
1628 
1629  pfree(hist);
1630 
1631  /* assuming bin_count wasn't set, force consistency */
1632  if (bin_count <= 0) bin_count = count;
1633  }
1634 
1635  /* next record */
1636  SPI_cursor_fetch(portal, TRUE, 1);
1637  }
1638 
1639  if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
1640  SPI_cursor_close(portal);
1641  SPI_finish();
1642 
1643  if (bin_width_count) pfree(bin_width);
1644 
1645  /* finish percent of histogram */
1646  if (sum > 0) {
1647  for (i = 0; i < count; i++)
1648  covhist[i].percent = covhist[i].count / (double) sum;
1649  }
1650 
1651  /* Store needed information */
1652  funcctx->user_fctx = covhist;
1653 
1654  /* total number of tuples to be returned */
1655  funcctx->max_calls = count;
1656 
1657  /* Build a tuple descriptor for our result type */
1658  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
1659  ereport(ERROR, (
1660  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1661  errmsg(
1662  "function returning record called in context "
1663  "that cannot accept type record"
1664  )
1665  ));
1666  }
1667 
1668  BlessTupleDesc(tupdesc);
1669  funcctx->tuple_desc = tupdesc;
1670 
1671  MemoryContextSwitchTo(oldcontext);
1672  }
1673 
1674  /* stuff done on every call of the function */
1675  funcctx = SRF_PERCALL_SETUP();
1676 
1677  call_cntr = funcctx->call_cntr;
1678  max_calls = funcctx->max_calls;
1679  tupdesc = funcctx->tuple_desc;
1680  covhist2 = funcctx->user_fctx;
1681 
1682  /* do when there is more left to send */
1683  if (call_cntr < max_calls) {
1684  Datum values[VALUES_LENGTH];
1685  bool nulls[VALUES_LENGTH];
1686  HeapTuple tuple;
1687  Datum result;
1688 
1689  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
1690 
1691  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
1692 
1693  values[0] = Float8GetDatum(covhist2[call_cntr].min);
1694  values[1] = Float8GetDatum(covhist2[call_cntr].max);
1695  values[2] = Int64GetDatum(covhist2[call_cntr].count);
1696  values[3] = Float8GetDatum(covhist2[call_cntr].percent);
1697 
1698  /* build a tuple */
1699  tuple = heap_form_tuple(tupdesc, values, nulls);
1700 
1701  /* make the tuple into a datum */
1702  result = HeapTupleGetDatum(tuple);
1703 
1704  SRF_RETURN_NEXT(funcctx, result);
1705  }
1706  /* do when there is no more left */
1707  else {
1708  pfree(covhist2);
1709  SRF_RETURN_DONE(funcctx);
1710  }
1711 }
1712 
1713 #undef VALUES_LENGTH
1714 #define VALUES_LENGTH 2
1715 
1719 PG_FUNCTION_INFO_V1(RASTER_quantile);
1720 Datum RASTER_quantile(PG_FUNCTION_ARGS)
1721 {
1722  FuncCallContext *funcctx;
1723  TupleDesc tupdesc;
1724 
1725  int i;
1726  rt_quantile quant;
1727  rt_quantile quant2;
1728  int call_cntr;
1729  int max_calls;
1730 
1731  /* first call of function */
1732  if (SRF_IS_FIRSTCALL()) {
1733  MemoryContext oldcontext;
1734 
1735  rt_pgraster *pgraster = NULL;
1736  rt_raster raster = NULL;
1737  rt_band band = NULL;
1738  int32_t bandindex = 0;
1739  int num_bands = 0;
1740  bool exclude_nodata_value = TRUE;
1741  double sample = 0;
1742  double *quantiles = NULL;
1743  uint32_t quantiles_count = 0;
1744  double quantile = 0;
1745  rt_bandstats stats = NULL;
1746  uint32_t count;
1747 
1748  int j;
1749  int n;
1750 
1751  ArrayType *array;
1752  Oid etype;
1753  Datum *e;
1754  bool *nulls;
1755  int16 typlen;
1756  bool typbyval;
1757  char typalign;
1758 
1759  /* create a function context for cross-call persistence */
1760  funcctx = SRF_FIRSTCALL_INIT();
1761 
1762  /* switch to memory context appropriate for multiple function calls */
1763  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
1764 
1765  /* pgraster is null, return nothing */
1766  if (PG_ARGISNULL(0)) {
1767  MemoryContextSwitchTo(oldcontext);
1768  SRF_RETURN_DONE(funcctx);
1769  }
1770  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1771 
1772  raster = rt_raster_deserialize(pgraster, FALSE);
1773  if (!raster) {
1774  PG_FREE_IF_COPY(pgraster, 0);
1775  MemoryContextSwitchTo(oldcontext);
1776  elog(ERROR, "RASTER_quantile: Cannot deserialize raster");
1777  SRF_RETURN_DONE(funcctx);
1778  }
1779 
1780  /* band index is 1-based */
1781  bandindex = PG_GETARG_INT32(1);
1782  num_bands = rt_raster_get_num_bands(raster);
1783  if (bandindex < 1 || bandindex > num_bands) {
1784  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
1785  rt_raster_destroy(raster);
1786  PG_FREE_IF_COPY(pgraster, 0);
1787  MemoryContextSwitchTo(oldcontext);
1788  SRF_RETURN_DONE(funcctx);
1789  }
1790 
1791  /* exclude_nodata_value flag */
1792  if (!PG_ARGISNULL(2))
1793  exclude_nodata_value = PG_GETARG_BOOL(2);
1794 
1795  /* sample % */
1796  if (!PG_ARGISNULL(3)) {
1797  sample = PG_GETARG_FLOAT8(3);
1798  if (sample < 0 || sample > 1) {
1799  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
1800  rt_raster_destroy(raster);
1801  PG_FREE_IF_COPY(pgraster, 0);
1802  MemoryContextSwitchTo(oldcontext);
1803  SRF_RETURN_DONE(funcctx);
1804  }
1805  else if (FLT_EQ(sample, 0.0))
1806  sample = 1;
1807  }
1808  else
1809  sample = 1;
1810 
1811  /* quantiles */
1812  if (!PG_ARGISNULL(4)) {
1813  array = PG_GETARG_ARRAYTYPE_P(4);
1814  etype = ARR_ELEMTYPE(array);
1815  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1816 
1817  switch (etype) {
1818  case FLOAT4OID:
1819  case FLOAT8OID:
1820  break;
1821  default:
1822  rt_raster_destroy(raster);
1823  PG_FREE_IF_COPY(pgraster, 0);
1824  MemoryContextSwitchTo(oldcontext);
1825  elog(ERROR, "RASTER_quantile: Invalid data type for quantiles");
1826  SRF_RETURN_DONE(funcctx);
1827  break;
1828  }
1829 
1830  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1831  &nulls, &n);
1832 
1833  quantiles = palloc(sizeof(double) * n);
1834  for (i = 0, j = 0; i < n; i++) {
1835  if (nulls[i]) continue;
1836 
1837  switch (etype) {
1838  case FLOAT4OID:
1839  quantile = (double) DatumGetFloat4(e[i]);
1840  break;
1841  case FLOAT8OID:
1842  quantile = (double) DatumGetFloat8(e[i]);
1843  break;
1844  }
1845 
1846  if (quantile < 0 || quantile > 1) {
1847  elog(NOTICE, "Invalid value for quantile (must be between 0 and 1). Returning NULL");
1848  pfree(quantiles);
1849  rt_raster_destroy(raster);
1850  PG_FREE_IF_COPY(pgraster, 0);
1851  MemoryContextSwitchTo(oldcontext);
1852  SRF_RETURN_DONE(funcctx);
1853  }
1854 
1855  quantiles[j] = quantile;
1856  POSTGIS_RT_DEBUGF(5, "quantiles[%d] = %f", j, quantiles[j]);
1857  j++;
1858  }
1859  quantiles_count = j;
1860 
1861  if (j < 1) {
1862  pfree(quantiles);
1863  quantiles = NULL;
1864  }
1865  }
1866 
1867  /* get band */
1868  band = rt_raster_get_band(raster, bandindex - 1);
1869  if (!band) {
1870  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
1871  rt_raster_destroy(raster);
1872  PG_FREE_IF_COPY(pgraster, 0);
1873  MemoryContextSwitchTo(oldcontext);
1874  SRF_RETURN_DONE(funcctx);
1875  }
1876 
1877  /* get stats */
1878  stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 1, NULL, NULL, NULL);
1879  rt_band_destroy(band);
1880  rt_raster_destroy(raster);
1881  PG_FREE_IF_COPY(pgraster, 0);
1882  if (NULL == stats || NULL == stats->values) {
1883  elog(NOTICE, "Cannot retrieve summary statistics for band at index %d", bandindex);
1884  MemoryContextSwitchTo(oldcontext);
1885  SRF_RETURN_DONE(funcctx);
1886  }
1887  else if (stats->count < 1) {
1888  elog(NOTICE, "Cannot compute quantiles for band at index %d as the band has no values", bandindex);
1889  MemoryContextSwitchTo(oldcontext);
1890  SRF_RETURN_DONE(funcctx);
1891  }
1892 
1893  /* get quantiles */
1894  quant = rt_band_get_quantiles(stats, quantiles, quantiles_count, &count);
1895  if (quantiles_count) pfree(quantiles);
1896  pfree(stats);
1897  if (NULL == quant || !count) {
1898  elog(NOTICE, "Cannot compute quantiles for band at index %d", bandindex);
1899  MemoryContextSwitchTo(oldcontext);
1900  SRF_RETURN_DONE(funcctx);
1901  }
1902 
1903  POSTGIS_RT_DEBUGF(3, "%d quantiles returned", count);
1904 
1905  /* Store needed information */
1906  funcctx->user_fctx = quant;
1907 
1908  /* total number of tuples to be returned */
1909  funcctx->max_calls = count;
1910 
1911  /* Build a tuple descriptor for our result type */
1912  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
1913  ereport(ERROR, (
1914  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1915  errmsg(
1916  "function returning record called in context "
1917  "that cannot accept type record"
1918  )
1919  ));
1920  }
1921 
1922  BlessTupleDesc(tupdesc);
1923  funcctx->tuple_desc = tupdesc;
1924 
1925  MemoryContextSwitchTo(oldcontext);
1926  }
1927 
1928  /* stuff done on every call of the function */
1929  funcctx = SRF_PERCALL_SETUP();
1930 
1931  call_cntr = funcctx->call_cntr;
1932  max_calls = funcctx->max_calls;
1933  tupdesc = funcctx->tuple_desc;
1934  quant2 = funcctx->user_fctx;
1935 
1936  /* do when there is more left to send */
1937  if (call_cntr < max_calls) {
1938  Datum values[VALUES_LENGTH];
1939  bool nulls[VALUES_LENGTH];
1940  HeapTuple tuple;
1941  Datum result;
1942 
1943  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
1944 
1945  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
1946 
1947  values[0] = Float8GetDatum(quant2[call_cntr].quantile);
1948  values[1] = Float8GetDatum(quant2[call_cntr].value);
1949 
1950  /* build a tuple */
1951  tuple = heap_form_tuple(tupdesc, values, nulls);
1952 
1953  /* make the tuple into a datum */
1954  result = HeapTupleGetDatum(tuple);
1955 
1956  SRF_RETURN_NEXT(funcctx, result);
1957  }
1958  /* do when there is no more left */
1959  else {
1960  pfree(quant2);
1961  SRF_RETURN_DONE(funcctx);
1962  }
1963 }
1964 
1968 PG_FUNCTION_INFO_V1(RASTER_quantileCoverage);
1969 Datum RASTER_quantileCoverage(PG_FUNCTION_ARGS)
1970 {
1971  FuncCallContext *funcctx;
1972  TupleDesc tupdesc;
1973 
1974  uint32_t i;
1975  rt_quantile covquant = NULL;
1976  rt_quantile covquant2;
1977  int call_cntr;
1978  int max_calls;
1979 
1980  POSTGIS_RT_DEBUG(3, "RASTER_quantileCoverage: Starting");
1981 
1982  /* first call of function */
1983  if (SRF_IS_FIRSTCALL()) {
1984  MemoryContext oldcontext;
1985 
1986  text *tablenametext = NULL;
1987  char *tablename = NULL;
1988  text *colnametext = NULL;
1989  char *colname = NULL;
1990  int32_t bandindex = 1;
1991  bool exclude_nodata_value = TRUE;
1992  double sample = 0;
1993  double *quantiles = NULL;
1994  uint32_t quantiles_count = 0;
1995  double quantile = 0;
1996  uint32_t count;
1997 
1998  int len = 0;
1999  char *sql = NULL;
2000  char *tmp = NULL;
2001  uint64_t cov_count = 0;
2002  int spi_result;
2003  Portal portal;
2004  SPITupleTable *tuptable = NULL;
2005  HeapTuple tuple;
2006  Datum datum;
2007  bool isNull = FALSE;
2008 
2009  rt_pgraster *pgraster = NULL;
2010  rt_raster raster = NULL;
2011  rt_band band = NULL;
2012  int num_bands = 0;
2013  struct quantile_llist *qlls = NULL;
2014  uint32_t qlls_count;
2015 
2016  int j;
2017  int n;
2018 
2019  ArrayType *array;
2020  Oid etype;
2021  Datum *e;
2022  bool *nulls;
2023  int16 typlen;
2024  bool typbyval;
2025  char typalign;
2026 
2027  POSTGIS_RT_DEBUG(3, "RASTER_quantileCoverage: first call of function");
2028 
2029  /* create a function context for cross-call persistence */
2030  funcctx = SRF_FIRSTCALL_INIT();
2031 
2032  /* switch to memory context appropriate for multiple function calls */
2033  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
2034 
2035  /* tablename is null, return null */
2036  if (PG_ARGISNULL(0)) {
2037  elog(NOTICE, "Table name must be provided");
2038  MemoryContextSwitchTo(oldcontext);
2039  SRF_RETURN_DONE(funcctx);
2040  }
2041  tablenametext = PG_GETARG_TEXT_P(0);
2042  tablename = text_to_cstring(tablenametext);
2043  if (!strlen(tablename)) {
2044  elog(NOTICE, "Table name must be provided");
2045  MemoryContextSwitchTo(oldcontext);
2046  SRF_RETURN_DONE(funcctx);
2047  }
2048  POSTGIS_RT_DEBUGF(3, "RASTER_quantileCoverage: tablename = %s", tablename);
2049 
2050  /* column name is null, return null */
2051  if (PG_ARGISNULL(1)) {
2052  elog(NOTICE, "Column name must be provided");
2053  MemoryContextSwitchTo(oldcontext);
2054  SRF_RETURN_DONE(funcctx);
2055  }
2056  colnametext = PG_GETARG_TEXT_P(1);
2057  colname = text_to_cstring(colnametext);
2058  if (!strlen(colname)) {
2059  elog(NOTICE, "Column name must be provided");
2060  MemoryContextSwitchTo(oldcontext);
2061  SRF_RETURN_DONE(funcctx);
2062  }
2063  POSTGIS_RT_DEBUGF(3, "RASTER_quantileCoverage: colname = %s", colname);
2064 
2065  /* band index is 1-based */
2066  if (!PG_ARGISNULL(2))
2067  bandindex = PG_GETARG_INT32(2);
2068 
2069  /* exclude_nodata_value flag */
2070  if (!PG_ARGISNULL(3))
2071  exclude_nodata_value = PG_GETARG_BOOL(3);
2072 
2073  /* sample % */
2074  if (!PG_ARGISNULL(4)) {
2075  sample = PG_GETARG_FLOAT8(4);
2076  if (sample < 0 || sample > 1) {
2077  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
2078  MemoryContextSwitchTo(oldcontext);
2079  SRF_RETURN_DONE(funcctx);
2080  }
2081  else if (FLT_EQ(sample, 0.0))
2082  sample = 1;
2083  }
2084  else
2085  sample = 1;
2086 
2087  /* quantiles */
2088  if (!PG_ARGISNULL(5)) {
2089  array = PG_GETARG_ARRAYTYPE_P(5);
2090  etype = ARR_ELEMTYPE(array);
2091  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
2092 
2093  switch (etype) {
2094  case FLOAT4OID:
2095  case FLOAT8OID:
2096  break;
2097  default:
2098  MemoryContextSwitchTo(oldcontext);
2099  elog(ERROR, "RASTER_quantileCoverage: Invalid data type for quantiles");
2100  SRF_RETURN_DONE(funcctx);
2101  break;
2102  }
2103 
2104  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
2105  &nulls, &n);
2106 
2107  quantiles = palloc(sizeof(double) * n);
2108  for (i = 0, j = 0; i < (uint32_t) n; i++) {
2109  if (nulls[i]) continue;
2110 
2111  switch (etype) {
2112  case FLOAT4OID:
2113  quantile = (double) DatumGetFloat4(e[i]);
2114  break;
2115  case FLOAT8OID:
2116  quantile = (double) DatumGetFloat8(e[i]);
2117  break;
2118  }
2119 
2120  if (quantile < 0 || quantile > 1) {
2121  elog(NOTICE, "Invalid value for quantile (must be between 0 and 1). Returning NULL");
2122  pfree(quantiles);
2123  MemoryContextSwitchTo(oldcontext);
2124  SRF_RETURN_DONE(funcctx);
2125  }
2126 
2127  quantiles[j] = quantile;
2128  POSTGIS_RT_DEBUGF(5, "quantiles[%d] = %f", j, quantiles[j]);
2129  j++;
2130  }
2131  quantiles_count = j;
2132 
2133  if (j < 1) {
2134  pfree(quantiles);
2135  quantiles = NULL;
2136  }
2137  }
2138 
2139  /* coverage stats */
2140  /* connect to database */
2141  spi_result = SPI_connect();
2142  if (spi_result != SPI_OK_CONNECT) {
2143  MemoryContextSwitchTo(oldcontext);
2144  elog(ERROR, "RASTER_quantileCoverage: Cannot connect to database using SPI");
2145  SRF_RETURN_DONE(funcctx);
2146  }
2147 
2148  len = sizeof(char) * (strlen("SELECT count FROM _st_summarystats('','',,::boolean,)") + strlen(tablename) + strlen(colname) + (MAX_INT_CHARLEN * 2) + MAX_DBL_CHARLEN + 1);
2149  sql = (char *) palloc(len);
2150  if (NULL == sql) {
2151 
2152  if (SPI_tuptable) SPI_freetuptable(tuptable);
2153  SPI_finish();
2154 
2155  MemoryContextSwitchTo(oldcontext);
2156  elog(ERROR, "RASTER_quantileCoverage: Cannot allocate memory for sql");
2157  SRF_RETURN_DONE(funcctx);
2158  }
2159 
2160  /* get stats */
2161  snprintf(sql, len, "SELECT count FROM _st_summarystats('%s','%s',%d,%d::boolean,%f)", tablename, colname, bandindex, (exclude_nodata_value ? 1 : 0), sample);
2162  POSTGIS_RT_DEBUGF(3, "stats sql: %s", sql);
2163  spi_result = SPI_execute(sql, TRUE, 0);
2164  pfree(sql);
2165  if (spi_result != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
2166 
2167  if (SPI_tuptable) SPI_freetuptable(tuptable);
2168  SPI_finish();
2169 
2170  MemoryContextSwitchTo(oldcontext);
2171  elog(ERROR, "RASTER_quantileCoverage: Cannot get summary stats of coverage");
2172  SRF_RETURN_DONE(funcctx);
2173  }
2174 
2175  tupdesc = SPI_tuptable->tupdesc;
2176  tuptable = SPI_tuptable;
2177  tuple = tuptable->vals[0];
2178 
2179  tmp = SPI_getvalue(tuple, tupdesc, 1);
2180  if (NULL == tmp || !strlen(tmp)) {
2181 
2182  if (SPI_tuptable) SPI_freetuptable(tuptable);
2183  SPI_finish();
2184 
2185  MemoryContextSwitchTo(oldcontext);
2186  elog(ERROR, "RASTER_quantileCoverage: Cannot get summary stats of coverage");
2187  SRF_RETURN_DONE(funcctx);
2188  }
2189  cov_count = strtol(tmp, NULL, 10);
2190  POSTGIS_RT_DEBUGF(3, "covcount = %d", (int) cov_count);
2191  pfree(tmp);
2192 
2193  /* iterate through rasters of coverage */
2194  /* create sql */
2195  len = sizeof(char) * (strlen("SELECT \"\" FROM \"\" WHERE \"\" IS NOT NULL") + (strlen(colname) * 2) + strlen(tablename) + 1);
2196  sql = (char *) palloc(len);
2197  if (NULL == sql) {
2198 
2199  if (SPI_tuptable) SPI_freetuptable(tuptable);
2200  SPI_finish();
2201 
2202  MemoryContextSwitchTo(oldcontext);
2203  elog(ERROR, "RASTER_quantileCoverage: Cannot allocate memory for sql");
2204  SRF_RETURN_DONE(funcctx);
2205  }
2206 
2207  /* get cursor */
2208  snprintf(sql, len, "SELECT \"%s\" FROM \"%s\" WHERE \"%s\" IS NOT NULL", colname, tablename, colname);
2209  POSTGIS_RT_DEBUGF(3, "coverage sql: %s", sql);
2210  portal = SPI_cursor_open_with_args(
2211  "coverage",
2212  sql,
2213  0, NULL,
2214  NULL, NULL,
2215  TRUE, 0
2216  );
2217  pfree(sql);
2218 
2219  /* process resultset */
2220  SPI_cursor_fetch(portal, TRUE, 1);
2221  while (SPI_processed == 1 && SPI_tuptable != NULL) {
2222  if (NULL != covquant) pfree(covquant);
2223 
2224  tupdesc = SPI_tuptable->tupdesc;
2225  tuple = SPI_tuptable->vals[0];
2226 
2227  datum = SPI_getbinval(tuple, tupdesc, 1, &isNull);
2228  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
2229  SPI_freetuptable(SPI_tuptable);
2230  SPI_cursor_close(portal);
2231  SPI_finish();
2232 
2233  MemoryContextSwitchTo(oldcontext);
2234  elog(ERROR, "RASTER_quantileCoverage: Cannot get raster of coverage");
2235  SRF_RETURN_DONE(funcctx);
2236  }
2237  else if (isNull) {
2238  SPI_cursor_fetch(portal, TRUE, 1);
2239  continue;
2240  }
2241 
2242  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(datum);
2243 
2244  raster = rt_raster_deserialize(pgraster, FALSE);
2245  if (!raster) {
2246 
2247  SPI_freetuptable(SPI_tuptable);
2248  SPI_cursor_close(portal);
2249  SPI_finish();
2250 
2251  MemoryContextSwitchTo(oldcontext);
2252  elog(ERROR, "RASTER_quantileCoverage: Cannot deserialize raster");
2253  SRF_RETURN_DONE(funcctx);
2254  }
2255 
2256  /* inspect number of bands*/
2257  num_bands = rt_raster_get_num_bands(raster);
2258  if (bandindex < 1 || bandindex > num_bands) {
2259  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
2260 
2261  rt_raster_destroy(raster);
2262 
2263  SPI_freetuptable(SPI_tuptable);
2264  SPI_cursor_close(portal);
2265  SPI_finish();
2266 
2267  MemoryContextSwitchTo(oldcontext);
2268  SRF_RETURN_DONE(funcctx);
2269  }
2270 
2271  /* get band */
2272  band = rt_raster_get_band(raster, bandindex - 1);
2273  if (!band) {
2274  elog(NOTICE, "Cannot find raster band of index %d. Returning NULL", bandindex);
2275 
2276  rt_raster_destroy(raster);
2277 
2278  SPI_freetuptable(SPI_tuptable);
2279  SPI_cursor_close(portal);
2280  SPI_finish();
2281 
2282  MemoryContextSwitchTo(oldcontext);
2283  SRF_RETURN_DONE(funcctx);
2284  }
2285 
2286  covquant = rt_band_get_quantiles_stream(
2287  band,
2288  exclude_nodata_value, sample, cov_count,
2289  &qlls, &qlls_count,
2290  quantiles, quantiles_count,
2291  &count
2292  );
2293 
2294  rt_band_destroy(band);
2295  rt_raster_destroy(raster);
2296 
2297  if (!covquant || !count) {
2298  elog(NOTICE, "Cannot compute quantiles for band at index %d", bandindex);
2299 
2300  SPI_freetuptable(SPI_tuptable);
2301  SPI_cursor_close(portal);
2302  SPI_finish();
2303 
2304  MemoryContextSwitchTo(oldcontext);
2305  SRF_RETURN_DONE(funcctx);
2306  }
2307 
2308  /* next record */
2309  SPI_cursor_fetch(portal, TRUE, 1);
2310  }
2311 
2312  covquant2 = SPI_palloc(sizeof(struct rt_quantile_t) * count);
2313  for (i = 0; i < count; i++) {
2314  covquant2[i].quantile = covquant[i].quantile;
2315  covquant2[i].has_value = covquant[i].has_value;
2316  if (covquant2[i].has_value)
2317  covquant2[i].value = covquant[i].value;
2318  }
2319 
2320  pfree(covquant);
2321  quantile_llist_destroy(&qlls, qlls_count);
2322 
2323  if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
2324  SPI_cursor_close(portal);
2325  SPI_finish();
2326 
2327  if (quantiles_count) pfree(quantiles);
2328 
2329  POSTGIS_RT_DEBUGF(3, "%d quantiles returned", count);
2330 
2331  /* Store needed information */
2332  funcctx->user_fctx = covquant2;
2333 
2334  /* total number of tuples to be returned */
2335  funcctx->max_calls = count;
2336 
2337  /* Build a tuple descriptor for our result type */
2338  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
2339  ereport(ERROR, (
2340  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
2341  errmsg(
2342  "function returning record called in context "
2343  "that cannot accept type record"
2344  )
2345  ));
2346  }
2347 
2348  BlessTupleDesc(tupdesc);
2349  funcctx->tuple_desc = tupdesc;
2350 
2351  MemoryContextSwitchTo(oldcontext);
2352  }
2353 
2354  /* stuff done on every call of the function */
2355  funcctx = SRF_PERCALL_SETUP();
2356 
2357  call_cntr = funcctx->call_cntr;
2358  max_calls = funcctx->max_calls;
2359  tupdesc = funcctx->tuple_desc;
2360  covquant2 = funcctx->user_fctx;
2361 
2362  /* do when there is more left to send */
2363  if (call_cntr < max_calls) {
2364  Datum values[VALUES_LENGTH];
2365  bool nulls[VALUES_LENGTH];
2366  HeapTuple tuple;
2367  Datum result;
2368 
2369  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
2370 
2371  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
2372 
2373  values[0] = Float8GetDatum(covquant2[call_cntr].quantile);
2374  if (covquant2[call_cntr].has_value)
2375  values[1] = Float8GetDatum(covquant2[call_cntr].value);
2376  else
2377  nulls[1] = TRUE;
2378 
2379  /* build a tuple */
2380  tuple = heap_form_tuple(tupdesc, values, nulls);
2381 
2382  /* make the tuple into a datum */
2383  result = HeapTupleGetDatum(tuple);
2384 
2385  SRF_RETURN_NEXT(funcctx, result);
2386  }
2387  /* do when there is no more left */
2388  else {
2389  POSTGIS_RT_DEBUG(3, "done");
2390  pfree(covquant2);
2391  SRF_RETURN_DONE(funcctx);
2392  }
2393 }
2394 
2395 #undef VALUES_LENGTH
2396 #define VALUES_LENGTH 3
2397 
2398 /* get counts of values */
2399 PG_FUNCTION_INFO_V1(RASTER_valueCount);
2400 Datum RASTER_valueCount(PG_FUNCTION_ARGS) {
2401  FuncCallContext *funcctx;
2402  TupleDesc tupdesc;
2403 
2404  int i;
2405  rt_valuecount vcnts;
2406  rt_valuecount vcnts2;
2407  int call_cntr;
2408  int max_calls;
2409 
2410  /* first call of function */
2411  if (SRF_IS_FIRSTCALL()) {
2412  MemoryContext oldcontext;
2413 
2414  rt_pgraster *pgraster = NULL;
2415  rt_raster raster = NULL;
2416  rt_band band = NULL;
2417  int32_t bandindex = 0;
2418  int num_bands = 0;
2419  bool exclude_nodata_value = TRUE;
2420  double *search_values = NULL;
2421  uint32_t search_values_count = 0;
2422  double roundto = 0;
2423  uint32_t count;
2424 
2425  int j;
2426  int n;
2427 
2428  ArrayType *array;
2429  Oid etype;
2430  Datum *e;
2431  bool *nulls;
2432  int16 typlen;
2433  bool typbyval;
2434  char typalign;
2435 
2436  /* create a function context for cross-call persistence */
2437  funcctx = SRF_FIRSTCALL_INIT();
2438 
2439  /* switch to memory context appropriate for multiple function calls */
2440  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
2441 
2442  /* pgraster is null, return nothing */
2443  if (PG_ARGISNULL(0)) {
2444  MemoryContextSwitchTo(oldcontext);
2445  SRF_RETURN_DONE(funcctx);
2446  }
2447  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2448 
2449  raster = rt_raster_deserialize(pgraster, FALSE);
2450  if (!raster) {
2451  PG_FREE_IF_COPY(pgraster, 0);
2452  MemoryContextSwitchTo(oldcontext);
2453  elog(ERROR, "RASTER_valueCount: Cannot deserialize raster");
2454  SRF_RETURN_DONE(funcctx);
2455  }
2456 
2457  /* band index is 1-based */
2458  bandindex = PG_GETARG_INT32(1);
2459  num_bands = rt_raster_get_num_bands(raster);
2460  if (bandindex < 1 || bandindex > num_bands) {
2461  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
2462  rt_raster_destroy(raster);
2463  PG_FREE_IF_COPY(pgraster, 0);
2464  MemoryContextSwitchTo(oldcontext);
2465  SRF_RETURN_DONE(funcctx);
2466  }
2467 
2468  /* exclude_nodata_value flag */
2469  if (!PG_ARGISNULL(2))
2470  exclude_nodata_value = PG_GETARG_BOOL(2);
2471 
2472  /* search values */
2473  if (!PG_ARGISNULL(3)) {
2474  array = PG_GETARG_ARRAYTYPE_P(3);
2475  etype = ARR_ELEMTYPE(array);
2476  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
2477 
2478  switch (etype) {
2479  case FLOAT4OID:
2480  case FLOAT8OID:
2481  break;
2482  default:
2483  rt_raster_destroy(raster);
2484  PG_FREE_IF_COPY(pgraster, 0);
2485  MemoryContextSwitchTo(oldcontext);
2486  elog(ERROR, "RASTER_valueCount: Invalid data type for values");
2487  SRF_RETURN_DONE(funcctx);
2488  break;
2489  }
2490 
2491  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
2492  &nulls, &n);
2493 
2494  search_values = palloc(sizeof(double) * n);
2495  for (i = 0, j = 0; i < n; i++) {
2496  if (nulls[i]) continue;
2497 
2498  switch (etype) {
2499  case FLOAT4OID:
2500  search_values[j] = (double) DatumGetFloat4(e[i]);
2501  break;
2502  case FLOAT8OID:
2503  search_values[j] = (double) DatumGetFloat8(e[i]);
2504  break;
2505  }
2506 
2507  POSTGIS_RT_DEBUGF(5, "search_values[%d] = %f", j, search_values[j]);
2508  j++;
2509  }
2510  search_values_count = j;
2511 
2512  if (j < 1) {
2513  pfree(search_values);
2514  search_values = NULL;
2515  }
2516  }
2517 
2518  /* roundto */
2519  if (!PG_ARGISNULL(4)) {
2520  roundto = PG_GETARG_FLOAT8(4);
2521  if (roundto < 0.) roundto = 0;
2522  }
2523 
2524  /* get band */
2525  band = rt_raster_get_band(raster, bandindex - 1);
2526  if (!band) {
2527  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
2528  rt_raster_destroy(raster);
2529  PG_FREE_IF_COPY(pgraster, 0);
2530  MemoryContextSwitchTo(oldcontext);
2531  SRF_RETURN_DONE(funcctx);
2532  }
2533 
2534  /* get counts of values */
2535  vcnts = rt_band_get_value_count(band, (int) exclude_nodata_value, search_values, search_values_count, roundto, NULL, &count);
2536  rt_band_destroy(band);
2537  rt_raster_destroy(raster);
2538  PG_FREE_IF_COPY(pgraster, 0);
2539  if (NULL == vcnts || !count) {
2540  elog(NOTICE, "Cannot count the values for band at index %d", bandindex);
2541  MemoryContextSwitchTo(oldcontext);
2542  SRF_RETURN_DONE(funcctx);
2543  }
2544 
2545  POSTGIS_RT_DEBUGF(3, "%d value counts returned", count);
2546 
2547  /* Store needed information */
2548  funcctx->user_fctx = vcnts;
2549 
2550  /* total number of tuples to be returned */
2551  funcctx->max_calls = count;
2552 
2553  /* Build a tuple descriptor for our result type */
2554  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
2555  ereport(ERROR, (
2556  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
2557  errmsg(
2558  "function returning record called in context "
2559  "that cannot accept type record"
2560  )
2561  ));
2562  }
2563 
2564  BlessTupleDesc(tupdesc);
2565  funcctx->tuple_desc = tupdesc;
2566 
2567  MemoryContextSwitchTo(oldcontext);
2568  }
2569 
2570  /* stuff done on every call of the function */
2571  funcctx = SRF_PERCALL_SETUP();
2572 
2573  call_cntr = funcctx->call_cntr;
2574  max_calls = funcctx->max_calls;
2575  tupdesc = funcctx->tuple_desc;
2576  vcnts2 = funcctx->user_fctx;
2577 
2578  /* do when there is more left to send */
2579  if (call_cntr < max_calls) {
2580  Datum values[VALUES_LENGTH];
2581  bool nulls[VALUES_LENGTH];
2582  HeapTuple tuple;
2583  Datum result;
2584 
2585  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
2586 
2587  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
2588 
2589  values[0] = Float8GetDatum(vcnts2[call_cntr].value);
2590  values[1] = UInt32GetDatum(vcnts2[call_cntr].count);
2591  values[2] = Float8GetDatum(vcnts2[call_cntr].percent);
2592 
2593  /* build a tuple */
2594  tuple = heap_form_tuple(tupdesc, values, nulls);
2595 
2596  /* make the tuple into a datum */
2597  result = HeapTupleGetDatum(tuple);
2598 
2599  SRF_RETURN_NEXT(funcctx, result);
2600  }
2601  /* do when there is no more left */
2602  else {
2603  pfree(vcnts2);
2604  SRF_RETURN_DONE(funcctx);
2605  }
2606 }
2607 
2608 /* get counts of values for a coverage */
2609 PG_FUNCTION_INFO_V1(RASTER_valueCountCoverage);
2610 Datum RASTER_valueCountCoverage(PG_FUNCTION_ARGS) {
2611  FuncCallContext *funcctx;
2612  TupleDesc tupdesc;
2613 
2614  uint32_t i;
2615  uint64_t covcount = 0;
2616  uint64_t covtotal = 0;
2617  rt_valuecount covvcnts = NULL;
2618  rt_valuecount covvcnts2;
2619  int call_cntr;
2620  int max_calls;
2621 
2622  POSTGIS_RT_DEBUG(3, "RASTER_valueCountCoverage: Starting");
2623 
2624  /* first call of function */
2625  if (SRF_IS_FIRSTCALL()) {
2626  MemoryContext oldcontext;
2627 
2628  text *tablenametext = NULL;
2629  char *tablename = NULL;
2630  text *colnametext = NULL;
2631  char *colname = NULL;
2632  int32_t bandindex = 1;
2633  bool exclude_nodata_value = TRUE;
2634  double *search_values = NULL;
2635  uint32_t search_values_count = 0;
2636  double roundto = 0;
2637 
2638  int len = 0;
2639  char *sql = NULL;
2640  int spi_result;
2641  Portal portal;
2642  HeapTuple tuple;
2643  Datum datum;
2644  bool isNull = FALSE;
2645  rt_pgraster *pgraster = NULL;
2646  rt_raster raster = NULL;
2647  rt_band band = NULL;
2648  int num_bands = 0;
2649  uint32_t count;
2650  uint32_t total;
2651  rt_valuecount vcnts = NULL;
2652  int exists = 0;
2653 
2654  uint32_t j;
2655  int n;
2656 
2657  ArrayType *array;
2658  Oid etype;
2659  Datum *e;
2660  bool *nulls;
2661  int16 typlen;
2662  bool typbyval;
2663  char typalign;
2664 
2665  /* create a function context for cross-call persistence */
2666  funcctx = SRF_FIRSTCALL_INIT();
2667 
2668  /* switch to memory context appropriate for multiple function calls */
2669  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
2670 
2671  /* tablename is null, return null */
2672  if (PG_ARGISNULL(0)) {
2673  elog(NOTICE, "Table name must be provided");
2674  MemoryContextSwitchTo(oldcontext);
2675  SRF_RETURN_DONE(funcctx);
2676  }
2677  tablenametext = PG_GETARG_TEXT_P(0);
2678  tablename = text_to_cstring(tablenametext);
2679  if (!strlen(tablename)) {
2680  elog(NOTICE, "Table name must be provided");
2681  MemoryContextSwitchTo(oldcontext);
2682  SRF_RETURN_DONE(funcctx);
2683  }
2684  POSTGIS_RT_DEBUGF(3, "tablename = %s", tablename);
2685 
2686  /* column name is null, return null */
2687  if (PG_ARGISNULL(1)) {
2688  elog(NOTICE, "Column name must be provided");
2689  MemoryContextSwitchTo(oldcontext);
2690  SRF_RETURN_DONE(funcctx);
2691  }
2692  colnametext = PG_GETARG_TEXT_P(1);
2693  colname = text_to_cstring(colnametext);
2694  if (!strlen(colname)) {
2695  elog(NOTICE, "Column name must be provided");
2696  MemoryContextSwitchTo(oldcontext);
2697  SRF_RETURN_DONE(funcctx);
2698  }
2699  POSTGIS_RT_DEBUGF(3, "colname = %s", colname);
2700 
2701  /* band index is 1-based */
2702  if (!PG_ARGISNULL(2))
2703  bandindex = PG_GETARG_INT32(2);
2704 
2705  /* exclude_nodata_value flag */
2706  if (!PG_ARGISNULL(3))
2707  exclude_nodata_value = PG_GETARG_BOOL(3);
2708 
2709  /* search values */
2710  if (!PG_ARGISNULL(4)) {
2711  array = PG_GETARG_ARRAYTYPE_P(4);
2712  etype = ARR_ELEMTYPE(array);
2713  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
2714 
2715  switch (etype) {
2716  case FLOAT4OID:
2717  case FLOAT8OID:
2718  break;
2719  default:
2720  MemoryContextSwitchTo(oldcontext);
2721  elog(ERROR, "RASTER_valueCountCoverage: Invalid data type for values");
2722  SRF_RETURN_DONE(funcctx);
2723  break;
2724  }
2725 
2726  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
2727  &nulls, &n);
2728 
2729  search_values = palloc(sizeof(double) * n);
2730  for (i = 0, j = 0; i < (uint32_t) n; i++) {
2731  if (nulls[i]) continue;
2732 
2733  switch (etype) {
2734  case FLOAT4OID:
2735  search_values[j] = (double) DatumGetFloat4(e[i]);
2736  break;
2737  case FLOAT8OID:
2738  search_values[j] = (double) DatumGetFloat8(e[i]);
2739  break;
2740  }
2741 
2742  POSTGIS_RT_DEBUGF(5, "search_values[%d] = %f", j, search_values[j]);
2743  j++;
2744  }
2745  search_values_count = j;
2746 
2747  if (j < 1) {
2748  pfree(search_values);
2749  search_values = NULL;
2750  }
2751  }
2752 
2753  /* roundto */
2754  if (!PG_ARGISNULL(5)) {
2755  roundto = PG_GETARG_FLOAT8(5);
2756  if (roundto < 0.) roundto = 0;
2757  }
2758 
2759  /* iterate through rasters of coverage */
2760  /* connect to database */
2761  spi_result = SPI_connect();
2762  if (spi_result != SPI_OK_CONNECT) {
2763 
2764  if (search_values_count) pfree(search_values);
2765 
2766  MemoryContextSwitchTo(oldcontext);
2767  elog(ERROR, "RASTER_valueCountCoverage: Cannot connect to database using SPI");
2768  SRF_RETURN_DONE(funcctx);
2769  }
2770 
2771  /* create sql */
2772  len = sizeof(char) * (strlen("SELECT \"\" FROM \"\" WHERE \"\" IS NOT NULL") + (strlen(colname) * 2) + strlen(tablename) + 1);
2773  sql = (char *) palloc(len);
2774  if (NULL == sql) {
2775 
2776  if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
2777  SPI_finish();
2778 
2779  if (search_values_count) pfree(search_values);
2780 
2781  MemoryContextSwitchTo(oldcontext);
2782  elog(ERROR, "RASTER_valueCountCoverage: Cannot allocate memory for sql");
2783  SRF_RETURN_DONE(funcctx);
2784  }
2785 
2786  /* get cursor */
2787  snprintf(sql, len, "SELECT \"%s\" FROM \"%s\" WHERE \"%s\" IS NOT NULL", colname, tablename, colname);
2788  POSTGIS_RT_DEBUGF(3, "RASTER_valueCountCoverage: %s", sql);
2789  portal = SPI_cursor_open_with_args(
2790  "coverage",
2791  sql,
2792  0, NULL,
2793  NULL, NULL,
2794  TRUE, 0
2795  );
2796  pfree(sql);
2797 
2798  /* process resultset */
2799  SPI_cursor_fetch(portal, TRUE, 1);
2800  while (SPI_processed == 1 && SPI_tuptable != NULL) {
2801  tupdesc = SPI_tuptable->tupdesc;
2802  tuple = SPI_tuptable->vals[0];
2803 
2804  datum = SPI_getbinval(tuple, tupdesc, 1, &isNull);
2805  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
2806 
2807  SPI_freetuptable(SPI_tuptable);
2808  SPI_cursor_close(portal);
2809  SPI_finish();
2810 
2811  if (NULL != covvcnts) pfree(covvcnts);
2812  if (search_values_count) pfree(search_values);
2813 
2814  MemoryContextSwitchTo(oldcontext);
2815  elog(ERROR, "RASTER_valueCountCoverage: Cannot get raster of coverage");
2816  SRF_RETURN_DONE(funcctx);
2817  }
2818  else if (isNull) {
2819  SPI_cursor_fetch(portal, TRUE, 1);
2820  continue;
2821  }
2822 
2823  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(datum);
2824 
2825  raster = rt_raster_deserialize(pgraster, FALSE);
2826  if (!raster) {
2827 
2828  SPI_freetuptable(SPI_tuptable);
2829  SPI_cursor_close(portal);
2830  SPI_finish();
2831 
2832  if (NULL != covvcnts) pfree(covvcnts);
2833  if (search_values_count) pfree(search_values);
2834 
2835  MemoryContextSwitchTo(oldcontext);
2836  elog(ERROR, "RASTER_valueCountCoverage: Cannot deserialize raster");
2837  SRF_RETURN_DONE(funcctx);
2838  }
2839 
2840  /* inspect number of bands*/
2841  num_bands = rt_raster_get_num_bands(raster);
2842  if (bandindex < 1 || bandindex > num_bands) {
2843  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
2844 
2845  rt_raster_destroy(raster);
2846 
2847  SPI_freetuptable(SPI_tuptable);
2848  SPI_cursor_close(portal);
2849  SPI_finish();
2850 
2851  if (NULL != covvcnts) pfree(covvcnts);
2852  if (search_values_count) pfree(search_values);
2853 
2854  MemoryContextSwitchTo(oldcontext);
2855  SRF_RETURN_DONE(funcctx);
2856  }
2857 
2858  /* get band */
2859  band = rt_raster_get_band(raster, bandindex - 1);
2860  if (!band) {
2861  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
2862 
2863  rt_raster_destroy(raster);
2864 
2865  SPI_freetuptable(SPI_tuptable);
2866  SPI_cursor_close(portal);
2867  SPI_finish();
2868 
2869  if (NULL != covvcnts) pfree(covvcnts);
2870  if (search_values_count) pfree(search_values);
2871 
2872  MemoryContextSwitchTo(oldcontext);
2873  SRF_RETURN_DONE(funcctx);
2874  }
2875 
2876  /* get counts of values */
2877  vcnts = rt_band_get_value_count(band, (int) exclude_nodata_value, search_values, search_values_count, roundto, &total, &count);
2878  rt_band_destroy(band);
2879  rt_raster_destroy(raster);
2880  if (NULL == vcnts || !count) {
2881  elog(NOTICE, "Cannot count the values for band at index %d", bandindex);
2882 
2883  SPI_freetuptable(SPI_tuptable);
2884  SPI_cursor_close(portal);
2885  SPI_finish();
2886 
2887  if (NULL != covvcnts) free(covvcnts);
2888  if (search_values_count) pfree(search_values);
2889 
2890  MemoryContextSwitchTo(oldcontext);
2891  SRF_RETURN_DONE(funcctx);
2892  }
2893 
2894  POSTGIS_RT_DEBUGF(3, "%d value counts returned", count);
2895 
2896  if (NULL == covvcnts) {
2897  covvcnts = (rt_valuecount) SPI_palloc(sizeof(struct rt_valuecount_t) * count);
2898  if (NULL == covvcnts) {
2899 
2900  SPI_freetuptable(SPI_tuptable);
2901  SPI_cursor_close(portal);
2902  SPI_finish();
2903 
2904  if (search_values_count) pfree(search_values);
2905 
2906  MemoryContextSwitchTo(oldcontext);
2907  elog(ERROR, "RASTER_valueCountCoverage: Cannot allocate memory for value counts of coverage");
2908  SRF_RETURN_DONE(funcctx);
2909  }
2910 
2911  for (i = 0; i < count; i++) {
2912  covvcnts[i].value = vcnts[i].value;
2913  covvcnts[i].count = vcnts[i].count;
2914  covvcnts[i].percent = -1;
2915  }
2916 
2917  covcount = count;
2918  }
2919  else {
2920  for (i = 0; i < count; i++) {
2921  exists = 0;
2922 
2923  for (j = 0; j < covcount; j++) {
2924  if (FLT_EQ(vcnts[i].value, covvcnts[j].value)) {
2925  exists = 1;
2926  break;
2927  }
2928  }
2929 
2930  if (exists) {
2931  covvcnts[j].count += vcnts[i].count;
2932  }
2933  else {
2934  covcount++;
2935  covvcnts = SPI_repalloc(covvcnts, sizeof(struct rt_valuecount_t) * covcount);
2936  if (!covvcnts) {
2937  SPI_freetuptable(SPI_tuptable);
2938  SPI_cursor_close(portal);
2939  SPI_finish();
2940 
2941  if (search_values_count) pfree(search_values);
2942 
2943  MemoryContextSwitchTo(oldcontext);
2944  elog(ERROR, "RASTER_valueCountCoverage: Cannot change allocated memory for value counts of coverage");
2945  SRF_RETURN_DONE(funcctx);
2946  }
2947 
2948  covvcnts[covcount - 1].value = vcnts[i].value;
2949  covvcnts[covcount - 1].count = vcnts[i].count;
2950  covvcnts[covcount - 1].percent = -1;
2951  }
2952  }
2953  }
2954 
2955  covtotal += total;
2956 
2957  pfree(vcnts);
2958 
2959  /* next record */
2960  SPI_cursor_fetch(portal, TRUE, 1);
2961  }
2962 
2963  if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
2964  SPI_cursor_close(portal);
2965  SPI_finish();
2966 
2967  if (search_values_count) pfree(search_values);
2968 
2969  /* compute percentages */
2970  for (i = 0; i < covcount; i++) {
2971  covvcnts[i].percent = (double) covvcnts[i].count / covtotal;
2972  }
2973 
2974  /* Store needed information */
2975  funcctx->user_fctx = covvcnts;
2976 
2977  /* total number of tuples to be returned */
2978  funcctx->max_calls = covcount;
2979 
2980  /* Build a tuple descriptor for our result type */
2981  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
2982  ereport(ERROR, (
2983  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
2984  errmsg(
2985  "function returning record called in context "
2986  "that cannot accept type record"
2987  )
2988  ));
2989  }
2990 
2991  BlessTupleDesc(tupdesc);
2992  funcctx->tuple_desc = tupdesc;
2993 
2994  MemoryContextSwitchTo(oldcontext);
2995  }
2996 
2997  /* stuff done on every call of the function */
2998  funcctx = SRF_PERCALL_SETUP();
2999 
3000  call_cntr = funcctx->call_cntr;
3001  max_calls = funcctx->max_calls;
3002  tupdesc = funcctx->tuple_desc;
3003  covvcnts2 = funcctx->user_fctx;
3004 
3005  /* do when there is more left to send */
3006  if (call_cntr < max_calls) {
3007  Datum values[VALUES_LENGTH];
3008  bool nulls[VALUES_LENGTH];
3009  HeapTuple tuple;
3010  Datum result;
3011 
3012  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
3013 
3014  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
3015 
3016  values[0] = Float8GetDatum(covvcnts2[call_cntr].value);
3017  values[1] = UInt32GetDatum(covvcnts2[call_cntr].count);
3018  values[2] = Float8GetDatum(covvcnts2[call_cntr].percent);
3019 
3020  /* build a tuple */
3021  tuple = heap_form_tuple(tupdesc, values, nulls);
3022 
3023  /* make the tuple into a datum */
3024  result = HeapTupleGetDatum(tuple);
3025 
3026  SRF_RETURN_NEXT(funcctx, result);
3027  }
3028  /* do when there is no more left */
3029  else {
3030  pfree(covvcnts2);
3031  SRF_RETURN_DONE(funcctx);
3032  }
3033 }
3034 
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:2234
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:61
genraster.count
int count
Definition: genraster.py:56
ovdump.band
band
Definition: ovdump.py:57
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:917
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:1191
RASTER_valueCountCoverage
Datum RASTER_valueCountCoverage(PG_FUNCTION_ARGS)
Definition: rtpg_statistics.c:2610
RASTER_valueCount
Datum RASTER_valueCount(PG_FUNCTION_ARGS)
Definition: rtpg_statistics.c:2400
RASTER_quantileCoverage
Datum RASTER_quantileCoverage(PG_FUNCTION_ARGS)
Definition: rtpg_statistics.c:1969
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:2396
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:1720
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:2394
quantile_llist
Definition: librtcore.h:2392
rt_band_t
Definition: librtcore.h:2313
rt_bandstats_t::mean
double mean
Definition: librtcore.h:2365
rt_bandstats_t::sorted
int sorted
Definition: librtcore.h:2369
rt_bandstats_t::sample
double sample
Definition: librtcore.h:2359
rt_bandstats_t::sum
double sum
Definition: librtcore.h:2364
rt_bandstats_t::max
double max
Definition: librtcore.h:2363
rt_bandstats_t::count
uint32_t count
Definition: librtcore.h:2360
rt_bandstats_t::min
double min
Definition: librtcore.h:2362
rt_bandstats_t::values
double * values
Definition: librtcore.h:2368
rt_bandstats_t::stddev
double stddev
Definition: librtcore.h:2366
rt_bandstats_t
Definition: librtcore.h:2358
rt_histogram_t::count
uint32_t count
Definition: librtcore.h:2374
rt_histogram_t::percent
double percent
Definition: librtcore.h:2375
rt_histogram_t::min
double min
Definition: librtcore.h:2377
rt_histogram_t::max
double max
Definition: librtcore.h:2378
rt_histogram_t
Definition: librtcore.h:2373
rt_quantile_t::value
double value
Definition: librtcore.h:2387
rt_quantile_t::quantile
double quantile
Definition: librtcore.h:2386
rt_quantile_t::has_value
uint32_t has_value
Definition: librtcore.h:2388
rt_quantile_t
Definition: librtcore.h:2385
rt_raster_serialized_t
Struct definitions.
Definition: librtcore.h:2250
rt_raster_t
Definition: librtcore.h:2284
rt_valuecount_t::count
uint32_t count
Definition: librtcore.h:2425
rt_valuecount_t::percent
double percent
Definition: librtcore.h:2426
rt_valuecount_t::value
double value
Definition: librtcore.h:2424
rt_valuecount_t
Definition: librtcore.h:2423
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
uint32_t
unsigned int uint32_t
Definition: uthash.h:78