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