PostGIS 3.7.0dev-r@@SVN_REVISION@@
Loading...
Searching...
No Matches

◆ RASTER_summaryStatsCoverage()

Datum RASTER_summaryStatsCoverage ( PG_FUNCTION_ARGS  )

Definition at line 196 of file rtpg_statistics.c.

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
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}
char result[OUT_DOUBLE_BUFFER_SIZE]
Definition cu_print.c:267
#define TRUE
Definition dbfopen.c:73
#define FALSE
Definition dbfopen.c:72
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition rt_raster.c:86
#define FLT_EQ(x, y)
Definition librtcore.h:2436
struct rt_bandstats_t * rt_bandstats
Definition librtcore.h:151
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition rt_band.c:499
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition rt_raster.c:376
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.
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition rt_raster.c:385
raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition rtrowdump.py:125
#define VALUES_LENGTH
uint32_t count
Definition librtcore.h:2562
double * values
Definition librtcore.h:2570
Struct definitions.
Definition librtcore.h:2452

References rt_bandstats_t::count, FALSE, FLT_EQ, rt_bandstats_t::max, rt_bandstats_t::mean, rt_bandstats_t::min, result, rt_band_destroy(), rt_band_get_summary_stats(), rt_raster_deserialize(), rt_raster_destroy(), rt_raster_get_band(), rt_raster_get_num_bands(), rt_bandstats_t::sample, rt_bandstats_t::sorted, rt_bandstats_t::stddev, rt_bandstats_t::sum, TRUE, rt_bandstats_t::values, and VALUES_LENGTH.

Here is the call graph for this function: