PostGIS 3.7.0dev-r@@SVN_REVISION@@
Loading...
Searching...
No Matches
rtpg_statistics.c
Go to the documentation of this file.
1/*
2 *
3 * WKTRaster - Raster Types for PostGIS
4 * http://trac.osgeo.org/postgis/wiki/WKTRaster
5 *
6 * Copyright (C) 2013 Bborie Park <dustymugs@gmail.com>
7 * Copyright (C) 2011-2013 Regents of the University of California
8 * <bkpark@ucdavis.edu>
9 * Copyright (C) 2010-2011 Jorge Arevalo <jorge.arevalo@deimos-space.com>
10 * Copyright (C) 2010-2011 David Zwarg <dzwarg@azavea.com>
11 * Copyright (C) 2009-2011 Pierre Racine <pierre.racine@sbf.ulaval.ca>
12 * Copyright (C) 2009-2011 Mateusz Loskot <mateusz@loskot.net>
13 * Copyright (C) 2008-2009 Sandro Santilli <strk@kbt.io>
14 *
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
19 *
20 * This program is distributed in the hope that it will be useful,
21 * but WITHOUT ANY WARRANTY; without even the implied warranty of
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 * GNU General Public License for more details.
24 *
25 * You should have received a copy of the GNU General Public License
26 * along with this program; if not, write to the Free Software Foundation,
27 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
28 *
29 */
30
31#include <postgres.h>
32#include <fmgr.h>
33#include <utils/builtins.h> /* for text_to_cstring() */
34#include "utils/lsyscache.h" /* for get_typlenbyvalalign */
35#include "utils/array.h" /* for ArrayType */
36#include "catalog/pg_type.h" /* for INT2OID, INT4OID, FLOAT4OID, FLOAT8OID and TEXTOID */
37#include <executor/spi.h>
38#include <funcapi.h> /* for SRF */
39
40#include "../../postgis_config.h"
41
42
43#include "access/htup_details.h" /* for heap_form_tuple() */
44
45
46#include "rtpostgis.h"
47
48/* Get summary stats */
49Datum RASTER_summaryStats(PG_FUNCTION_ARGS);
50Datum RASTER_summaryStatsCoverage(PG_FUNCTION_ARGS);
51
52Datum RASTER_summaryStats_transfn(PG_FUNCTION_ARGS);
53Datum RASTER_summaryStats_finalfn(PG_FUNCTION_ARGS);
54
55/* get histogram */
56Datum RASTER_histogram(PG_FUNCTION_ARGS);
57
58/* get quantiles */
59Datum RASTER_quantile(PG_FUNCTION_ARGS);
60
61/* get counts of values */
62Datum RASTER_valueCount(PG_FUNCTION_ARGS);
63Datum RASTER_valueCountCoverage(PG_FUNCTION_ARGS);
64
65#define VALUES_LENGTH 6
66
70PG_FUNCTION_INFO_V1(RASTER_summaryStats);
71Datum 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 Datum values[VALUES_LENGTH];
84 bool nulls[VALUES_LENGTH];
85 HeapTuple tuple;
86 Datum result;
87
88 /* pgraster is null, return null */
89 if (PG_ARGISNULL(0))
90 PG_RETURN_NULL();
91 pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
92
93 raster = rt_raster_deserialize(pgraster, FALSE);
94 if (!raster) {
95 PG_FREE_IF_COPY(pgraster, 0);
96 elog(ERROR, "RASTER_summaryStats: Cannot deserialize raster");
97 PG_RETURN_NULL();
98 }
99
100 /* band index is 1-based */
101 if (!PG_ARGISNULL(1))
102 bandindex = PG_GETARG_INT32(1);
103 num_bands = rt_raster_get_num_bands(raster);
104 if (bandindex < 1 || bandindex > num_bands) {
105 elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
106 rt_raster_destroy(raster);
107 PG_FREE_IF_COPY(pgraster, 0);
108 PG_RETURN_NULL();
109 }
110
111 /* exclude_nodata_value flag */
112 if (!PG_ARGISNULL(2))
113 exclude_nodata_value = PG_GETARG_BOOL(2);
114
115 /* sample % */
116 if (!PG_ARGISNULL(3)) {
117 sample = PG_GETARG_FLOAT8(3);
118 if (sample < 0 || sample > 1) {
119 elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
120 rt_raster_destroy(raster);
121 PG_FREE_IF_COPY(pgraster, 0);
122 PG_RETURN_NULL();
123 }
124 else if (FLT_EQ(sample, 0.0))
125 sample = 1;
126 }
127 else
128 sample = 1;
129
130 /* get band */
131 band = rt_raster_get_band(raster, bandindex - 1);
132 if (!band) {
133 elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
134 rt_raster_destroy(raster);
135 PG_FREE_IF_COPY(pgraster, 0);
136 PG_RETURN_NULL();
137 }
138
139 /* we don't need the raw values, hence the zero parameter */
140 stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 0, NULL, NULL, NULL);
141 rt_band_destroy(band);
142 rt_raster_destroy(raster);
143 PG_FREE_IF_COPY(pgraster, 0);
144 if (NULL == stats) {
145 elog(NOTICE, "Cannot compute summary statistics for band at index %d. Returning NULL", bandindex);
146 PG_RETURN_NULL();
147 }
148
149 /* Build a tuple descriptor for our result type */
150 if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
151 ereport(ERROR, (
152 errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
153 errmsg(
154 "function returning record called in context "
155 "that cannot accept type record"
156 )
157 ));
158 }
159
160 BlessTupleDesc(tupdesc);
161
162 memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
163
164 values[0] = Int64GetDatum(stats->count);
165 if (stats->count > 0) {
166 values[1] = Float8GetDatum(stats->sum);
167 values[2] = Float8GetDatum(stats->mean);
168 values[3] = Float8GetDatum(stats->stddev);
169 values[4] = Float8GetDatum(stats->min);
170 values[5] = Float8GetDatum(stats->max);
171 }
172 else {
173 nulls[1] = TRUE;
174 nulls[2] = TRUE;
175 nulls[3] = TRUE;
176 nulls[4] = TRUE;
177 nulls[5] = TRUE;
178 }
179
180 /* build a tuple */
181 tuple = heap_form_tuple(tupdesc, values, nulls);
182
183 /* make the tuple into a datum */
184 result = HeapTupleGetDatum(tuple);
185
186 /* clean up */
187 pfree(stats);
188
189 PG_RETURN_DATUM(result);
190}
191
195PG_FUNCTION_INFO_V1(RASTER_summaryStatsCoverage);
196Datum RASTER_summaryStatsCoverage(PG_FUNCTION_ARGS)
197{
198 text *tablenametext = NULL;
199 char *tablename = NULL;
200 text *colnametext = NULL;
201 char *colname = NULL;
202 int32_t bandindex = 1;
203 bool exclude_nodata_value = TRUE;
204 double sample = 0;
205
206 int len = 0;
207 char *sql = NULL;
208 int spi_result;
209 Portal portal;
210 TupleDesc tupdesc;
211 SPITupleTable *tuptable = NULL;
212 HeapTuple tuple;
213 Datum datum;
214 bool isNull = FALSE;
215
216 rt_pgraster *pgraster = NULL;
217 rt_raster raster = NULL;
218 rt_band band = NULL;
219 int num_bands = 0;
220 uint64_t cK = 0;
221 double cM = 0;
222 double cQ = 0;
223 rt_bandstats stats = NULL;
224 rt_bandstats rtn = NULL;
225
226 Datum values[VALUES_LENGTH];
227 bool nulls[VALUES_LENGTH];
228 Datum result;
229
230 /* tablename is null, return null */
231 if (PG_ARGISNULL(0)) {
232 elog(NOTICE, "Table name must be provided");
233 PG_RETURN_NULL();
234 }
235 tablenametext = PG_GETARG_TEXT_P(0);
236 tablename = text_to_cstring(tablenametext);
237 if (!strlen(tablename)) {
238 elog(NOTICE, "Table name must be provided");
239 PG_RETURN_NULL();
240 }
241
242 /* column name is null, return null */
243 if (PG_ARGISNULL(1)) {
244 elog(NOTICE, "Column name must be provided");
245 PG_RETURN_NULL();
246 }
247 colnametext = PG_GETARG_TEXT_P(1);
248 colname = text_to_cstring(colnametext);
249 if (!strlen(colname)) {
250 elog(NOTICE, "Column name must be provided");
251 PG_RETURN_NULL();
252 }
253
254 /* band index is 1-based */
255 if (!PG_ARGISNULL(2))
256 bandindex = PG_GETARG_INT32(2);
257
258 /* exclude_nodata_value flag */
259 if (!PG_ARGISNULL(3))
260 exclude_nodata_value = PG_GETARG_BOOL(3);
261
262 /* sample % */
263 if (!PG_ARGISNULL(4)) {
264 sample = PG_GETARG_FLOAT8(4);
265 if (sample < 0 || sample > 1) {
266 elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
267 rt_raster_destroy(raster);
268 PG_RETURN_NULL();
269 }
270 else if (FLT_EQ(sample, 0.0))
271 sample = 1;
272 }
273 else
274 sample = 1;
275
276 /* iterate through rasters of coverage */
277 /* connect to database */
278 spi_result = SPI_connect();
279 if (spi_result != SPI_OK_CONNECT) {
280 pfree(sql);
281 elog(ERROR, "RASTER_summaryStatsCoverage: Cannot connect to database using SPI");
282 PG_RETURN_NULL();
283 }
284
285 /* create sql */
286 len = sizeof(char) * (strlen("SELECT \"\" FROM \"\" WHERE \"\" IS NOT NULL") + (strlen(colname) * 2) + strlen(tablename) + 1);
287 sql = (char *) palloc(len);
288 if (NULL == sql) {
289 if (SPI_tuptable) SPI_freetuptable(tuptable);
290 SPI_finish();
291 elog(ERROR, "RASTER_summaryStatsCoverage: Cannot allocate memory for sql");
292 PG_RETURN_NULL();
293 }
294
295 /* get cursor */
296 snprintf(sql, len, "SELECT \"%s\" FROM \"%s\" WHERE \"%s\" IS NOT NULL", colname, tablename, colname);
297 portal = SPI_cursor_open_with_args(
298 "coverage",
299 sql,
300 0, NULL,
301 NULL, NULL,
302 TRUE, 0
303 );
304 pfree(sql);
305
306 /* process resultset */
307 SPI_cursor_fetch(portal, TRUE, 1);
308 while (SPI_processed == 1 && SPI_tuptable != NULL) {
309 tupdesc = SPI_tuptable->tupdesc;
310 tuple = SPI_tuptable->vals[0];
311
312 datum = SPI_getbinval(tuple, tupdesc, 1, &isNull);
313 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
314 SPI_freetuptable(SPI_tuptable);
315 SPI_cursor_close(portal);
316 SPI_finish();
317
318 if (NULL != rtn) pfree(rtn);
319 elog(ERROR, "RASTER_summaryStatsCoverage: Cannot get raster of coverage");
320 PG_RETURN_NULL();
321 }
322 else if (isNull) {
323 SPI_cursor_fetch(portal, TRUE, 1);
324 continue;
325 }
326
327 pgraster = (rt_pgraster *) PG_DETOAST_DATUM(datum);
328
329 raster = rt_raster_deserialize(pgraster, FALSE);
330 if (!raster) {
331 SPI_freetuptable(SPI_tuptable);
332 SPI_cursor_close(portal);
333 SPI_finish();
334
335 if (NULL != rtn) pfree(rtn);
336 elog(ERROR, "RASTER_summaryStatsCoverage: Cannot deserialize raster");
337 PG_RETURN_NULL();
338 }
339
340 /* inspect number of bands */
341 num_bands = rt_raster_get_num_bands(raster);
342 if (bandindex < 1 || bandindex > num_bands) {
343 elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
344
345 rt_raster_destroy(raster);
346
347 SPI_freetuptable(SPI_tuptable);
348 SPI_cursor_close(portal);
349 SPI_finish();
350
351 if (NULL != rtn) pfree(rtn);
352 PG_RETURN_NULL();
353 }
354
355 /* get band */
356 band = rt_raster_get_band(raster, bandindex - 1);
357 if (!band) {
358 elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
359
360 rt_raster_destroy(raster);
361
362 SPI_freetuptable(SPI_tuptable);
363 SPI_cursor_close(portal);
364 SPI_finish();
365
366 if (NULL != rtn) pfree(rtn);
367 PG_RETURN_NULL();
368 }
369
370 /* we don't need the raw values, hence the zero parameter */
371 stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 0, &cK, &cM, &cQ);
372
373 rt_band_destroy(band);
374 rt_raster_destroy(raster);
375
376 if (NULL == stats) {
377 elog(NOTICE, "Cannot compute summary statistics for band at index %d. Returning NULL", bandindex);
378
379 SPI_freetuptable(SPI_tuptable);
380 SPI_cursor_close(portal);
381 SPI_finish();
382
383 if (NULL != rtn) pfree(rtn);
384 PG_RETURN_NULL();
385 }
386
387 /* initialize rtn */
388 if (stats->count > 0) {
389 if (NULL == rtn) {
390 rtn = (rt_bandstats) SPI_palloc(sizeof(struct rt_bandstats_t));
391 if (NULL == rtn) {
392 SPI_freetuptable(SPI_tuptable);
393 SPI_cursor_close(portal);
394 SPI_finish();
395
396 elog(ERROR, "RASTER_summaryStatsCoverage: Cannot allocate memory for summary stats of coverage");
397 PG_RETURN_NULL();
398 }
399
400 rtn->sample = stats->sample;
401 rtn->count = stats->count;
402 rtn->min = stats->min;
403 rtn->max = stats->max;
404 rtn->sum = stats->sum;
405 rtn->mean = stats->mean;
406 rtn->stddev = -1;
407
408 rtn->values = NULL;
409 rtn->sorted = 0;
410 }
411 else {
412 rtn->count += stats->count;
413 rtn->sum += stats->sum;
414
415 if (stats->min < rtn->min)
416 rtn->min = stats->min;
417 if (stats->max > rtn->max)
418 rtn->max = stats->max;
419 }
420 }
421
422 pfree(stats);
423
424 /* next record */
425 SPI_cursor_fetch(portal, TRUE, 1);
426 }
427
428 if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
429 SPI_cursor_close(portal);
430 SPI_finish();
431
432 if (NULL == rtn) {
433 elog(ERROR, "RASTER_summaryStatsCoverage: Cannot compute coverage summary stats");
434 PG_RETURN_NULL();
435 }
436
437 /* coverage mean and deviation */
438 rtn->mean = rtn->sum / rtn->count;
439 /* sample deviation */
440 if (rtn->sample > 0 && rtn->sample < 1)
441 rtn->stddev = sqrt(cQ / (rtn->count - 1));
442 /* standard deviation */
443 else
444 rtn->stddev = sqrt(cQ / rtn->count);
445
446 /* Build a tuple descriptor for our result type */
447 if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
448 ereport(ERROR, (
449 errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
450 errmsg(
451 "function returning record called in context "
452 "that cannot accept type record"
453 )
454 ));
455 }
456
457 BlessTupleDesc(tupdesc);
458
459 memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
460
461 values[0] = Int64GetDatum(rtn->count);
462 if (rtn->count > 0) {
463 values[1] = Float8GetDatum(rtn->sum);
464 values[2] = Float8GetDatum(rtn->mean);
465 values[3] = Float8GetDatum(rtn->stddev);
466 values[4] = Float8GetDatum(rtn->min);
467 values[5] = Float8GetDatum(rtn->max);
468 }
469 else {
470 nulls[1] = TRUE;
471 nulls[2] = TRUE;
472 nulls[3] = TRUE;
473 nulls[4] = TRUE;
474 nulls[5] = TRUE;
475 }
476
477 /* build a tuple */
478 tuple = heap_form_tuple(tupdesc, values, nulls);
479
480 /* make the tuple into a datum */
481 result = HeapTupleGetDatum(tuple);
482
483 /* clean up */
484 pfree(rtn);
485
486 PG_RETURN_DATUM(result);
487}
488
489/* ---------------------------------------------------------------- */
490/* Aggregate ST_SummaryStats */
491/* ---------------------------------------------------------------- */
492
493typedef struct rtpg_summarystats_arg_t *rtpg_summarystats_arg;
494struct rtpg_summarystats_arg_t {
495 rt_bandstats stats;
496
497 /* coefficients for one-pass standard deviation */
498 uint64_t cK;
499 double cM;
500 double cQ;
501
502 int32_t band_index; /* one-based */
503 bool exclude_nodata_value;
504 double sample; /* value between 0 and 1 */
505};
506
507static void
508rtpg_summarystats_arg_destroy(rtpg_summarystats_arg arg) {
509 if (arg->stats != NULL)
510 pfree(arg->stats);
511
512 pfree(arg);
513}
514
515static rtpg_summarystats_arg
516rtpg_summarystats_arg_init() {
517 rtpg_summarystats_arg arg = NULL;
518
519 arg = palloc(sizeof(struct rtpg_summarystats_arg_t));
520 if (arg == NULL) {
521 elog(
522 ERROR,
523 "rtpg_summarystats_arg_init: Cannot allocate memory for function arguments"
524 );
525 return NULL;
526 }
527
528 arg->stats = (rt_bandstats) palloc(sizeof(struct rt_bandstats_t));
529 if (arg->stats == NULL) {
530 rtpg_summarystats_arg_destroy(arg);
531 elog(
532 ERROR,
533 "rtpg_summarystats_arg_init: Cannot allocate memory for stats function argument"
534 );
535 return NULL;
536 }
537
538 arg->stats->sample = 0;
539 arg->stats->count = 0;
540 arg->stats->min = 0;
541 arg->stats->max = 0;
542 arg->stats->sum = 0;
543 arg->stats->mean = 0;
544 arg->stats->stddev = -1;
545 arg->stats->values = NULL;
546 arg->stats->sorted = 0;
547
548 arg->cK = 0;
549 arg->cM = 0;
550 arg->cQ = 0;
551
552 arg->band_index = 1;
553 arg->exclude_nodata_value = TRUE;
554 arg->sample = 1;
555
556 return arg;
557}
558
559PG_FUNCTION_INFO_V1(RASTER_summaryStats_transfn);
560Datum RASTER_summaryStats_transfn(PG_FUNCTION_ARGS)
561{
562 MemoryContext aggcontext;
563 MemoryContext oldcontext;
564 rtpg_summarystats_arg state = NULL;
565 bool skiparg = FALSE;
566
567 int i = 0;
568
569 rt_pgraster *pgraster = NULL;
570 rt_raster raster = NULL;
571 rt_band band = NULL;
572 int num_bands = 0;
573 rt_bandstats stats = NULL;
574
575 POSTGIS_RT_DEBUG(3, "Starting...");
576
577 /* cannot be called directly as this is exclusive aggregate function */
578 if (!AggCheckCallContext(fcinfo, &aggcontext)) {
579 elog(
580 ERROR,
581 "RASTER_summaryStats_transfn: Cannot be called in a non-aggregate context"
582 );
583 PG_RETURN_NULL();
584 }
585
586 /* switch to aggcontext */
587 oldcontext = MemoryContextSwitchTo(aggcontext);
588
589 if (PG_ARGISNULL(0)) {
590 POSTGIS_RT_DEBUG(3, "Creating state variable");
591
592 state = rtpg_summarystats_arg_init();
593 if (state == NULL) {
594 MemoryContextSwitchTo(oldcontext);
595 elog(
596 ERROR,
597 "RASTER_summaryStats_transfn: Cannot allocate memory for state variable"
598 );
599 PG_RETURN_NULL();
600 }
601
602 skiparg = FALSE;
603 }
604 else {
605 POSTGIS_RT_DEBUG(3, "State variable already exists");
606 state = (rtpg_summarystats_arg) PG_GETARG_POINTER(0);
607 skiparg = TRUE;
608 }
609
610 /* raster arg is NOT NULL */
611 if (!PG_ARGISNULL(1)) {
612 /* deserialize raster */
613 pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
614
615 /* Get raster object */
616 raster = rt_raster_deserialize(pgraster, FALSE);
617 if (raster == NULL) {
618
619 rtpg_summarystats_arg_destroy(state);
620 PG_FREE_IF_COPY(pgraster, 1);
621
622 MemoryContextSwitchTo(oldcontext);
623 elog(ERROR, "RASTER_summaryStats_transfn: Cannot deserialize raster");
624 PG_RETURN_NULL();
625 }
626 }
627
628 do {
629 Oid calltype;
630 int nargs = 0;
631
632 if (skiparg)
633 break;
634
635 /* 4 or 5 total possible args */
636 nargs = PG_NARGS();
637 POSTGIS_RT_DEBUGF(4, "nargs = %d", nargs);
638
639 for (i = 2; i < nargs; i++) {
640 if (PG_ARGISNULL(i))
641 continue;
642
643 calltype = get_fn_expr_argtype(fcinfo->flinfo, i);
644
645 /* band index */
646 if (
647 (calltype == INT2OID || calltype == INT4OID) &&
648 i == 2
649 ) {
650 if (calltype == INT2OID)
651 state->band_index = PG_GETARG_INT16(i);
652 else
653 state->band_index = PG_GETARG_INT32(i);
654
655 /* basic check, > 0 */
656 if (state->band_index < 1) {
657
658 rtpg_summarystats_arg_destroy(state);
659 if (raster != NULL) {
660 rt_raster_destroy(raster);
661 PG_FREE_IF_COPY(pgraster, 1);
662 }
663
664 MemoryContextSwitchTo(oldcontext);
665 elog(
666 ERROR,
667 "RASTER_summaryStats_transfn: Invalid band index (must use 1-based). Returning NULL"
668 );
669 PG_RETURN_NULL();
670 }
671 }
672 /* exclude_nodata_value */
673 else if (
674 calltype == BOOLOID && (
675 i == 2 || i == 3
676 )
677 ) {
678 state->exclude_nodata_value = PG_GETARG_BOOL(i);
679 }
680 /* sample rate */
681 else if (
682 (calltype == FLOAT4OID || calltype == FLOAT8OID) &&
683 (i == 3 || i == 4)
684 ) {
685 if (calltype == FLOAT4OID)
686 state->sample = PG_GETARG_FLOAT4(i);
687 else
688 state->sample = PG_GETARG_FLOAT8(i);
689
690 /* basic check, 0 <= sample <= 1 */
691 if (state->sample < 0. || state->sample > 1.) {
692
693 rtpg_summarystats_arg_destroy(state);
694 if (raster != NULL) {
695 rt_raster_destroy(raster);
696 PG_FREE_IF_COPY(pgraster, 1);
697 }
698
699 MemoryContextSwitchTo(oldcontext);
700 elog(
701 ERROR,
702 "Invalid sample percentage (must be between 0 and 1). Returning NULL"
703 );
704
705 PG_RETURN_NULL();
706 }
707 else if (FLT_EQ(state->sample, 0.0))
708 state->sample = 1;
709 }
710 /* unknown arg */
711 else {
712 rtpg_summarystats_arg_destroy(state);
713 if (raster != NULL) {
714 rt_raster_destroy(raster);
715 PG_FREE_IF_COPY(pgraster, 1);
716 }
717
718 MemoryContextSwitchTo(oldcontext);
719 elog(
720 ERROR,
721 "RASTER_summaryStats_transfn: Unknown function parameter at index %d",
722 i
723 );
724 PG_RETURN_NULL();
725 }
726 }
727 }
728 while (0);
729
730 /* null raster, return */
731 if (PG_ARGISNULL(1)) {
732 POSTGIS_RT_DEBUG(4, "NULL raster so processing required");
733 MemoryContextSwitchTo(oldcontext);
734 PG_RETURN_POINTER(state);
735 }
736
737 /* inspect number of bands */
738 num_bands = rt_raster_get_num_bands(raster);
739 if (state->band_index > num_bands) {
740 elog(
741 NOTICE,
742 "Raster does not have band at index %d. Skipping raster",
743 state->band_index
744 );
745
746 rt_raster_destroy(raster);
747 PG_FREE_IF_COPY(pgraster, 1);
748
749 MemoryContextSwitchTo(oldcontext);
750 PG_RETURN_POINTER(state);
751 }
752
753 /* get band */
754 band = rt_raster_get_band(raster, state->band_index - 1);
755 if (!band) {
756 elog(
757 NOTICE, "Cannot find band at index %d. Skipping raster",
758 state->band_index
759 );
760
761 rt_raster_destroy(raster);
762 PG_FREE_IF_COPY(pgraster, 1);
763
764 MemoryContextSwitchTo(oldcontext);
765 PG_RETURN_POINTER(state);
766 }
767
768 /* we don't need the raw values, hence the zero parameter */
769 stats = rt_band_get_summary_stats(
770 band, (int) state->exclude_nodata_value,
771 state->sample, 0,
772 &(state->cK), &(state->cM), &(state->cQ)
773 );
774
775 rt_band_destroy(band);
776 rt_raster_destroy(raster);
777 PG_FREE_IF_COPY(pgraster, 1);
778
779 if (NULL == stats) {
780 elog(
781 NOTICE,
782 "Cannot compute summary statistics for band at index %d. Returning NULL",
783 state->band_index
784 );
785
786 rtpg_summarystats_arg_destroy(state);
787
788 MemoryContextSwitchTo(oldcontext);
789 PG_RETURN_NULL();
790 }
791
792 if (stats->count > 0) {
793 if (state->stats->count < 1) {
794 state->stats->sample = stats->sample;
795 state->stats->count = stats->count;
796 state->stats->min = stats->min;
797 state->stats->max = stats->max;
798 state->stats->sum = stats->sum;
799 state->stats->mean = stats->mean;
800 state->stats->stddev = -1;
801 }
802 else {
803 state->stats->count += stats->count;
804 state->stats->sum += stats->sum;
805
806 if (stats->min < state->stats->min)
807 state->stats->min = stats->min;
808 if (stats->max > state->stats->max)
809 state->stats->max = stats->max;
810 }
811 }
812
813 pfree(stats);
814
815 /* switch back to local context */
816 MemoryContextSwitchTo(oldcontext);
817
818 POSTGIS_RT_DEBUG(3, "Finished");
819
820 PG_RETURN_POINTER(state);
821}
822
823PG_FUNCTION_INFO_V1(RASTER_summaryStats_finalfn);
824Datum RASTER_summaryStats_finalfn(PG_FUNCTION_ARGS)
825{
826 rtpg_summarystats_arg state = NULL;
827
828 TupleDesc tupdesc;
829 HeapTuple tuple;
830 Datum values[VALUES_LENGTH];
831 bool nulls[VALUES_LENGTH];
832 Datum result;
833
834 POSTGIS_RT_DEBUG(3, "Starting...");
835
836 /* cannot be called directly as this is exclusive aggregate function */
837 if (!AggCheckCallContext(fcinfo, NULL)) {
838 elog(ERROR, "RASTER_summaryStats_finalfn: Cannot be called in a non-aggregate context");
839 PG_RETURN_NULL();
840 }
841
842 /* NULL, return null */
843 if (PG_ARGISNULL(0))
844 PG_RETURN_NULL();
845
846 state = (rtpg_summarystats_arg) PG_GETARG_POINTER(0);
847
848 if (NULL == state) {
849 elog(ERROR, "RASTER_summaryStats_finalfn: Cannot compute coverage summary stats");
850 PG_RETURN_NULL();
851 }
852
853 /* coverage mean and deviation */
854 if (state->stats->count > 0) {
855 state->stats->mean = state->stats->sum / state->stats->count;
856 /* sample deviation */
857 if (state->stats->sample > 0 && state->stats->sample < 1)
858 state->stats->stddev = sqrt(state->cQ / (state->stats->count - 1));
859 /* standard deviation */
860 else
861 state->stats->stddev = sqrt(state->cQ / state->stats->count);
862 }
863
864 /* Build a tuple descriptor for our result type */
865 if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
866 rtpg_summarystats_arg_destroy(state);
867 ereport(ERROR, (
868 errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
869 errmsg(
870 "function returning record called in context "
871 "that cannot accept type record"
872 )
873 ));
874 }
875
876 BlessTupleDesc(tupdesc);
877
878 memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
879
880 values[0] = Int64GetDatum(state->stats->count);
881 if (state->stats->count > 0) {
882 values[1] = Float8GetDatum(state->stats->sum);
883 values[2] = Float8GetDatum(state->stats->mean);
884 values[3] = Float8GetDatum(state->stats->stddev);
885 values[4] = Float8GetDatum(state->stats->min);
886 values[5] = Float8GetDatum(state->stats->max);
887 }
888 else {
889 nulls[1] = TRUE;
890 nulls[2] = TRUE;
891 nulls[3] = TRUE;
892 nulls[4] = TRUE;
893 nulls[5] = TRUE;
894 }
895
896 /* build a tuple */
897 tuple = heap_form_tuple(tupdesc, values, nulls);
898
899 /* make the tuple into a datum */
900 result = HeapTupleGetDatum(tuple);
901
902 /* clean up */
903 /* For Windowing functions, it is important to leave */
904 /* the state intact, knowing that the aggcontext will be */
905 /* freed by PgSQL when the statement is complete. */
906 /* https://trac.osgeo.org/postgis/ticket/4770 */
907 // rtpg_summarystats_arg_destroy(state);
908
909 PG_RETURN_DATUM(result);
910}
911
912#undef VALUES_LENGTH
913#define VALUES_LENGTH 4
914
918PG_FUNCTION_INFO_V1(RASTER_histogram);
919Datum 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, (uint32_t)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 Datum values[VALUES_LENGTH];
1161 bool nulls[VALUES_LENGTH];
1162 HeapTuple tuple;
1163 Datum result;
1164
1165 POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
1166
1167 memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
1168
1169 values[0] = Float8GetDatum(hist2[call_cntr].min);
1170 values[1] = Float8GetDatum(hist2[call_cntr].max);
1171 values[2] = Int64GetDatum(hist2[call_cntr].count);
1172 values[3] = Float8GetDatum(hist2[call_cntr].percent);
1173
1174 /* build a tuple */
1175 tuple = heap_form_tuple(tupdesc, values, nulls);
1176
1177 /* make the tuple into a datum */
1178 result = HeapTupleGetDatum(tuple);
1179
1180 SRF_RETURN_NEXT(funcctx, result);
1181 }
1182 /* do when there is no more left */
1183 else {
1184 pfree(hist2);
1185 SRF_RETURN_DONE(funcctx);
1186 }
1187}
1188
1189#undef VALUES_LENGTH
1190#define VALUES_LENGTH 2
1191
1195PG_FUNCTION_INFO_V1(RASTER_quantile);
1196Datum RASTER_quantile(PG_FUNCTION_ARGS)
1197{
1198 FuncCallContext *funcctx;
1199 TupleDesc tupdesc;
1200
1201 int i;
1202 rt_quantile quant;
1203 rt_quantile quant2;
1204 int call_cntr;
1205 int max_calls;
1206
1207 /* first call of function */
1208 if (SRF_IS_FIRSTCALL()) {
1209 MemoryContext oldcontext;
1210
1211 rt_pgraster *pgraster = NULL;
1212 rt_raster raster = NULL;
1213 rt_band band = NULL;
1214 int32_t bandindex = 0;
1215 int num_bands = 0;
1216 bool exclude_nodata_value = TRUE;
1217 double sample = 0;
1218 double *quantiles = NULL;
1219 uint32_t quantiles_count = 0;
1220 double quantile = 0;
1221 rt_bandstats stats = NULL;
1222 uint32_t count;
1223
1224 int j;
1225 int n;
1226
1227 ArrayType *array;
1228 Oid etype;
1229 Datum *e;
1230 bool *nulls;
1231 int16 typlen;
1232 bool typbyval;
1233 char typalign;
1234
1235 /* create a function context for cross-call persistence */
1236 funcctx = SRF_FIRSTCALL_INIT();
1237
1238 /* switch to memory context appropriate for multiple function calls */
1239 oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
1240
1241 /* pgraster is null, return nothing */
1242 if (PG_ARGISNULL(0)) {
1243 MemoryContextSwitchTo(oldcontext);
1244 SRF_RETURN_DONE(funcctx);
1245 }
1246 pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1247
1248 raster = rt_raster_deserialize(pgraster, FALSE);
1249 if (!raster) {
1250 PG_FREE_IF_COPY(pgraster, 0);
1251 MemoryContextSwitchTo(oldcontext);
1252 elog(ERROR, "RASTER_quantile: Cannot deserialize raster");
1253 SRF_RETURN_DONE(funcctx);
1254 }
1255
1256 /* band index is 1-based */
1257 bandindex = PG_GETARG_INT32(1);
1258 num_bands = rt_raster_get_num_bands(raster);
1259 if (bandindex < 1 || bandindex > num_bands) {
1260 elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
1261 rt_raster_destroy(raster);
1262 PG_FREE_IF_COPY(pgraster, 0);
1263 MemoryContextSwitchTo(oldcontext);
1264 SRF_RETURN_DONE(funcctx);
1265 }
1266
1267 /* exclude_nodata_value flag */
1268 if (!PG_ARGISNULL(2))
1269 exclude_nodata_value = PG_GETARG_BOOL(2);
1270
1271 /* sample % */
1272 if (!PG_ARGISNULL(3)) {
1273 sample = PG_GETARG_FLOAT8(3);
1274 if (sample < 0 || sample > 1) {
1275 elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
1276 rt_raster_destroy(raster);
1277 PG_FREE_IF_COPY(pgraster, 0);
1278 MemoryContextSwitchTo(oldcontext);
1279 SRF_RETURN_DONE(funcctx);
1280 }
1281 else if (FLT_EQ(sample, 0.0))
1282 sample = 1;
1283 }
1284 else
1285 sample = 1;
1286
1287 /* quantiles */
1288 if (!PG_ARGISNULL(4)) {
1289 array = PG_GETARG_ARRAYTYPE_P(4);
1290 etype = ARR_ELEMTYPE(array);
1291 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1292
1293 switch (etype) {
1294 case FLOAT4OID:
1295 case FLOAT8OID:
1296 break;
1297 default:
1298 rt_raster_destroy(raster);
1299 PG_FREE_IF_COPY(pgraster, 0);
1300 MemoryContextSwitchTo(oldcontext);
1301 elog(ERROR, "RASTER_quantile: Invalid data type for quantiles");
1302 SRF_RETURN_DONE(funcctx);
1303 break;
1304 }
1305
1306 deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1307 &nulls, &n);
1308
1309 quantiles = palloc(sizeof(double) * n);
1310 for (i = 0, j = 0; i < n; i++) {
1311 if (nulls[i]) continue;
1312
1313 switch (etype) {
1314 case FLOAT4OID:
1315 quantile = (double) DatumGetFloat4(e[i]);
1316 break;
1317 case FLOAT8OID:
1318 quantile = (double) DatumGetFloat8(e[i]);
1319 break;
1320 }
1321
1322 if (quantile < 0 || quantile > 1) {
1323 elog(NOTICE, "Invalid value for quantile (must be between 0 and 1). Returning NULL");
1324 pfree(quantiles);
1325 rt_raster_destroy(raster);
1326 PG_FREE_IF_COPY(pgraster, 0);
1327 MemoryContextSwitchTo(oldcontext);
1328 SRF_RETURN_DONE(funcctx);
1329 }
1330
1331 quantiles[j] = quantile;
1332 POSTGIS_RT_DEBUGF(5, "quantiles[%d] = %f", j, quantiles[j]);
1333 j++;
1334 }
1335 quantiles_count = j;
1336
1337 if (j < 1) {
1338 pfree(quantiles);
1339 quantiles = NULL;
1340 }
1341 }
1342
1343 /* get band */
1344 band = rt_raster_get_band(raster, bandindex - 1);
1345 if (!band) {
1346 elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
1347 rt_raster_destroy(raster);
1348 PG_FREE_IF_COPY(pgraster, 0);
1349 MemoryContextSwitchTo(oldcontext);
1350 SRF_RETURN_DONE(funcctx);
1351 }
1352
1353 /* get stats */
1354 stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 1, NULL, NULL, NULL);
1355 rt_band_destroy(band);
1356 rt_raster_destroy(raster);
1357 PG_FREE_IF_COPY(pgraster, 0);
1358 if (NULL == stats || NULL == stats->values) {
1359 elog(NOTICE, "Cannot retrieve summary statistics for band at index %d", bandindex);
1360 MemoryContextSwitchTo(oldcontext);
1361 SRF_RETURN_DONE(funcctx);
1362 }
1363 else if (stats->count < 1) {
1364 elog(NOTICE, "Cannot compute quantiles for band at index %d as the band has no values", bandindex);
1365 MemoryContextSwitchTo(oldcontext);
1366 SRF_RETURN_DONE(funcctx);
1367 }
1368
1369 /* get quantiles */
1370 quant = rt_band_get_quantiles(stats, quantiles, quantiles_count, &count);
1371 if (quantiles_count) pfree(quantiles);
1372 pfree(stats);
1373 if (NULL == quant || !count) {
1374 elog(NOTICE, "Cannot compute quantiles for band at index %d", bandindex);
1375 MemoryContextSwitchTo(oldcontext);
1376 SRF_RETURN_DONE(funcctx);
1377 }
1378
1379 POSTGIS_RT_DEBUGF(3, "%d quantiles returned", count);
1380
1381 /* Store needed information */
1382 funcctx->user_fctx = quant;
1383
1384 /* total number of tuples to be returned */
1385 funcctx->max_calls = count;
1386
1387 /* Build a tuple descriptor for our result type */
1388 if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
1389 ereport(ERROR, (
1390 errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1391 errmsg(
1392 "function returning record called in context "
1393 "that cannot accept type record"
1394 )
1395 ));
1396 }
1397
1398 BlessTupleDesc(tupdesc);
1399 funcctx->tuple_desc = tupdesc;
1400
1401 MemoryContextSwitchTo(oldcontext);
1402 }
1403
1404 /* stuff done on every call of the function */
1405 funcctx = SRF_PERCALL_SETUP();
1406
1407 call_cntr = funcctx->call_cntr;
1408 max_calls = funcctx->max_calls;
1409 tupdesc = funcctx->tuple_desc;
1410 quant2 = funcctx->user_fctx;
1411
1412 /* do when there is more left to send */
1413 if (call_cntr < max_calls) {
1414 Datum values[VALUES_LENGTH];
1415 bool nulls[VALUES_LENGTH];
1416 HeapTuple tuple;
1417 Datum result;
1418
1419 POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
1420
1421 memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
1422
1423 values[0] = Float8GetDatum(quant2[call_cntr].quantile);
1424 values[1] = Float8GetDatum(quant2[call_cntr].value);
1425
1426 /* build a tuple */
1427 tuple = heap_form_tuple(tupdesc, values, nulls);
1428
1429 /* make the tuple into a datum */
1430 result = HeapTupleGetDatum(tuple);
1431
1432 SRF_RETURN_NEXT(funcctx, result);
1433 }
1434 /* do when there is no more left */
1435 else {
1436 pfree(quant2);
1437 SRF_RETURN_DONE(funcctx);
1438 }
1439}
1440
1441#undef VALUES_LENGTH
1442#define VALUES_LENGTH 3
1443
1444/* get counts of values */
1445PG_FUNCTION_INFO_V1(RASTER_valueCount);
1446Datum RASTER_valueCount(PG_FUNCTION_ARGS) {
1447 FuncCallContext *funcctx;
1448 TupleDesc tupdesc;
1449
1450 int i;
1451 rt_valuecount vcnts;
1452 rt_valuecount vcnts2;
1453 int call_cntr;
1454 int max_calls;
1455
1456 /* first call of function */
1457 if (SRF_IS_FIRSTCALL()) {
1458 MemoryContext oldcontext;
1459
1460 rt_pgraster *pgraster = NULL;
1461 rt_raster raster = NULL;
1462 rt_band band = NULL;
1463 int32_t bandindex = 0;
1464 int num_bands = 0;
1465 bool exclude_nodata_value = TRUE;
1466 double *search_values = NULL;
1467 uint32_t search_values_count = 0;
1468 double roundto = 0;
1469 uint32_t count;
1470
1471 int j;
1472 int n;
1473
1474 ArrayType *array;
1475 Oid etype;
1476 Datum *e;
1477 bool *nulls;
1478 int16 typlen;
1479 bool typbyval;
1480 char typalign;
1481
1482 /* create a function context for cross-call persistence */
1483 funcctx = SRF_FIRSTCALL_INIT();
1484
1485 /* switch to memory context appropriate for multiple function calls */
1486 oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
1487
1488 /* pgraster is null, return nothing */
1489 if (PG_ARGISNULL(0)) {
1490 MemoryContextSwitchTo(oldcontext);
1491 SRF_RETURN_DONE(funcctx);
1492 }
1493 pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1494
1495 raster = rt_raster_deserialize(pgraster, FALSE);
1496 if (!raster) {
1497 PG_FREE_IF_COPY(pgraster, 0);
1498 MemoryContextSwitchTo(oldcontext);
1499 elog(ERROR, "RASTER_valueCount: Cannot deserialize raster");
1500 SRF_RETURN_DONE(funcctx);
1501 }
1502
1503 /* band index is 1-based */
1504 bandindex = PG_GETARG_INT32(1);
1505 num_bands = rt_raster_get_num_bands(raster);
1506 if (bandindex < 1 || bandindex > num_bands) {
1507 elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
1508 rt_raster_destroy(raster);
1509 PG_FREE_IF_COPY(pgraster, 0);
1510 MemoryContextSwitchTo(oldcontext);
1511 SRF_RETURN_DONE(funcctx);
1512 }
1513
1514 /* exclude_nodata_value flag */
1515 if (!PG_ARGISNULL(2))
1516 exclude_nodata_value = PG_GETARG_BOOL(2);
1517
1518 /* search values */
1519 if (!PG_ARGISNULL(3)) {
1520 array = PG_GETARG_ARRAYTYPE_P(3);
1521 etype = ARR_ELEMTYPE(array);
1522 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1523
1524 switch (etype) {
1525 case FLOAT4OID:
1526 case FLOAT8OID:
1527 break;
1528 default:
1529 rt_raster_destroy(raster);
1530 PG_FREE_IF_COPY(pgraster, 0);
1531 MemoryContextSwitchTo(oldcontext);
1532 elog(ERROR, "RASTER_valueCount: Invalid data type for values");
1533 SRF_RETURN_DONE(funcctx);
1534 break;
1535 }
1536
1537 deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1538 &nulls, &n);
1539
1540 search_values = palloc(sizeof(double) * n);
1541 for (i = 0, j = 0; i < n; i++) {
1542 if (nulls[i]) continue;
1543
1544 switch (etype) {
1545 case FLOAT4OID:
1546 search_values[j] = (double) DatumGetFloat4(e[i]);
1547 break;
1548 case FLOAT8OID:
1549 search_values[j] = (double) DatumGetFloat8(e[i]);
1550 break;
1551 }
1552
1553 POSTGIS_RT_DEBUGF(5, "search_values[%d] = %f", j, search_values[j]);
1554 j++;
1555 }
1556 search_values_count = j;
1557
1558 if (j < 1) {
1559 pfree(search_values);
1560 search_values = NULL;
1561 }
1562 }
1563
1564 /* roundto */
1565 if (!PG_ARGISNULL(4)) {
1566 roundto = PG_GETARG_FLOAT8(4);
1567 if (roundto < 0.) roundto = 0;
1568 }
1569
1570 /* get band */
1571 band = rt_raster_get_band(raster, bandindex - 1);
1572 if (!band) {
1573 elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
1574 rt_raster_destroy(raster);
1575 PG_FREE_IF_COPY(pgraster, 0);
1576 MemoryContextSwitchTo(oldcontext);
1577 SRF_RETURN_DONE(funcctx);
1578 }
1579
1580 /* get counts of values */
1581 vcnts = rt_band_get_value_count(band, (int) exclude_nodata_value, search_values, search_values_count, roundto, NULL, &count);
1582 rt_band_destroy(band);
1583 rt_raster_destroy(raster);
1584 PG_FREE_IF_COPY(pgraster, 0);
1585 if (NULL == vcnts || !count) {
1586 elog(NOTICE, "Cannot count the values for band at index %d", bandindex);
1587 MemoryContextSwitchTo(oldcontext);
1588 SRF_RETURN_DONE(funcctx);
1589 }
1590
1591 POSTGIS_RT_DEBUGF(3, "%d value counts returned", count);
1592
1593 /* Store needed information */
1594 funcctx->user_fctx = vcnts;
1595
1596 /* total number of tuples to be returned */
1597 funcctx->max_calls = count;
1598
1599 /* Build a tuple descriptor for our result type */
1600 if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
1601 ereport(ERROR, (
1602 errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1603 errmsg(
1604 "function returning record called in context "
1605 "that cannot accept type record"
1606 )
1607 ));
1608 }
1609
1610 BlessTupleDesc(tupdesc);
1611 funcctx->tuple_desc = tupdesc;
1612
1613 MemoryContextSwitchTo(oldcontext);
1614 }
1615
1616 /* stuff done on every call of the function */
1617 funcctx = SRF_PERCALL_SETUP();
1618
1619 call_cntr = funcctx->call_cntr;
1620 max_calls = funcctx->max_calls;
1621 tupdesc = funcctx->tuple_desc;
1622 vcnts2 = funcctx->user_fctx;
1623
1624 /* do when there is more left to send */
1625 if (call_cntr < max_calls) {
1626 Datum values[VALUES_LENGTH];
1627 bool nulls[VALUES_LENGTH];
1628 HeapTuple tuple;
1629 Datum result;
1630
1631 POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
1632
1633 memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
1634
1635 values[0] = Float8GetDatum(vcnts2[call_cntr].value);
1636 values[1] = UInt32GetDatum(vcnts2[call_cntr].count);
1637 values[2] = Float8GetDatum(vcnts2[call_cntr].percent);
1638
1639 /* build a tuple */
1640 tuple = heap_form_tuple(tupdesc, values, nulls);
1641
1642 /* make the tuple into a datum */
1643 result = HeapTupleGetDatum(tuple);
1644
1645 SRF_RETURN_NEXT(funcctx, result);
1646 }
1647 /* do when there is no more left */
1648 else {
1649 pfree(vcnts2);
1650 SRF_RETURN_DONE(funcctx);
1651 }
1652}
1653
1654/* get counts of values for a coverage */
1655PG_FUNCTION_INFO_V1(RASTER_valueCountCoverage);
1656Datum RASTER_valueCountCoverage(PG_FUNCTION_ARGS) {
1657 FuncCallContext *funcctx;
1658 TupleDesc tupdesc;
1659
1660 uint32_t i;
1661 uint64_t covcount = 0;
1662 uint64_t covtotal = 0;
1663 rt_valuecount covvcnts = NULL;
1664 rt_valuecount covvcnts2;
1665 int call_cntr;
1666 int max_calls;
1667
1668 POSTGIS_RT_DEBUG(3, "RASTER_valueCountCoverage: Starting");
1669
1670 /* first call of function */
1671 if (SRF_IS_FIRSTCALL()) {
1672 MemoryContext oldcontext;
1673
1674 text *tablenametext = NULL;
1675 char *tablename = NULL;
1676 text *colnametext = NULL;
1677 char *colname = NULL;
1678 int32_t bandindex = 1;
1679 bool exclude_nodata_value = TRUE;
1680 double *search_values = NULL;
1681 uint32_t search_values_count = 0;
1682 double roundto = 0;
1683
1684 int len = 0;
1685 char *sql = NULL;
1686 int spi_result;
1687 Portal portal;
1688 HeapTuple tuple;
1689 Datum datum;
1690 bool isNull = FALSE;
1691 rt_pgraster *pgraster = NULL;
1692 rt_raster raster = NULL;
1693 rt_band band = NULL;
1694 int num_bands = 0;
1695 uint32_t count;
1696 uint32_t total;
1697 rt_valuecount vcnts = NULL;
1698 int exists = 0;
1699
1700 uint32_t j;
1701 int n;
1702
1703 ArrayType *array;
1704 Oid etype;
1705 Datum *e;
1706 bool *nulls;
1707 int16 typlen;
1708 bool typbyval;
1709 char typalign;
1710
1711 /* create a function context for cross-call persistence */
1712 funcctx = SRF_FIRSTCALL_INIT();
1713
1714 /* switch to memory context appropriate for multiple function calls */
1715 oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
1716
1717 /* tablename is null, return null */
1718 if (PG_ARGISNULL(0)) {
1719 elog(NOTICE, "Table name must be provided");
1720 MemoryContextSwitchTo(oldcontext);
1721 SRF_RETURN_DONE(funcctx);
1722 }
1723 tablenametext = PG_GETARG_TEXT_P(0);
1724 tablename = text_to_cstring(tablenametext);
1725 if (!strlen(tablename)) {
1726 elog(NOTICE, "Table name must be provided");
1727 MemoryContextSwitchTo(oldcontext);
1728 SRF_RETURN_DONE(funcctx);
1729 }
1730 POSTGIS_RT_DEBUGF(3, "tablename = %s", tablename);
1731
1732 /* column name is null, return null */
1733 if (PG_ARGISNULL(1)) {
1734 elog(NOTICE, "Column name must be provided");
1735 MemoryContextSwitchTo(oldcontext);
1736 SRF_RETURN_DONE(funcctx);
1737 }
1738 colnametext = PG_GETARG_TEXT_P(1);
1739 colname = text_to_cstring(colnametext);
1740 if (!strlen(colname)) {
1741 elog(NOTICE, "Column name must be provided");
1742 MemoryContextSwitchTo(oldcontext);
1743 SRF_RETURN_DONE(funcctx);
1744 }
1745 POSTGIS_RT_DEBUGF(3, "colname = %s", colname);
1746
1747 /* band index is 1-based */
1748 if (!PG_ARGISNULL(2))
1749 bandindex = PG_GETARG_INT32(2);
1750
1751 /* exclude_nodata_value flag */
1752 if (!PG_ARGISNULL(3))
1753 exclude_nodata_value = PG_GETARG_BOOL(3);
1754
1755 /* search values */
1756 if (!PG_ARGISNULL(4)) {
1757 array = PG_GETARG_ARRAYTYPE_P(4);
1758 etype = ARR_ELEMTYPE(array);
1759 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1760
1761 switch (etype) {
1762 case FLOAT4OID:
1763 case FLOAT8OID:
1764 break;
1765 default:
1766 MemoryContextSwitchTo(oldcontext);
1767 elog(ERROR, "RASTER_valueCountCoverage: Invalid data type for values");
1768 SRF_RETURN_DONE(funcctx);
1769 break;
1770 }
1771
1772 deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1773 &nulls, &n);
1774
1775 search_values = palloc(sizeof(double) * n);
1776 for (i = 0, j = 0; i < (uint32_t) n; i++) {
1777 if (nulls[i]) continue;
1778
1779 switch (etype) {
1780 case FLOAT4OID:
1781 search_values[j] = (double) DatumGetFloat4(e[i]);
1782 break;
1783 case FLOAT8OID:
1784 search_values[j] = (double) DatumGetFloat8(e[i]);
1785 break;
1786 }
1787
1788 POSTGIS_RT_DEBUGF(5, "search_values[%d] = %f", j, search_values[j]);
1789 j++;
1790 }
1791 search_values_count = j;
1792
1793 if (j < 1) {
1794 pfree(search_values);
1795 search_values = NULL;
1796 }
1797 }
1798
1799 /* roundto */
1800 if (!PG_ARGISNULL(5)) {
1801 roundto = PG_GETARG_FLOAT8(5);
1802 if (roundto < 0.) roundto = 0;
1803 }
1804
1805 /* iterate through rasters of coverage */
1806 /* connect to database */
1807 spi_result = SPI_connect();
1808 if (spi_result != SPI_OK_CONNECT) {
1809
1810 if (search_values_count) pfree(search_values);
1811
1812 MemoryContextSwitchTo(oldcontext);
1813 elog(ERROR, "RASTER_valueCountCoverage: Cannot connect to database using SPI");
1814 SRF_RETURN_DONE(funcctx);
1815 }
1816
1817 /* create sql */
1818 len = sizeof(char) * (strlen("SELECT \"\" FROM \"\" WHERE \"\" IS NOT NULL") + (strlen(colname) * 2) + strlen(tablename) + 1);
1819 sql = (char *) palloc(len);
1820 if (NULL == sql) {
1821
1822 if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
1823 SPI_finish();
1824
1825 if (search_values_count) pfree(search_values);
1826
1827 MemoryContextSwitchTo(oldcontext);
1828 elog(ERROR, "RASTER_valueCountCoverage: Cannot allocate memory for sql");
1829 SRF_RETURN_DONE(funcctx);
1830 }
1831
1832 /* get cursor */
1833 snprintf(sql, len, "SELECT \"%s\" FROM \"%s\" WHERE \"%s\" IS NOT NULL", colname, tablename, colname);
1834 POSTGIS_RT_DEBUGF(3, "RASTER_valueCountCoverage: %s", sql);
1835 portal = SPI_cursor_open_with_args(
1836 "coverage",
1837 sql,
1838 0, NULL,
1839 NULL, NULL,
1840 TRUE, 0
1841 );
1842 pfree(sql);
1843
1844 /* process resultset */
1845 SPI_cursor_fetch(portal, TRUE, 1);
1846 while (SPI_processed == 1 && SPI_tuptable != NULL) {
1847 tupdesc = SPI_tuptable->tupdesc;
1848 tuple = SPI_tuptable->vals[0];
1849
1850 datum = SPI_getbinval(tuple, tupdesc, 1, &isNull);
1851 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1852
1853 SPI_freetuptable(SPI_tuptable);
1854 SPI_cursor_close(portal);
1855 SPI_finish();
1856
1857 if (NULL != covvcnts) pfree(covvcnts);
1858 if (search_values_count) pfree(search_values);
1859
1860 MemoryContextSwitchTo(oldcontext);
1861 elog(ERROR, "RASTER_valueCountCoverage: Cannot get raster of coverage");
1862 SRF_RETURN_DONE(funcctx);
1863 }
1864 else if (isNull) {
1865 SPI_cursor_fetch(portal, TRUE, 1);
1866 continue;
1867 }
1868
1869 pgraster = (rt_pgraster *) PG_DETOAST_DATUM(datum);
1870
1871 raster = rt_raster_deserialize(pgraster, FALSE);
1872 if (!raster) {
1873
1874 SPI_freetuptable(SPI_tuptable);
1875 SPI_cursor_close(portal);
1876 SPI_finish();
1877
1878 if (NULL != covvcnts) pfree(covvcnts);
1879 if (search_values_count) pfree(search_values);
1880
1881 MemoryContextSwitchTo(oldcontext);
1882 elog(ERROR, "RASTER_valueCountCoverage: Cannot deserialize raster");
1883 SRF_RETURN_DONE(funcctx);
1884 }
1885
1886 /* inspect number of bands*/
1887 num_bands = rt_raster_get_num_bands(raster);
1888 if (bandindex < 1 || bandindex > num_bands) {
1889 elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
1890
1891 rt_raster_destroy(raster);
1892
1893 SPI_freetuptable(SPI_tuptable);
1894 SPI_cursor_close(portal);
1895 SPI_finish();
1896
1897 if (NULL != covvcnts) pfree(covvcnts);
1898 if (search_values_count) pfree(search_values);
1899
1900 MemoryContextSwitchTo(oldcontext);
1901 SRF_RETURN_DONE(funcctx);
1902 }
1903
1904 /* get band */
1905 band = rt_raster_get_band(raster, bandindex - 1);
1906 if (!band) {
1907 elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
1908
1909 rt_raster_destroy(raster);
1910
1911 SPI_freetuptable(SPI_tuptable);
1912 SPI_cursor_close(portal);
1913 SPI_finish();
1914
1915 if (NULL != covvcnts) pfree(covvcnts);
1916 if (search_values_count) pfree(search_values);
1917
1918 MemoryContextSwitchTo(oldcontext);
1919 SRF_RETURN_DONE(funcctx);
1920 }
1921
1922 /* get counts of values */
1923 vcnts = rt_band_get_value_count(band, (int) exclude_nodata_value, search_values, search_values_count, roundto, &total, &count);
1924 rt_band_destroy(band);
1925 rt_raster_destroy(raster);
1926 if (NULL == vcnts || !count) {
1927 elog(NOTICE, "Cannot count the values for band at index %d", bandindex);
1928
1929 SPI_freetuptable(SPI_tuptable);
1930 SPI_cursor_close(portal);
1931 SPI_finish();
1932
1933 if (NULL != covvcnts) free(covvcnts);
1934 if (search_values_count) pfree(search_values);
1935
1936 MemoryContextSwitchTo(oldcontext);
1937 SRF_RETURN_DONE(funcctx);
1938 }
1939
1940 POSTGIS_RT_DEBUGF(3, "%d value counts returned", count);
1941
1942 if (NULL == covvcnts) {
1943 covvcnts = (rt_valuecount) SPI_palloc(sizeof(struct rt_valuecount_t) * count);
1944 if (NULL == covvcnts) {
1945
1946 SPI_freetuptable(SPI_tuptable);
1947 SPI_cursor_close(portal);
1948 SPI_finish();
1949
1950 if (search_values_count) pfree(search_values);
1951
1952 MemoryContextSwitchTo(oldcontext);
1953 elog(ERROR, "RASTER_valueCountCoverage: Cannot allocate memory for value counts of coverage");
1954 SRF_RETURN_DONE(funcctx);
1955 }
1956
1957 for (i = 0; i < count; i++) {
1958 covvcnts[i].value = vcnts[i].value;
1959 covvcnts[i].count = vcnts[i].count;
1960 covvcnts[i].percent = -1;
1961 }
1962
1963 covcount = count;
1964 }
1965 else {
1966 for (i = 0; i < count; i++) {
1967 exists = 0;
1968
1969 for (j = 0; j < covcount; j++) {
1970 if (FLT_EQ(vcnts[i].value, covvcnts[j].value)) {
1971 exists = 1;
1972 break;
1973 }
1974 }
1975
1976 if (exists) {
1977 covvcnts[j].count += vcnts[i].count;
1978 }
1979 else {
1980 covcount++;
1981 covvcnts = SPI_repalloc(covvcnts, sizeof(struct rt_valuecount_t) * covcount);
1982 if (!covvcnts) {
1983 SPI_freetuptable(SPI_tuptable);
1984 SPI_cursor_close(portal);
1985 SPI_finish();
1986
1987 if (search_values_count) pfree(search_values);
1988
1989 MemoryContextSwitchTo(oldcontext);
1990 elog(ERROR, "RASTER_valueCountCoverage: Cannot change allocated memory for value counts of coverage");
1991 SRF_RETURN_DONE(funcctx);
1992 }
1993
1994 covvcnts[covcount - 1].value = vcnts[i].value;
1995 covvcnts[covcount - 1].count = vcnts[i].count;
1996 covvcnts[covcount - 1].percent = -1;
1997 }
1998 }
1999 }
2000
2001 covtotal += total;
2002
2003 pfree(vcnts);
2004
2005 /* next record */
2006 SPI_cursor_fetch(portal, TRUE, 1);
2007 }
2008
2009 if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
2010 SPI_cursor_close(portal);
2011 SPI_finish();
2012
2013 if (search_values_count) pfree(search_values);
2014
2015 /* compute percentages */
2016 for (i = 0; i < covcount; i++) {
2017 covvcnts[i].percent = (double) covvcnts[i].count / covtotal;
2018 }
2019
2020 /* Store needed information */
2021 funcctx->user_fctx = covvcnts;
2022
2023 /* total number of tuples to be returned */
2024 funcctx->max_calls = covcount;
2025
2026 /* Build a tuple descriptor for our result type */
2027 if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
2028 ereport(ERROR, (
2029 errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
2030 errmsg(
2031 "function returning record called in context "
2032 "that cannot accept type record"
2033 )
2034 ));
2035 }
2036
2037 BlessTupleDesc(tupdesc);
2038 funcctx->tuple_desc = tupdesc;
2039
2040 MemoryContextSwitchTo(oldcontext);
2041 }
2042
2043 /* stuff done on every call of the function */
2044 funcctx = SRF_PERCALL_SETUP();
2045
2046 call_cntr = funcctx->call_cntr;
2047 max_calls = funcctx->max_calls;
2048 tupdesc = funcctx->tuple_desc;
2049 covvcnts2 = funcctx->user_fctx;
2050
2051 /* do when there is more left to send */
2052 if (call_cntr < max_calls) {
2053 Datum values[VALUES_LENGTH];
2054 bool nulls[VALUES_LENGTH];
2055 HeapTuple tuple;
2056 Datum result;
2057
2058 POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
2059
2060 memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
2061
2062 values[0] = Float8GetDatum(covvcnts2[call_cntr].value);
2063 values[1] = UInt32GetDatum(covvcnts2[call_cntr].count);
2064 values[2] = Float8GetDatum(covvcnts2[call_cntr].percent);
2065
2066 /* build a tuple */
2067 tuple = heap_form_tuple(tupdesc, values, nulls);
2068
2069 /* make the tuple into a datum */
2070 result = HeapTupleGetDatum(tuple);
2071
2072 SRF_RETURN_NEXT(funcctx, result);
2073 }
2074 /* do when there is no more left */
2075 else {
2076 if (covvcnts2) pfree(covvcnts2);
2077 SRF_RETURN_DONE(funcctx);
2078 }
2079}
2080
result
char result[OUT_DOUBLE_BUFFER_SIZE]
Definition cu_print.c:267
TRUE
#define TRUE
Definition dbfopen.c:73
FALSE
#define FALSE
Definition dbfopen.c:72
rt_raster_destroy
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition rt_raster.c:86
FLT_EQ
#define FLT_EQ(x, y)
Definition librtcore.h:2436
rt_bandstats
struct rt_bandstats_t * rt_bandstats
Definition librtcore.h:151
rt_band_destroy
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition rt_band.c:499
rt_raster_get_num_bands
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition rt_raster.c:376
rt_valuecount
struct rt_valuecount_t * rt_valuecount
Definition librtcore.h:154
rt_band_get_histogram
rt_histogram rt_band_get_histogram(rt_bandstats stats, uint32_t bin_count, double *bin_widths, uint32_t bin_widths_count, int right, double min, double max, uint32_t *rtn_count)
Count the distribution of data.
Definition rt_statistics.c:413
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:1609
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:673
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:111
rt_raster_deserialize
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition rt_serialize.c:732
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:385
free
void free(void *)
RASTER_summaryStats
Datum RASTER_summaryStats(PG_FUNCTION_ARGS)
Definition rtpg_statistics.c:71
RASTER_summaryStatsCoverage
Datum RASTER_summaryStatsCoverage(PG_FUNCTION_ARGS)
Definition rtpg_statistics.c:196
rtpg_summarystats_arg
struct rtpg_summarystats_arg_t * rtpg_summarystats_arg
Definition rtpg_statistics.c:493
RASTER_histogram
Datum RASTER_histogram(PG_FUNCTION_ARGS)
Definition rtpg_statistics.c:919
rtpg_summarystats_arg_destroy
static void rtpg_summarystats_arg_destroy(rtpg_summarystats_arg arg)
Definition rtpg_statistics.c:508
RASTER_valueCountCoverage
Datum RASTER_valueCountCoverage(PG_FUNCTION_ARGS)
Definition rtpg_statistics.c:1656
RASTER_valueCount
Datum RASTER_valueCount(PG_FUNCTION_ARGS)
Definition rtpg_statistics.c:1446
rtpg_summarystats_arg_init
static rtpg_summarystats_arg rtpg_summarystats_arg_init()
Definition rtpg_statistics.c:516
VALUES_LENGTH
#define VALUES_LENGTH
Definition rtpg_statistics.c:65
RASTER_summaryStats_transfn
Datum RASTER_summaryStats_transfn(PG_FUNCTION_ARGS)
Definition rtpg_statistics.c:560
RASTER_quantile
Datum RASTER_quantile(PG_FUNCTION_ARGS)
Definition rtpg_statistics.c:1196
RASTER_summaryStats_finalfn
Datum RASTER_summaryStats_finalfn(PG_FUNCTION_ARGS)
Definition rtpg_statistics.c:824
PG_FUNCTION_INFO_V1
PG_FUNCTION_INFO_V1(RASTER_summaryStats)
Get summary stats of a band.
POSTGIS_RT_DEBUG
#define POSTGIS_RT_DEBUG(level, msg)
Definition rtpostgis.h:65
POSTGIS_RT_DEBUGF
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition rtpostgis.h:69
rt_band_t
Definition librtcore.h:2515
rt_bandstats_t::mean
double mean
Definition librtcore.h:2567
rt_bandstats_t::sorted
int sorted
Definition librtcore.h:2571
rt_bandstats_t::sample
double sample
Definition librtcore.h:2561
rt_bandstats_t::sum
double sum
Definition librtcore.h:2566
rt_bandstats_t::max
double max
Definition librtcore.h:2565
rt_bandstats_t::count
uint32_t count
Definition librtcore.h:2562
rt_bandstats_t::min
double min
Definition librtcore.h:2564
rt_bandstats_t::values
double * values
Definition librtcore.h:2570
rt_bandstats_t::stddev
double stddev
Definition librtcore.h:2568
rt_bandstats_t
Definition librtcore.h:2560
rt_histogram_t
Definition librtcore.h:2575
rt_quantile_t
Definition librtcore.h:2587
rt_raster_serialized_t
Struct definitions.
Definition librtcore.h:2452
rt_raster_t
Definition librtcore.h:2486
rt_valuecount_t::count
uint32_t count
Definition librtcore.h:2627
rt_valuecount_t::percent
double percent
Definition librtcore.h:2628
rt_valuecount_t::value
double value
Definition librtcore.h:2626
rt_valuecount_t
Definition librtcore.h:2625
rtpg_summarystats_arg_t::cQ
double cQ
Definition rtpg_statistics.c:500
rtpg_summarystats_arg_t::cM
double cM
Definition rtpg_statistics.c:499
rtpg_summarystats_arg_t::band_index
int32_t band_index
Definition rtpg_statistics.c:502
rtpg_summarystats_arg_t::cK
uint64_t cK
Definition rtpg_statistics.c:498
rtpg_summarystats_arg_t::stats
rt_bandstats stats
Definition rtpg_statistics.c:495
rtpg_summarystats_arg_t::sample
double sample
Definition rtpg_statistics.c:504
rtpg_summarystats_arg_t::exclude_nodata_value
bool exclude_nodata_value
Definition rtpg_statistics.c:503
rtpg_summarystats_arg_t
Definition rtpg_statistics.c:494