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

Definition at line 8678 of file rt_pg.c.

References ovdump::band, genraster::count, rt_bandstats_t::count, rt_histogram_t::count, FALSE, FLT_EQ, rt_histogram_t::max, MAX_DBL_CHARLEN, MAX_INT_CHARLEN, rt_histogram_t::min, rt_histogram_t::percent, 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(), rtgdalraster::sql, TRUE, and rt_raster_serialized_t::width.

8679 {
8680  FuncCallContext *funcctx;
8681  TupleDesc tupdesc;
8682 
8683  int i;
8684  rt_histogram covhist = NULL;
8685  rt_histogram covhist2;
8686  int call_cntr;
8687  int max_calls;
8688 
8689  POSTGIS_RT_DEBUG(3, "RASTER_histogramCoverage: Starting");
8690 
8691  /* first call of function */
8692  if (SRF_IS_FIRSTCALL()) {
8693  MemoryContext oldcontext;
8694 
8695  text *tablenametext = NULL;
8696  char *tablename = NULL;
8697  text *colnametext = NULL;
8698  char *colname = NULL;
8699  int32_t bandindex = 1;
8700  bool exclude_nodata_value = TRUE;
8701  double sample = 0;
8702  uint32_t bin_count = 0;
8703  double *bin_width = NULL;
8704  uint32_t bin_width_count = 0;
8705  double width = 0;
8706  bool right = FALSE;
8707  uint32_t count;
8708 
8709  int len = 0;
8710  char *sql = NULL;
8711  char *tmp = NULL;
8712  double min = 0;
8713  double max = 0;
8714  int spi_result;
8715  Portal portal;
8716  SPITupleTable *tuptable = NULL;
8717  HeapTuple tuple;
8718  Datum datum;
8719  bool isNull = FALSE;
8720 
8721  rt_pgraster *pgraster = NULL;
8722  rt_raster raster = NULL;
8723  rt_band band = NULL;
8724  int num_bands = 0;
8725  rt_bandstats stats = NULL;
8726  rt_histogram hist;
8727  uint64_t sum = 0;
8728 
8729  int j;
8730  int n;
8731 
8732  ArrayType *array;
8733  Oid etype;
8734  Datum *e;
8735  bool *nulls;
8736  int16 typlen;
8737  bool typbyval;
8738  char typalign;
8739 
8740  POSTGIS_RT_DEBUG(3, "RASTER_histogramCoverage: first call of function");
8741 
8742  /* create a function context for cross-call persistence */
8743  funcctx = SRF_FIRSTCALL_INIT();
8744 
8745  /* switch to memory context appropriate for multiple function calls */
8746  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
8747 
8748  /* tablename is null, return null */
8749  if (PG_ARGISNULL(0)) {
8750  elog(NOTICE, "Table name must be provided");
8751  MemoryContextSwitchTo(oldcontext);
8752  SRF_RETURN_DONE(funcctx);
8753  }
8754  tablenametext = PG_GETARG_TEXT_P(0);
8755  tablename = text_to_cstring(tablenametext);
8756  if (!strlen(tablename)) {
8757  elog(NOTICE, "Table name must be provided");
8758  MemoryContextSwitchTo(oldcontext);
8759  SRF_RETURN_DONE(funcctx);
8760  }
8761  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: tablename = %s", tablename);
8762 
8763  /* column name is null, return null */
8764  if (PG_ARGISNULL(1)) {
8765  elog(NOTICE, "Column name must be provided");
8766  MemoryContextSwitchTo(oldcontext);
8767  SRF_RETURN_DONE(funcctx);
8768  }
8769  colnametext = PG_GETARG_TEXT_P(1);
8770  colname = text_to_cstring(colnametext);
8771  if (!strlen(colname)) {
8772  elog(NOTICE, "Column name must be provided");
8773  MemoryContextSwitchTo(oldcontext);
8774  SRF_RETURN_DONE(funcctx);
8775  }
8776  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: colname = %s", colname);
8777 
8778  /* band index is 1-based */
8779  if (!PG_ARGISNULL(2))
8780  bandindex = PG_GETARG_INT32(2);
8781 
8782  /* exclude_nodata_value flag */
8783  if (!PG_ARGISNULL(3))
8784  exclude_nodata_value = PG_GETARG_BOOL(3);
8785 
8786  /* sample % */
8787  if (!PG_ARGISNULL(4)) {
8788  sample = PG_GETARG_FLOAT8(4);
8789  if (sample < 0 || sample > 1) {
8790  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
8791  MemoryContextSwitchTo(oldcontext);
8792  SRF_RETURN_DONE(funcctx);
8793  }
8794  else if (FLT_EQ(sample, 0.0))
8795  sample = 1;
8796  }
8797  else
8798  sample = 1;
8799 
8800  /* bin_count */
8801  if (!PG_ARGISNULL(5)) {
8802  bin_count = PG_GETARG_INT32(5);
8803  if (bin_count < 1) bin_count = 0;
8804  }
8805 
8806  /* bin_width */
8807  if (!PG_ARGISNULL(6)) {
8808  array = PG_GETARG_ARRAYTYPE_P(6);
8809  etype = ARR_ELEMTYPE(array);
8810  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
8811 
8812  switch (etype) {
8813  case FLOAT4OID:
8814  case FLOAT8OID:
8815  break;
8816  default:
8817  MemoryContextSwitchTo(oldcontext);
8818  elog(ERROR, "RASTER_histogramCoverage: Invalid data type for width");
8819  SRF_RETURN_DONE(funcctx);
8820  break;
8821  }
8822 
8823  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
8824  &nulls, &n);
8825 
8826  bin_width = palloc(sizeof(double) * n);
8827  for (i = 0, j = 0; i < n; i++) {
8828  if (nulls[i]) continue;
8829 
8830  switch (etype) {
8831  case FLOAT4OID:
8832  width = (double) DatumGetFloat4(e[i]);
8833  break;
8834  case FLOAT8OID:
8835  width = (double) DatumGetFloat8(e[i]);
8836  break;
8837  }
8838 
8839  if (width < 0 || FLT_EQ(width, 0.0)) {
8840  elog(NOTICE, "Invalid value for width (must be greater than 0). Returning NULL");
8841  pfree(bin_width);
8842  MemoryContextSwitchTo(oldcontext);
8843  SRF_RETURN_DONE(funcctx);
8844  }
8845 
8846  bin_width[j] = width;
8847  POSTGIS_RT_DEBUGF(5, "bin_width[%d] = %f", j, bin_width[j]);
8848  j++;
8849  }
8850  bin_width_count = j;
8851 
8852  if (j < 1) {
8853  pfree(bin_width);
8854  bin_width = NULL;
8855  }
8856  }
8857 
8858  /* right */
8859  if (!PG_ARGISNULL(7))
8860  right = PG_GETARG_BOOL(7);
8861 
8862  /* connect to database */
8863  spi_result = SPI_connect();
8864  if (spi_result != SPI_OK_CONNECT) {
8865 
8866  if (bin_width_count) pfree(bin_width);
8867 
8868  MemoryContextSwitchTo(oldcontext);
8869  elog(ERROR, "RASTER_histogramCoverage: Could not connect to database using SPI");
8870  SRF_RETURN_DONE(funcctx);
8871  }
8872 
8873  /* coverage stats */
8874  len = sizeof(char) * (strlen("SELECT min, max FROM _st_summarystats('','',,::boolean,)") + strlen(tablename) + strlen(colname) + (MAX_INT_CHARLEN * 2) + MAX_DBL_CHARLEN + 1);
8875  sql = (char *) palloc(len);
8876  if (NULL == sql) {
8877 
8878  if (SPI_tuptable) SPI_freetuptable(tuptable);
8879  SPI_finish();
8880 
8881  if (bin_width_count) pfree(bin_width);
8882 
8883  MemoryContextSwitchTo(oldcontext);
8884  elog(ERROR, "RASTER_histogramCoverage: Could not allocate memory for sql");
8885  SRF_RETURN_DONE(funcctx);
8886  }
8887 
8888  /* get stats */
8889  snprintf(sql, len, "SELECT min, max FROM _st_summarystats('%s','%s',%d,%d::boolean,%f)", tablename, colname, bandindex, (exclude_nodata_value ? 1 : 0), sample);
8890  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: %s", sql);
8891  spi_result = SPI_execute(sql, TRUE, 0);
8892  pfree(sql);
8893  if (spi_result != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
8894 
8895  if (SPI_tuptable) SPI_freetuptable(tuptable);
8896  SPI_finish();
8897 
8898  if (bin_width_count) pfree(bin_width);
8899 
8900  MemoryContextSwitchTo(oldcontext);
8901  elog(ERROR, "RASTER_histogramCoverage: Could not get summary stats of coverage");
8902  SRF_RETURN_DONE(funcctx);
8903  }
8904 
8905  tupdesc = SPI_tuptable->tupdesc;
8906  tuptable = SPI_tuptable;
8907  tuple = tuptable->vals[0];
8908 
8909  tmp = SPI_getvalue(tuple, tupdesc, 1);
8910  if (NULL == tmp || !strlen(tmp)) {
8911 
8912  if (SPI_tuptable) SPI_freetuptable(tuptable);
8913  SPI_finish();
8914 
8915  if (bin_width_count) pfree(bin_width);
8916 
8917  MemoryContextSwitchTo(oldcontext);
8918  elog(ERROR, "RASTER_histogramCoverage: Could not get summary stats of coverage");
8919  SRF_RETURN_DONE(funcctx);
8920  }
8921  min = strtod(tmp, NULL);
8922  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: min = %f", min);
8923  pfree(tmp);
8924 
8925  tmp = SPI_getvalue(tuple, tupdesc, 2);
8926  if (NULL == tmp || !strlen(tmp)) {
8927 
8928  if (SPI_tuptable) SPI_freetuptable(tuptable);
8929  SPI_finish();
8930 
8931  if (bin_width_count) pfree(bin_width);
8932 
8933  MemoryContextSwitchTo(oldcontext);
8934  elog(ERROR, "RASTER_histogramCoverage: Could not get summary stats of coverage");
8935  SRF_RETURN_DONE(funcctx);
8936  }
8937  max = strtod(tmp, NULL);
8938  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: max = %f", max);
8939  pfree(tmp);
8940 
8941  /* iterate through rasters of coverage */
8942  /* create sql */
8943  len = sizeof(char) * (strlen("SELECT \"\" FROM \"\" WHERE \"\" IS NOT NULL") + (strlen(colname) * 2) + strlen(tablename) + 1);
8944  sql = (char *) palloc(len);
8945  if (NULL == sql) {
8946 
8947  if (SPI_tuptable) SPI_freetuptable(tuptable);
8948  SPI_finish();
8949 
8950  if (bin_width_count) pfree(bin_width);
8951 
8952  MemoryContextSwitchTo(oldcontext);
8953  elog(ERROR, "RASTER_histogramCoverage: Could not allocate memory for sql");
8954  SRF_RETURN_DONE(funcctx);
8955  }
8956 
8957  /* get cursor */
8958  snprintf(sql, len, "SELECT \"%s\" FROM \"%s\" WHERE \"%s\" IS NOT NULL", colname, tablename, colname);
8959  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: %s", sql);
8960  portal = SPI_cursor_open_with_args(
8961  "coverage",
8962  sql,
8963  0, NULL,
8964  NULL, NULL,
8965  TRUE, 0
8966  );
8967  pfree(sql);
8968 
8969  /* process resultset */
8970  SPI_cursor_fetch(portal, TRUE, 1);
8971  while (SPI_processed == 1 && SPI_tuptable != NULL) {
8972  tupdesc = SPI_tuptable->tupdesc;
8973  tuptable = SPI_tuptable;
8974  tuple = tuptable->vals[0];
8975 
8976  datum = SPI_getbinval(tuple, tupdesc, 1, &isNull);
8977  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
8978 
8979  if (SPI_tuptable) SPI_freetuptable(tuptable);
8980  SPI_cursor_close(portal);
8981  SPI_finish();
8982 
8983  if (NULL != covhist) pfree(covhist);
8984  if (bin_width_count) pfree(bin_width);
8985 
8986  MemoryContextSwitchTo(oldcontext);
8987  elog(ERROR, "RASTER_histogramCoverage: Could not get raster of coverage");
8988  SRF_RETURN_DONE(funcctx);
8989  }
8990  else if (isNull) {
8991  SPI_cursor_fetch(portal, TRUE, 1);
8992  continue;
8993  }
8994 
8995  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(datum);
8996 
8997  raster = rt_raster_deserialize(pgraster, FALSE);
8998  if (!raster) {
8999 
9000  if (SPI_tuptable) SPI_freetuptable(tuptable);
9001  SPI_cursor_close(portal);
9002  SPI_finish();
9003 
9004  if (NULL != covhist) pfree(covhist);
9005  if (bin_width_count) pfree(bin_width);
9006 
9007  MemoryContextSwitchTo(oldcontext);
9008  elog(ERROR, "RASTER_histogramCoverage: Could not deserialize raster");
9009  SRF_RETURN_DONE(funcctx);
9010  }
9011 
9012  /* inspect number of bands*/
9013  num_bands = rt_raster_get_num_bands(raster);
9014  if (bandindex < 1 || bandindex > num_bands) {
9015  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
9016 
9017  rt_raster_destroy(raster);
9018 
9019  if (SPI_tuptable) SPI_freetuptable(tuptable);
9020  SPI_cursor_close(portal);
9021  SPI_finish();
9022 
9023  if (NULL != covhist) pfree(covhist);
9024  if (bin_width_count) pfree(bin_width);
9025 
9026  MemoryContextSwitchTo(oldcontext);
9027  SRF_RETURN_DONE(funcctx);
9028  }
9029 
9030  /* get band */
9031  band = rt_raster_get_band(raster, bandindex - 1);
9032  if (!band) {
9033  elog(NOTICE, "Could not find band at index %d. Returning NULL", bandindex);
9034 
9035  rt_raster_destroy(raster);
9036 
9037  if (SPI_tuptable) SPI_freetuptable(tuptable);
9038  SPI_cursor_close(portal);
9039  SPI_finish();
9040 
9041  if (NULL != covhist) pfree(covhist);
9042  if (bin_width_count) pfree(bin_width);
9043 
9044  MemoryContextSwitchTo(oldcontext);
9045  SRF_RETURN_DONE(funcctx);
9046  }
9047 
9048  /* we need the raw values, hence the non-zero parameter */
9049  stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 1, NULL, NULL, NULL);
9050 
9051  rt_band_destroy(band);
9052  rt_raster_destroy(raster);
9053 
9054  if (NULL == stats) {
9055  elog(NOTICE, "Could not compute summary statistics for band at index %d. Returning NULL", bandindex);
9056 
9057  if (SPI_tuptable) SPI_freetuptable(tuptable);
9058  SPI_cursor_close(portal);
9059  SPI_finish();
9060 
9061  if (NULL != covhist) pfree(covhist);
9062  if (bin_width_count) pfree(bin_width);
9063 
9064  MemoryContextSwitchTo(oldcontext);
9065  SRF_RETURN_DONE(funcctx);
9066  }
9067 
9068  /* get histogram */
9069  if (stats->count > 0) {
9070  hist = rt_band_get_histogram(stats, bin_count, bin_width, bin_width_count, right, min, max, &count);
9071  pfree(stats);
9072  if (NULL == hist || !count) {
9073  elog(NOTICE, "Could not compute histogram for band at index %d", bandindex);
9074 
9075  if (SPI_tuptable) SPI_freetuptable(tuptable);
9076  SPI_cursor_close(portal);
9077  SPI_finish();
9078 
9079  if (NULL != covhist) pfree(covhist);
9080  if (bin_width_count) pfree(bin_width);
9081 
9082  MemoryContextSwitchTo(oldcontext);
9083  SRF_RETURN_DONE(funcctx);
9084  }
9085 
9086  POSTGIS_RT_DEBUGF(3, "%d bins returned", count);
9087 
9088  /* coverage histogram */
9089  if (NULL == covhist) {
9090  covhist = (rt_histogram) SPI_palloc(sizeof(struct rt_histogram_t) * count);
9091  if (NULL == covhist) {
9092 
9093  pfree(hist);
9094  if (SPI_tuptable) SPI_freetuptable(tuptable);
9095  SPI_cursor_close(portal);
9096  SPI_finish();
9097 
9098  if (bin_width_count) pfree(bin_width);
9099 
9100  MemoryContextSwitchTo(oldcontext);
9101  elog(ERROR, "RASTER_histogramCoverage: Could not allocate memory for histogram of coverage");
9102  SRF_RETURN_DONE(funcctx);
9103  }
9104 
9105  for (i = 0; i < count; i++) {
9106  sum += hist[i].count;
9107  covhist[i].count = hist[i].count;
9108  covhist[i].percent = 0;
9109  covhist[i].min = hist[i].min;
9110  covhist[i].max = hist[i].max;
9111  }
9112  }
9113  else {
9114  for (i = 0; i < count; i++) {
9115  sum += hist[i].count;
9116  covhist[i].count += hist[i].count;
9117  }
9118  }
9119 
9120  pfree(hist);
9121 
9122  /* assuming bin_count wasn't set, force consistency */
9123  if (bin_count <= 0) bin_count = count;
9124  }
9125 
9126  /* next record */
9127  SPI_cursor_fetch(portal, TRUE, 1);
9128  }
9129 
9130  if (SPI_tuptable) SPI_freetuptable(tuptable);
9131  SPI_cursor_close(portal);
9132  SPI_finish();
9133 
9134  if (bin_width_count) pfree(bin_width);
9135 
9136  /* finish percent of histogram */
9137  if (sum > 0) {
9138  for (i = 0; i < count; i++)
9139  covhist[i].percent = covhist[i].count / (double) sum;
9140  }
9141 
9142  /* Store needed information */
9143  funcctx->user_fctx = covhist;
9144 
9145  /* total number of tuples to be returned */
9146  funcctx->max_calls = count;
9147 
9148  /* Build a tuple descriptor for our result type */
9149  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
9150  ereport(ERROR, (
9151  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
9152  errmsg(
9153  "function returning record called in context "
9154  "that cannot accept type record"
9155  )
9156  ));
9157  }
9158 
9159  BlessTupleDesc(tupdesc);
9160  funcctx->tuple_desc = tupdesc;
9161 
9162  MemoryContextSwitchTo(oldcontext);
9163  }
9164 
9165  /* stuff done on every call of the function */
9166  funcctx = SRF_PERCALL_SETUP();
9167 
9168  call_cntr = funcctx->call_cntr;
9169  max_calls = funcctx->max_calls;
9170  tupdesc = funcctx->tuple_desc;
9171  covhist2 = funcctx->user_fctx;
9172 
9173  /* do when there is more left to send */
9174  if (call_cntr < max_calls) {
9175  int values_length = 4;
9176  Datum values[values_length];
9177  bool nulls[values_length];
9178  HeapTuple tuple;
9179  Datum result;
9180 
9181  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
9182 
9183  memset(nulls, FALSE, sizeof(bool) * values_length);
9184 
9185  values[0] = Float8GetDatum(covhist2[call_cntr].min);
9186  values[1] = Float8GetDatum(covhist2[call_cntr].max);
9187  values[2] = Int64GetDatum(covhist2[call_cntr].count);
9188  values[3] = Float8GetDatum(covhist2[call_cntr].percent);
9189 
9190  /* build a tuple */
9191  tuple = heap_form_tuple(tupdesc, values, nulls);
9192 
9193  /* make the tuple into a datum */
9194  result = HeapTupleGetDatum(tuple);
9195 
9196  SRF_RETURN_NEXT(funcctx, result);
9197  }
9198  /* do when there is no more left */
9199  else {
9200  pfree(covhist2);
9201  SRF_RETURN_DONE(funcctx);
9202  }
9203 }
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
#define MAX_DBL_CHARLEN
Definition: rt_pg.c:64
double max
Definition: rt_api.h:2295
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 percent
Definition: rt_api.h:2292
#define POSTGIS_RT_DEBUG(level, msg)
Definition: rt_pg.h:58
int count
Definition: genraster.py:57
struct rt_histogram_t * rt_histogram
Definition: rt_api.h:138
#define FLT_EQ(x, y)
Definition: rt_api.h:2159
double min
Definition: rt_api.h:2294
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 MAX_INT_CHARLEN
Definition: rt_pg.c:65
#define TRUE
Definition: dbfopen.c:170
uint32_t count
Definition: rt_api.h:2291

Here is the call graph for this function: