PostGIS  2.5.7dev-r@@SVN_REVISION@@

◆ 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");
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 
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 
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 
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:2234
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.
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:381
band
Definition: ovdump.py:57
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:2360
double * values
Definition: librtcore.h:2368
Struct definitions.
Definition: librtcore.h:2250

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

Here is the call graph for this function: