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

Definition at line 9209 of file rt_pg.c.

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

9210 {
9211  FuncCallContext *funcctx;
9212  TupleDesc tupdesc;
9213 
9214  int i;
9215  rt_quantile quant;
9216  rt_quantile quant2;
9217  int call_cntr;
9218  int max_calls;
9219 
9220  /* first call of function */
9221  if (SRF_IS_FIRSTCALL()) {
9222  MemoryContext oldcontext;
9223 
9224  rt_pgraster *pgraster = NULL;
9225  rt_raster raster = NULL;
9226  rt_band band = NULL;
9227  int32_t bandindex = 0;
9228  int num_bands = 0;
9229  bool exclude_nodata_value = TRUE;
9230  double sample = 0;
9231  double *quantiles = NULL;
9232  uint32_t quantiles_count = 0;
9233  double quantile = 0;
9234  rt_bandstats stats = NULL;
9235  uint32_t count;
9236 
9237  int j;
9238  int n;
9239 
9240  ArrayType *array;
9241  Oid etype;
9242  Datum *e;
9243  bool *nulls;
9244  int16 typlen;
9245  bool typbyval;
9246  char typalign;
9247 
9248  /* create a function context for cross-call persistence */
9249  funcctx = SRF_FIRSTCALL_INIT();
9250 
9251  /* switch to memory context appropriate for multiple function calls */
9252  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
9253 
9254  /* pgraster is null, return nothing */
9255  if (PG_ARGISNULL(0)) {
9256  MemoryContextSwitchTo(oldcontext);
9257  SRF_RETURN_DONE(funcctx);
9258  }
9259  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
9260 
9261  raster = rt_raster_deserialize(pgraster, FALSE);
9262  if (!raster) {
9263  PG_FREE_IF_COPY(pgraster, 0);
9264  MemoryContextSwitchTo(oldcontext);
9265  elog(ERROR, "RASTER_quantile: Could not deserialize raster");
9266  SRF_RETURN_DONE(funcctx);
9267  }
9268 
9269  /* band index is 1-based */
9270  bandindex = PG_GETARG_INT32(1);
9271  num_bands = rt_raster_get_num_bands(raster);
9272  if (bandindex < 1 || bandindex > num_bands) {
9273  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
9274  rt_raster_destroy(raster);
9275  PG_FREE_IF_COPY(pgraster, 0);
9276  MemoryContextSwitchTo(oldcontext);
9277  SRF_RETURN_DONE(funcctx);
9278  }
9279 
9280  /* exclude_nodata_value flag */
9281  if (!PG_ARGISNULL(2))
9282  exclude_nodata_value = PG_GETARG_BOOL(2);
9283 
9284  /* sample % */
9285  if (!PG_ARGISNULL(3)) {
9286  sample = PG_GETARG_FLOAT8(3);
9287  if (sample < 0 || sample > 1) {
9288  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
9289  rt_raster_destroy(raster);
9290  PG_FREE_IF_COPY(pgraster, 0);
9291  MemoryContextSwitchTo(oldcontext);
9292  SRF_RETURN_DONE(funcctx);
9293  }
9294  else if (FLT_EQ(sample, 0.0))
9295  sample = 1;
9296  }
9297  else
9298  sample = 1;
9299 
9300  /* quantiles */
9301  if (!PG_ARGISNULL(4)) {
9302  array = PG_GETARG_ARRAYTYPE_P(4);
9303  etype = ARR_ELEMTYPE(array);
9304  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
9305 
9306  switch (etype) {
9307  case FLOAT4OID:
9308  case FLOAT8OID:
9309  break;
9310  default:
9311  rt_raster_destroy(raster);
9312  PG_FREE_IF_COPY(pgraster, 0);
9313  MemoryContextSwitchTo(oldcontext);
9314  elog(ERROR, "RASTER_quantile: Invalid data type for quantiles");
9315  SRF_RETURN_DONE(funcctx);
9316  break;
9317  }
9318 
9319  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
9320  &nulls, &n);
9321 
9322  quantiles = palloc(sizeof(double) * n);
9323  for (i = 0, j = 0; i < n; i++) {
9324  if (nulls[i]) continue;
9325 
9326  switch (etype) {
9327  case FLOAT4OID:
9328  quantile = (double) DatumGetFloat4(e[i]);
9329  break;
9330  case FLOAT8OID:
9331  quantile = (double) DatumGetFloat8(e[i]);
9332  break;
9333  }
9334 
9335  if (quantile < 0 || quantile > 1) {
9336  elog(NOTICE, "Invalid value for quantile (must be between 0 and 1). Returning NULL");
9337  pfree(quantiles);
9338  rt_raster_destroy(raster);
9339  PG_FREE_IF_COPY(pgraster, 0);
9340  MemoryContextSwitchTo(oldcontext);
9341  SRF_RETURN_DONE(funcctx);
9342  }
9343 
9344  quantiles[j] = quantile;
9345  POSTGIS_RT_DEBUGF(5, "quantiles[%d] = %f", j, quantiles[j]);
9346  j++;
9347  }
9348  quantiles_count = j;
9349 
9350  if (j < 1) {
9351  pfree(quantiles);
9352  quantiles = NULL;
9353  }
9354  }
9355 
9356  /* get band */
9357  band = rt_raster_get_band(raster, bandindex - 1);
9358  if (!band) {
9359  elog(NOTICE, "Could not find band at index %d. Returning NULL", bandindex);
9360  rt_raster_destroy(raster);
9361  PG_FREE_IF_COPY(pgraster, 0);
9362  MemoryContextSwitchTo(oldcontext);
9363  SRF_RETURN_DONE(funcctx);
9364  }
9365 
9366  /* get stats */
9367  stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 1, NULL, NULL, NULL);
9368  rt_band_destroy(band);
9369  rt_raster_destroy(raster);
9370  PG_FREE_IF_COPY(pgraster, 0);
9371  if (NULL == stats || NULL == stats->values) {
9372  elog(NOTICE, "Could not retrieve summary statistics for band at index %d", bandindex);
9373  MemoryContextSwitchTo(oldcontext);
9374  SRF_RETURN_DONE(funcctx);
9375  }
9376  else if (stats->count < 1) {
9377  elog(NOTICE, "Could not compute quantiles for band at index %d as the band has no values", bandindex);
9378  MemoryContextSwitchTo(oldcontext);
9379  SRF_RETURN_DONE(funcctx);
9380  }
9381 
9382  /* get quantiles */
9383  quant = rt_band_get_quantiles(stats, quantiles, quantiles_count, &count);
9384  if (quantiles_count) pfree(quantiles);
9385  pfree(stats);
9386  if (NULL == quant || !count) {
9387  elog(NOTICE, "Could not compute quantiles for band at index %d", bandindex);
9388  MemoryContextSwitchTo(oldcontext);
9389  SRF_RETURN_DONE(funcctx);
9390  }
9391 
9392  POSTGIS_RT_DEBUGF(3, "%d quantiles returned", count);
9393 
9394  /* Store needed information */
9395  funcctx->user_fctx = quant;
9396 
9397  /* total number of tuples to be returned */
9398  funcctx->max_calls = count;
9399 
9400  /* Build a tuple descriptor for our result type */
9401  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
9402  ereport(ERROR, (
9403  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
9404  errmsg(
9405  "function returning record called in context "
9406  "that cannot accept type record"
9407  )
9408  ));
9409  }
9410 
9411  BlessTupleDesc(tupdesc);
9412  funcctx->tuple_desc = tupdesc;
9413 
9414  MemoryContextSwitchTo(oldcontext);
9415  }
9416 
9417  /* stuff done on every call of the function */
9418  funcctx = SRF_PERCALL_SETUP();
9419 
9420  call_cntr = funcctx->call_cntr;
9421  max_calls = funcctx->max_calls;
9422  tupdesc = funcctx->tuple_desc;
9423  quant2 = funcctx->user_fctx;
9424 
9425  /* do when there is more left to send */
9426  if (call_cntr < max_calls) {
9427  int values_length = 2;
9428  Datum values[values_length];
9429  bool nulls[values_length];
9430  HeapTuple tuple;
9431  Datum result;
9432 
9433  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
9434 
9435  memset(nulls, FALSE, sizeof(bool) * values_length);
9436 
9437  values[0] = Float8GetDatum(quant2[call_cntr].quantile);
9438  values[1] = Float8GetDatum(quant2[call_cntr].value);
9439 
9440  /* build a tuple */
9441  tuple = heap_form_tuple(tupdesc, values, nulls);
9442 
9443  /* make the tuple into a datum */
9444  result = HeapTupleGetDatum(tuple);
9445 
9446  SRF_RETURN_NEXT(funcctx, result);
9447  }
9448  /* do when there is no more left */
9449  else {
9450  pfree(quant2);
9451  SRF_RETURN_DONE(funcctx);
9452  }
9453 }
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
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
rt_quantile rt_band_get_quantiles(rt_bandstats stats, double *quantiles, int quantiles_count, uint32_t *rtn_count)
Compute the default set of or requested quantiles for a set of data the quantile formula used is same...
Definition: rt_api.c:3823
char ** result
Definition: liblwgeom.h:218
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: