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