PostGIS  2.1.10dev-r@@SVN_REVISION@@
Datum RASTER_histogram ( PG_FUNCTION_ARGS  )

Definition at line 8403 of file rt_pg.c.

References ovdump::band, genraster::count, rt_bandstats_t::count, FALSE, FLT_EQ, POSTGIS_RT_DEBUG, POSTGIS_RT_DEBUGF, rtrowdump::raster, result, rt_band_destroy(), rt_band_get_histogram(), rt_band_get_summary_stats(), rt_raster_deserialize(), rt_raster_destroy(), rt_raster_get_band(), rt_raster_get_num_bands(), TRUE, rt_bandstats_t::values, and rt_raster_serialized_t::width.

8404 {
8405  FuncCallContext *funcctx;
8406  TupleDesc tupdesc;
8407 
8408  int i;
8409  rt_histogram hist;
8410  rt_histogram hist2;
8411  int call_cntr;
8412  int max_calls;
8413 
8414  /* first call of function */
8415  if (SRF_IS_FIRSTCALL()) {
8416  MemoryContext oldcontext;
8417 
8418  rt_pgraster *pgraster = NULL;
8419  rt_raster raster = NULL;
8420  rt_band band = NULL;
8421  int32_t bandindex = 1;
8422  int num_bands = 0;
8423  bool exclude_nodata_value = TRUE;
8424  double sample = 0;
8425  uint32_t bin_count = 0;
8426  double *bin_width = NULL;
8427  uint32_t bin_width_count = 0;
8428  double width = 0;
8429  bool right = FALSE;
8430  double min = 0;
8431  double max = 0;
8432  rt_bandstats stats = NULL;
8433  uint32_t count;
8434 
8435  int j;
8436  int n;
8437 
8438  ArrayType *array;
8439  Oid etype;
8440  Datum *e;
8441  bool *nulls;
8442  int16 typlen;
8443  bool typbyval;
8444  char typalign;
8445 
8446  POSTGIS_RT_DEBUG(3, "RASTER_histogram: Starting");
8447 
8448  /* create a function context for cross-call persistence */
8449  funcctx = SRF_FIRSTCALL_INIT();
8450 
8451  /* switch to memory context appropriate for multiple function calls */
8452  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
8453 
8454  /* pgraster is null, return nothing */
8455  if (PG_ARGISNULL(0)) {
8456  MemoryContextSwitchTo(oldcontext);
8457  SRF_RETURN_DONE(funcctx);
8458  }
8459  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
8460 
8461  raster = rt_raster_deserialize(pgraster, FALSE);
8462  if (!raster) {
8463  PG_FREE_IF_COPY(pgraster, 0);
8464  MemoryContextSwitchTo(oldcontext);
8465  elog(ERROR, "RASTER_histogram: Could not deserialize raster");
8466  SRF_RETURN_DONE(funcctx);
8467  }
8468 
8469  /* band index is 1-based */
8470  if (!PG_ARGISNULL(1))
8471  bandindex = PG_GETARG_INT32(1);
8472  num_bands = rt_raster_get_num_bands(raster);
8473  if (bandindex < 1 || bandindex > num_bands) {
8474  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
8475  rt_raster_destroy(raster);
8476  PG_FREE_IF_COPY(pgraster, 0);
8477  MemoryContextSwitchTo(oldcontext);
8478  SRF_RETURN_DONE(funcctx);
8479  }
8480 
8481  /* exclude_nodata_value flag */
8482  if (!PG_ARGISNULL(2))
8483  exclude_nodata_value = PG_GETARG_BOOL(2);
8484 
8485  /* sample % */
8486  if (!PG_ARGISNULL(3)) {
8487  sample = PG_GETARG_FLOAT8(3);
8488  if (sample < 0 || sample > 1) {
8489  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
8490  rt_raster_destroy(raster);
8491  PG_FREE_IF_COPY(pgraster, 0);
8492  MemoryContextSwitchTo(oldcontext);
8493  SRF_RETURN_DONE(funcctx);
8494  }
8495  else if (FLT_EQ(sample, 0.0))
8496  sample = 1;
8497  }
8498  else
8499  sample = 1;
8500 
8501  /* bin_count */
8502  if (!PG_ARGISNULL(4)) {
8503  bin_count = PG_GETARG_INT32(4);
8504  if (bin_count < 1) bin_count = 0;
8505  }
8506 
8507  /* bin_width */
8508  if (!PG_ARGISNULL(5)) {
8509  array = PG_GETARG_ARRAYTYPE_P(5);
8510  etype = ARR_ELEMTYPE(array);
8511  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
8512 
8513  switch (etype) {
8514  case FLOAT4OID:
8515  case FLOAT8OID:
8516  break;
8517  default:
8518  rt_raster_destroy(raster);
8519  PG_FREE_IF_COPY(pgraster, 0);
8520  MemoryContextSwitchTo(oldcontext);
8521  elog(ERROR, "RASTER_histogram: Invalid data type for width");
8522  SRF_RETURN_DONE(funcctx);
8523  break;
8524  }
8525 
8526  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
8527  &nulls, &n);
8528 
8529  bin_width = palloc(sizeof(double) * n);
8530  for (i = 0, j = 0; i < n; i++) {
8531  if (nulls[i]) continue;
8532 
8533  switch (etype) {
8534  case FLOAT4OID:
8535  width = (double) DatumGetFloat4(e[i]);
8536  break;
8537  case FLOAT8OID:
8538  width = (double) DatumGetFloat8(e[i]);
8539  break;
8540  }
8541 
8542  if (width < 0 || FLT_EQ(width, 0.0)) {
8543  elog(NOTICE, "Invalid value for width (must be greater than 0). Returning NULL");
8544  pfree(bin_width);
8545  rt_raster_destroy(raster);
8546  PG_FREE_IF_COPY(pgraster, 0);
8547  MemoryContextSwitchTo(oldcontext);
8548  SRF_RETURN_DONE(funcctx);
8549  }
8550 
8551  bin_width[j] = width;
8552  POSTGIS_RT_DEBUGF(5, "bin_width[%d] = %f", j, bin_width[j]);
8553  j++;
8554  }
8555  bin_width_count = j;
8556 
8557  if (j < 1) {
8558  pfree(bin_width);
8559  bin_width = NULL;
8560  }
8561  }
8562 
8563  /* right */
8564  if (!PG_ARGISNULL(6))
8565  right = PG_GETARG_BOOL(6);
8566 
8567  /* min */
8568  if (!PG_ARGISNULL(7)) min = PG_GETARG_FLOAT8(7);
8569 
8570  /* max */
8571  if (!PG_ARGISNULL(8)) max = PG_GETARG_FLOAT8(8);
8572 
8573  /* get band */
8574  band = rt_raster_get_band(raster, bandindex - 1);
8575  if (!band) {
8576  elog(NOTICE, "Could not find band at index %d. Returning NULL", bandindex);
8577  rt_raster_destroy(raster);
8578  PG_FREE_IF_COPY(pgraster, 0);
8579  MemoryContextSwitchTo(oldcontext);
8580  SRF_RETURN_DONE(funcctx);
8581  }
8582 
8583  /* get stats */
8584  stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 1, NULL, NULL, NULL);
8585  rt_band_destroy(band);
8586  rt_raster_destroy(raster);
8587  PG_FREE_IF_COPY(pgraster, 0);
8588  if (NULL == stats || NULL == stats->values) {
8589  elog(NOTICE, "Could not compute summary statistics for band at index %d", bandindex);
8590  MemoryContextSwitchTo(oldcontext);
8591  SRF_RETURN_DONE(funcctx);
8592  }
8593  else if (stats->count < 1) {
8594  elog(NOTICE, "Could not compute histogram for band at index %d as the band has no values", bandindex);
8595  MemoryContextSwitchTo(oldcontext);
8596  SRF_RETURN_DONE(funcctx);
8597  }
8598 
8599  /* get histogram */
8600  hist = rt_band_get_histogram(stats, bin_count, bin_width, bin_width_count, right, min, max, &count);
8601  if (bin_width_count) pfree(bin_width);
8602  pfree(stats);
8603  if (NULL == hist || !count) {
8604  elog(NOTICE, "Could not compute histogram for band at index %d", bandindex);
8605  MemoryContextSwitchTo(oldcontext);
8606  SRF_RETURN_DONE(funcctx);
8607  }
8608 
8609  POSTGIS_RT_DEBUGF(3, "%d bins returned", count);
8610 
8611  /* Store needed information */
8612  funcctx->user_fctx = hist;
8613 
8614  /* total number of tuples to be returned */
8615  funcctx->max_calls = count;
8616 
8617  /* Build a tuple descriptor for our result type */
8618  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
8619  ereport(ERROR, (
8620  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
8621  errmsg(
8622  "function returning record called in context "
8623  "that cannot accept type record"
8624  )
8625  ));
8626  }
8627 
8628  BlessTupleDesc(tupdesc);
8629  funcctx->tuple_desc = tupdesc;
8630 
8631  MemoryContextSwitchTo(oldcontext);
8632  }
8633 
8634  /* stuff done on every call of the function */
8635  funcctx = SRF_PERCALL_SETUP();
8636 
8637  call_cntr = funcctx->call_cntr;
8638  max_calls = funcctx->max_calls;
8639  tupdesc = funcctx->tuple_desc;
8640  hist2 = funcctx->user_fctx;
8641 
8642  /* do when there is more left to send */
8643  if (call_cntr < max_calls) {
8644  int values_length = 4;
8645  Datum values[values_length];
8646  bool nulls[values_length];
8647  HeapTuple tuple;
8648  Datum result;
8649 
8650  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
8651 
8652  memset(nulls, FALSE, sizeof(bool) * values_length);
8653 
8654  values[0] = Float8GetDatum(hist2[call_cntr].min);
8655  values[1] = Float8GetDatum(hist2[call_cntr].max);
8656  values[2] = Int64GetDatum(hist2[call_cntr].count);
8657  values[3] = Float8GetDatum(hist2[call_cntr].percent);
8658 
8659  /* build a tuple */
8660  tuple = heap_form_tuple(tupdesc, values, nulls);
8661 
8662  /* make the tuple into a datum */
8663  result = HeapTupleGetDatum(tuple);
8664 
8665  SRF_RETURN_NEXT(funcctx, result);
8666  }
8667  /* do when there is no more left */
8668  else {
8669  pfree(hist2);
8670  SRF_RETURN_DONE(funcctx);
8671  }
8672 }
int rt_raster_get_num_bands(rt_raster raster)
Definition: rt_api.c:5677
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_api.c:5387
rt_histogram rt_band_get_histogram(rt_bandstats stats, int bin_count, double *bin_width, int bin_width_count, int right, double min, double max, uint32_t *rtn_count)
Count the distribution of data.
Definition: rt_api.c:3567
uint32_t count
Definition: rt_api.h:2277
tuple band
Definition: ovdump.py:57
tuple raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:123
char ** result
Definition: liblwgeom.h:218
#define POSTGIS_RT_DEBUG(level, msg)
Definition: rt_pg.h:58
int count
Definition: genraster.py:57
double * values
Definition: rt_api.h:2285
#define FLT_EQ(x, y)
Definition: rt_api.h:2159
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition: rt_api.c:1650
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_api.c:3269
rt_band rt_raster_get_band(rt_raster raster, int n)
Return Nth band, or NULL if unavailable.
Definition: rt_api.c:5686
#define FALSE
Definition: dbfopen.c:169
Struct definitions.
Definition: rt_api.h:2175
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition: rt_pg.h:62
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition: rt_api.c:8350
#define TRUE
Definition: dbfopen.c:170

Here is the call graph for this function: