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