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

Definition at line 9459 of file rt_pg.c.

References ovdump::band, genraster::count, FALSE, FLT_EQ, rt_quantile_t::has_value, MAX_DBL_CHARLEN, MAX_INT_CHARLEN, POSTGIS_RT_DEBUG, POSTGIS_RT_DEBUGF, rt_quantile_t::quantile, quantile_llist::quantile, quantile_llist_destroy(), rtrowdump::raster, result, rt_band_destroy(), rt_band_get_quantiles_stream(), rt_raster_deserialize(), rt_raster_destroy(), rt_raster_get_band(), rt_raster_get_num_bands(), rtgdalraster::sql, TRUE, genraster::value, and rt_quantile_t::value.

9460 {
9461  FuncCallContext *funcctx;
9462  TupleDesc tupdesc;
9463 
9464  int i;
9465  rt_quantile covquant = NULL;
9466  rt_quantile covquant2;
9467  int call_cntr;
9468  int max_calls;
9469 
9470  POSTGIS_RT_DEBUG(3, "RASTER_quantileCoverage: Starting");
9471 
9472  /* first call of function */
9473  if (SRF_IS_FIRSTCALL()) {
9474  MemoryContext oldcontext;
9475 
9476  text *tablenametext = NULL;
9477  char *tablename = NULL;
9478  text *colnametext = NULL;
9479  char *colname = NULL;
9480  int32_t bandindex = 1;
9481  bool exclude_nodata_value = TRUE;
9482  double sample = 0;
9483  double *quantiles = NULL;
9484  uint32_t quantiles_count = 0;
9485  double quantile = 0;
9486  uint32_t count;
9487 
9488  int len = 0;
9489  char *sql = NULL;
9490  char *tmp = NULL;
9491  uint64_t cov_count = 0;
9492  int spi_result;
9493  Portal portal;
9494  SPITupleTable *tuptable = NULL;
9495  HeapTuple tuple;
9496  Datum datum;
9497  bool isNull = FALSE;
9498 
9499  rt_pgraster *pgraster = NULL;
9500  rt_raster raster = NULL;
9501  rt_band band = NULL;
9502  int num_bands = 0;
9503  struct quantile_llist *qlls = NULL;
9504  uint32_t qlls_count;
9505 
9506  int j;
9507  int n;
9508 
9509  ArrayType *array;
9510  Oid etype;
9511  Datum *e;
9512  bool *nulls;
9513  int16 typlen;
9514  bool typbyval;
9515  char typalign;
9516 
9517  POSTGIS_RT_DEBUG(3, "RASTER_quantileCoverage: first call of function");
9518 
9519  /* create a function context for cross-call persistence */
9520  funcctx = SRF_FIRSTCALL_INIT();
9521 
9522  /* switch to memory context appropriate for multiple function calls */
9523  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
9524 
9525  /* tablename is null, return null */
9526  if (PG_ARGISNULL(0)) {
9527  elog(NOTICE, "Table name must be provided");
9528  MemoryContextSwitchTo(oldcontext);
9529  SRF_RETURN_DONE(funcctx);
9530  }
9531  tablenametext = PG_GETARG_TEXT_P(0);
9532  tablename = text_to_cstring(tablenametext);
9533  if (!strlen(tablename)) {
9534  elog(NOTICE, "Table name must be provided");
9535  MemoryContextSwitchTo(oldcontext);
9536  SRF_RETURN_DONE(funcctx);
9537  }
9538  POSTGIS_RT_DEBUGF(3, "RASTER_quantileCoverage: tablename = %s", tablename);
9539 
9540  /* column name is null, return null */
9541  if (PG_ARGISNULL(1)) {
9542  elog(NOTICE, "Column name must be provided");
9543  MemoryContextSwitchTo(oldcontext);
9544  SRF_RETURN_DONE(funcctx);
9545  }
9546  colnametext = PG_GETARG_TEXT_P(1);
9547  colname = text_to_cstring(colnametext);
9548  if (!strlen(colname)) {
9549  elog(NOTICE, "Column name must be provided");
9550  MemoryContextSwitchTo(oldcontext);
9551  SRF_RETURN_DONE(funcctx);
9552  }
9553  POSTGIS_RT_DEBUGF(3, "RASTER_quantileCoverage: colname = %s", colname);
9554 
9555  /* band index is 1-based */
9556  if (!PG_ARGISNULL(2))
9557  bandindex = PG_GETARG_INT32(2);
9558 
9559  /* exclude_nodata_value flag */
9560  if (!PG_ARGISNULL(3))
9561  exclude_nodata_value = PG_GETARG_BOOL(3);
9562 
9563  /* sample % */
9564  if (!PG_ARGISNULL(4)) {
9565  sample = PG_GETARG_FLOAT8(4);
9566  if (sample < 0 || sample > 1) {
9567  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
9568  MemoryContextSwitchTo(oldcontext);
9569  SRF_RETURN_DONE(funcctx);
9570  }
9571  else if (FLT_EQ(sample, 0.0))
9572  sample = 1;
9573  }
9574  else
9575  sample = 1;
9576 
9577  /* quantiles */
9578  if (!PG_ARGISNULL(5)) {
9579  array = PG_GETARG_ARRAYTYPE_P(5);
9580  etype = ARR_ELEMTYPE(array);
9581  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
9582 
9583  switch (etype) {
9584  case FLOAT4OID:
9585  case FLOAT8OID:
9586  break;
9587  default:
9588  MemoryContextSwitchTo(oldcontext);
9589  elog(ERROR, "RASTER_quantileCoverage: Invalid data type for quantiles");
9590  SRF_RETURN_DONE(funcctx);
9591  break;
9592  }
9593 
9594  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
9595  &nulls, &n);
9596 
9597  quantiles = palloc(sizeof(double) * n);
9598  for (i = 0, j = 0; i < n; i++) {
9599  if (nulls[i]) continue;
9600 
9601  switch (etype) {
9602  case FLOAT4OID:
9603  quantile = (double) DatumGetFloat4(e[i]);
9604  break;
9605  case FLOAT8OID:
9606  quantile = (double) DatumGetFloat8(e[i]);
9607  break;
9608  }
9609 
9610  if (quantile < 0 || quantile > 1) {
9611  elog(NOTICE, "Invalid value for quantile (must be between 0 and 1). Returning NULL");
9612  pfree(quantiles);
9613  MemoryContextSwitchTo(oldcontext);
9614  SRF_RETURN_DONE(funcctx);
9615  }
9616 
9617  quantiles[j] = quantile;
9618  POSTGIS_RT_DEBUGF(5, "quantiles[%d] = %f", j, quantiles[j]);
9619  j++;
9620  }
9621  quantiles_count = j;
9622 
9623  if (j < 1) {
9624  pfree(quantiles);
9625  quantiles = NULL;
9626  }
9627  }
9628 
9629  /* coverage stats */
9630  /* connect to database */
9631  spi_result = SPI_connect();
9632  if (spi_result != SPI_OK_CONNECT) {
9633  MemoryContextSwitchTo(oldcontext);
9634  elog(ERROR, "RASTER_quantileCoverage: Could not connect to database using SPI");
9635  SRF_RETURN_DONE(funcctx);
9636  }
9637 
9638  len = sizeof(char) * (strlen("SELECT count FROM _st_summarystats('','',,::boolean,)") + strlen(tablename) + strlen(colname) + (MAX_INT_CHARLEN * 2) + MAX_DBL_CHARLEN + 1);
9639  sql = (char *) palloc(len);
9640  if (NULL == sql) {
9641 
9642  if (SPI_tuptable) SPI_freetuptable(tuptable);
9643  SPI_finish();
9644 
9645  MemoryContextSwitchTo(oldcontext);
9646  elog(ERROR, "RASTER_quantileCoverage: Could not allocate memory for sql");
9647  SRF_RETURN_DONE(funcctx);
9648  }
9649 
9650  /* get stats */
9651  snprintf(sql, len, "SELECT count FROM _st_summarystats('%s','%s',%d,%d::boolean,%f)", tablename, colname, bandindex, (exclude_nodata_value ? 1 : 0), sample);
9652  POSTGIS_RT_DEBUGF(3, "stats sql: %s", sql);
9653  spi_result = SPI_execute(sql, TRUE, 0);
9654  pfree(sql);
9655  if (spi_result != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
9656 
9657  if (SPI_tuptable) SPI_freetuptable(tuptable);
9658  SPI_finish();
9659 
9660  MemoryContextSwitchTo(oldcontext);
9661  elog(ERROR, "RASTER_quantileCoverage: Could not get summary stats of coverage");
9662  SRF_RETURN_DONE(funcctx);
9663  }
9664 
9665  tupdesc = SPI_tuptable->tupdesc;
9666  tuptable = SPI_tuptable;
9667  tuple = tuptable->vals[0];
9668 
9669  tmp = SPI_getvalue(tuple, tupdesc, 1);
9670  if (NULL == tmp || !strlen(tmp)) {
9671 
9672  if (SPI_tuptable) SPI_freetuptable(tuptable);
9673  SPI_finish();
9674 
9675  MemoryContextSwitchTo(oldcontext);
9676  elog(ERROR, "RASTER_quantileCoverage: Could not get summary stats of coverage");
9677  SRF_RETURN_DONE(funcctx);
9678  }
9679  cov_count = strtol(tmp, NULL, 10);
9680  POSTGIS_RT_DEBUGF(3, "covcount = %d", (int) cov_count);
9681  pfree(tmp);
9682 
9683  /* iterate through rasters of coverage */
9684  /* create sql */
9685  len = sizeof(char) * (strlen("SELECT \"\" FROM \"\" WHERE \"\" IS NOT NULL") + (strlen(colname) * 2) + strlen(tablename) + 1);
9686  sql = (char *) palloc(len);
9687  if (NULL == sql) {
9688 
9689  if (SPI_tuptable) SPI_freetuptable(tuptable);
9690  SPI_finish();
9691 
9692  MemoryContextSwitchTo(oldcontext);
9693  elog(ERROR, "RASTER_quantileCoverage: Could not allocate memory for sql");
9694  SRF_RETURN_DONE(funcctx);
9695  }
9696 
9697  /* get cursor */
9698  snprintf(sql, len, "SELECT \"%s\" FROM \"%s\" WHERE \"%s\" IS NOT NULL", colname, tablename, colname);
9699  POSTGIS_RT_DEBUGF(3, "coverage sql: %s", sql);
9700  portal = SPI_cursor_open_with_args(
9701  "coverage",
9702  sql,
9703  0, NULL,
9704  NULL, NULL,
9705  TRUE, 0
9706  );
9707  pfree(sql);
9708 
9709  /* process resultset */
9710  SPI_cursor_fetch(portal, TRUE, 1);
9711  while (SPI_processed == 1 && SPI_tuptable != NULL) {
9712  if (NULL != covquant) pfree(covquant);
9713 
9714  tupdesc = SPI_tuptable->tupdesc;
9715  tuptable = SPI_tuptable;
9716  tuple = tuptable->vals[0];
9717 
9718  datum = SPI_getbinval(tuple, tupdesc, 1, &isNull);
9719  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
9720 
9721  if (SPI_tuptable) SPI_freetuptable(tuptable);
9722  SPI_cursor_close(portal);
9723  SPI_finish();
9724 
9725  MemoryContextSwitchTo(oldcontext);
9726  elog(ERROR, "RASTER_quantileCoverage: Could not get raster of coverage");
9727  SRF_RETURN_DONE(funcctx);
9728  }
9729  else if (isNull) {
9730  SPI_cursor_fetch(portal, TRUE, 1);
9731  continue;
9732  }
9733 
9734  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(datum);
9735 
9736  raster = rt_raster_deserialize(pgraster, FALSE);
9737  if (!raster) {
9738 
9739  if (SPI_tuptable) SPI_freetuptable(tuptable);
9740  SPI_cursor_close(portal);
9741  SPI_finish();
9742 
9743  MemoryContextSwitchTo(oldcontext);
9744  elog(ERROR, "RASTER_quantileCoverage: Could not deserialize raster");
9745  SRF_RETURN_DONE(funcctx);
9746  }
9747 
9748  /* inspect number of bands*/
9749  num_bands = rt_raster_get_num_bands(raster);
9750  if (bandindex < 1 || bandindex > num_bands) {
9751  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
9752 
9753  rt_raster_destroy(raster);
9754 
9755  if (SPI_tuptable) SPI_freetuptable(tuptable);
9756  SPI_cursor_close(portal);
9757  SPI_finish();
9758 
9759  MemoryContextSwitchTo(oldcontext);
9760  SRF_RETURN_DONE(funcctx);
9761  }
9762 
9763  /* get band */
9764  band = rt_raster_get_band(raster, bandindex - 1);
9765  if (!band) {
9766  elog(NOTICE, "Could not find raster band of index %d. Returning NULL", bandindex);
9767 
9768  rt_raster_destroy(raster);
9769 
9770  if (SPI_tuptable) SPI_freetuptable(tuptable);
9771  SPI_cursor_close(portal);
9772  SPI_finish();
9773 
9774  MemoryContextSwitchTo(oldcontext);
9775  SRF_RETURN_DONE(funcctx);
9776  }
9777 
9778  covquant = rt_band_get_quantiles_stream(
9779  band,
9780  exclude_nodata_value, sample, cov_count,
9781  &qlls, &qlls_count,
9782  quantiles, quantiles_count,
9783  &count
9784  );
9785 
9786  rt_band_destroy(band);
9787  rt_raster_destroy(raster);
9788 
9789  if (NULL == covquant || !count) {
9790  elog(NOTICE, "Could not compute quantiles for band at index %d", bandindex);
9791 
9792  if (SPI_tuptable) SPI_freetuptable(tuptable);
9793  SPI_cursor_close(portal);
9794  SPI_finish();
9795 
9796  MemoryContextSwitchTo(oldcontext);
9797  SRF_RETURN_DONE(funcctx);
9798  }
9799 
9800  /* next record */
9801  SPI_cursor_fetch(portal, TRUE, 1);
9802  }
9803 
9804  covquant2 = SPI_palloc(sizeof(struct rt_quantile_t) * count);
9805  for (i = 0; i < count; i++) {
9806  covquant2[i].quantile = covquant[i].quantile;
9807  covquant2[i].has_value = covquant[i].has_value;
9808  if (covquant2[i].has_value)
9809  covquant2[i].value = covquant[i].value;
9810  }
9811 
9812  if (NULL != covquant) pfree(covquant);
9813  quantile_llist_destroy(&qlls, qlls_count);
9814 
9815  if (SPI_tuptable) SPI_freetuptable(tuptable);
9816  SPI_cursor_close(portal);
9817  SPI_finish();
9818 
9819  if (quantiles_count) pfree(quantiles);
9820 
9821  POSTGIS_RT_DEBUGF(3, "%d quantiles returned", count);
9822 
9823  /* Store needed information */
9824  funcctx->user_fctx = covquant2;
9825 
9826  /* total number of tuples to be returned */
9827  funcctx->max_calls = count;
9828 
9829  /* Build a tuple descriptor for our result type */
9830  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
9831  ereport(ERROR, (
9832  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
9833  errmsg(
9834  "function returning record called in context "
9835  "that cannot accept type record"
9836  )
9837  ));
9838  }
9839 
9840  BlessTupleDesc(tupdesc);
9841  funcctx->tuple_desc = tupdesc;
9842 
9843  MemoryContextSwitchTo(oldcontext);
9844  }
9845 
9846  /* stuff done on every call of the function */
9847  funcctx = SRF_PERCALL_SETUP();
9848 
9849  call_cntr = funcctx->call_cntr;
9850  max_calls = funcctx->max_calls;
9851  tupdesc = funcctx->tuple_desc;
9852  covquant2 = funcctx->user_fctx;
9853 
9854  /* do when there is more left to send */
9855  if (call_cntr < max_calls) {
9856  int values_length = 2;
9857  Datum values[values_length];
9858  bool nulls[values_length];
9859  HeapTuple tuple;
9860  Datum result;
9861 
9862  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
9863 
9864  memset(nulls, FALSE, sizeof(bool) * values_length);
9865 
9866  values[0] = Float8GetDatum(covquant2[call_cntr].quantile);
9867  if (covquant2[call_cntr].has_value)
9868  values[1] = Float8GetDatum(covquant2[call_cntr].value);
9869  else
9870  nulls[1] = TRUE;
9871 
9872  /* build a tuple */
9873  tuple = heap_form_tuple(tupdesc, values, nulls);
9874 
9875  /* make the tuple into a datum */
9876  result = HeapTupleGetDatum(tuple);
9877 
9878  SRF_RETURN_NEXT(funcctx, result);
9879  }
9880  /* do when there is no more left */
9881  else {
9882  POSTGIS_RT_DEBUG(3, "done");
9883  pfree(covquant2);
9884  SRF_RETURN_DONE(funcctx);
9885  }
9886 }
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
double quantile
Definition: rt_api.h:2311
#define MAX_DBL_CHARLEN
Definition: rt_pg.c:64
double quantile
Definition: rt_api.h:2303
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
double value
Definition: rt_api.h:2304
#define POSTGIS_RT_DEBUG(level, msg)
Definition: rt_pg.h:58
int count
Definition: genraster.py:57
#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
int quantile_llist_destroy(struct quantile_llist **list, uint32_t list_count)
Definition: rt_api.c:4025
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
uint32_t has_value
Definition: rt_api.h:2305
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition: rt_api.c:8350
#define MAX_INT_CHARLEN
Definition: rt_pg.c:65
#define TRUE
Definition: dbfopen.c:170
rt_quantile rt_band_get_quantiles_stream(rt_band band, int exclude_nodata_value, double sample, uint64_t cov_count, struct quantile_llist **qlls, uint32_t *qlls_count, double *quantiles, int quantiles_count, uint32_t *rtn_count)
Compute the default set of or requested quantiles for a coverage.
Definition: rt_api.c:4156

Here is the call graph for this function: