PostGIS  3.4.0dev-r@@SVN_REVISION@@

◆ 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");
268  PG_RETURN_NULL();
269  }
270  else if (FLT_EQ(sample, 0.0))
271  sample = 1;
272  }
273  else
274  sample = 1;
275 
276  /* iterate through rasters of coverage */
277  /* connect to database */
278  spi_result = SPI_connect();
279  if (spi_result != SPI_OK_CONNECT) {
280  pfree(sql);
281  elog(ERROR, "RASTER_summaryStatsCoverage: Cannot connect to database using SPI");
282  PG_RETURN_NULL();
283  }
284 
285  /* create sql */
286  len = sizeof(char) * (strlen("SELECT \"\" FROM \"\" WHERE \"\" IS NOT NULL") + (strlen(colname) * 2) + strlen(tablename) + 1);
287  sql = (char *) palloc(len);
288  if (NULL == sql) {
289  if (SPI_tuptable) SPI_freetuptable(tuptable);
290  SPI_finish();
291  elog(ERROR, "RASTER_summaryStatsCoverage: Cannot allocate memory for sql");
292  PG_RETURN_NULL();
293  }
294 
295  /* get cursor */
296  snprintf(sql, len, "SELECT \"%s\" FROM \"%s\" WHERE \"%s\" IS NOT NULL", colname, tablename, colname);
297  portal = SPI_cursor_open_with_args(
298  "coverage",
299  sql,
300  0, NULL,
301  NULL, NULL,
302  TRUE, 0
303  );
304  pfree(sql);
305 
306  /* process resultset */
307  SPI_cursor_fetch(portal, TRUE, 1);
308  while (SPI_processed == 1 && SPI_tuptable != NULL) {
309  tupdesc = SPI_tuptable->tupdesc;
310  tuple = SPI_tuptable->vals[0];
311 
312  datum = SPI_getbinval(tuple, tupdesc, 1, &isNull);
313  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
314  SPI_freetuptable(SPI_tuptable);
315  SPI_cursor_close(portal);
316  SPI_finish();
317 
318  if (NULL != rtn) pfree(rtn);
319  elog(ERROR, "RASTER_summaryStatsCoverage: Cannot get raster of coverage");
320  PG_RETURN_NULL();
321  }
322  else if (isNull) {
323  SPI_cursor_fetch(portal, TRUE, 1);
324  continue;
325  }
326 
327  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(datum);
328 
329  raster = rt_raster_deserialize(pgraster, FALSE);
330  if (!raster) {
331  SPI_freetuptable(SPI_tuptable);
332  SPI_cursor_close(portal);
333  SPI_finish();
334 
335  if (NULL != rtn) pfree(rtn);
336  elog(ERROR, "RASTER_summaryStatsCoverage: Cannot deserialize raster");
337  PG_RETURN_NULL();
338  }
339 
340  /* inspect number of bands */
341  num_bands = rt_raster_get_num_bands(raster);
342  if (bandindex < 1 || bandindex > num_bands) {
343  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
344 
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 
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 
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:262
#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:2387
struct rt_bandstats_t * rt_bandstats
Definition: librtcore.h:152
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: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.
Definition: rt_serialize.c:725
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:385
band
Definition: ovdump.py:58
raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
#define VALUES_LENGTH
uint32_t count
Definition: librtcore.h:2513
double * values
Definition: librtcore.h:2521
Struct definitions.
Definition: librtcore.h:2403

References ovdump::band, rt_bandstats_t::count, FALSE, FLT_EQ, rt_bandstats_t::max, rt_bandstats_t::mean, rt_bandstats_t::min, rtrowdump::raster, 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, rtgdalraster::sql, rt_bandstats_t::stddev, rt_bandstats_t::sum, TRUE, rt_bandstats_t::values, and VALUES_LENGTH.

Here is the call graph for this function: