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

◆ RASTER_summaryStatsCoverage()

Datum RASTER_summaryStatsCoverage ( PG_FUNCTION_ARGS  )

Definition at line 198 of file rtpg_statistics.c.

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
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}
#define TRUE
Definition dbfopen.c:169
#define FALSE
Definition dbfopen.c:168
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition rt_raster.c:82
#define FLT_EQ(x, y)
Definition librtcore.h:2235
struct rt_bandstats_t * rt_bandstats
Definition librtcore.h:150
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition rt_band.c:340
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition rt_raster.c:372
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:381
raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition rtrowdump.py:121
char * text_to_cstring(const text *textptr)
#define VALUES_LENGTH
uint32_t count
Definition librtcore.h:2361
double * values
Definition librtcore.h:2369
Struct definitions.
Definition librtcore.h:2251

References rt_bandstats_t::count, FALSE, FLT_EQ, rt_bandstats_t::max, rt_bandstats_t::mean, rt_bandstats_t::min, 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, text_to_cstring(), TRUE, rt_bandstats_t::values, and VALUES_LENGTH.

Here is the call graph for this function: