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

Definition at line 10101 of file rt_pg.c.

References ovdump::band, genraster::count, rt_valuecount_t::count, FALSE, FLT_EQ, rt_valuecount_t::percent, POSTGIS_RT_DEBUG, POSTGIS_RT_DEBUGF, rtrowdump::raster, result, rt_band_destroy(), rt_band_get_value_count(), rt_raster_deserialize(), rt_raster_destroy(), rt_raster_get_band(), rt_raster_get_num_bands(), rtgdalraster::sql, TRUE, genraster::value, and rt_valuecount_t::value.

10101  {
10102  FuncCallContext *funcctx;
10103  TupleDesc tupdesc;
10104 
10105  int i;
10106  uint64_t covcount = 0;
10107  uint64_t covtotal = 0;
10108  rt_valuecount covvcnts = NULL;
10109  rt_valuecount covvcnts2;
10110  int call_cntr;
10111  int max_calls;
10112 
10113  POSTGIS_RT_DEBUG(3, "RASTER_valueCountCoverage: Starting");
10114 
10115  /* first call of function */
10116  if (SRF_IS_FIRSTCALL()) {
10117  MemoryContext oldcontext;
10118 
10119  text *tablenametext = NULL;
10120  char *tablename = NULL;
10121  text *colnametext = NULL;
10122  char *colname = NULL;
10123  int32_t bandindex = 1;
10124  bool exclude_nodata_value = TRUE;
10125  double *search_values = NULL;
10126  uint32_t search_values_count = 0;
10127  double roundto = 0;
10128 
10129  int len = 0;
10130  char *sql = NULL;
10131  int spi_result;
10132  Portal portal;
10133  SPITupleTable *tuptable = NULL;
10134  HeapTuple tuple;
10135  Datum datum;
10136  bool isNull = FALSE;
10137  rt_pgraster *pgraster = NULL;
10138  rt_raster raster = NULL;
10139  rt_band band = NULL;
10140  int num_bands = 0;
10141  uint32_t count;
10142  uint32_t total;
10143  rt_valuecount vcnts = NULL;
10144  int exists = 0;
10145 
10146  int j;
10147  int n;
10148 
10149  ArrayType *array;
10150  Oid etype;
10151  Datum *e;
10152  bool *nulls;
10153  int16 typlen;
10154  bool typbyval;
10155  char typalign;
10156 
10157  /* create a function context for cross-call persistence */
10158  funcctx = SRF_FIRSTCALL_INIT();
10159 
10160  /* switch to memory context appropriate for multiple function calls */
10161  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
10162 
10163  /* tablename is null, return null */
10164  if (PG_ARGISNULL(0)) {
10165  elog(NOTICE, "Table name must be provided");
10166  MemoryContextSwitchTo(oldcontext);
10167  SRF_RETURN_DONE(funcctx);
10168  }
10169  tablenametext = PG_GETARG_TEXT_P(0);
10170  tablename = text_to_cstring(tablenametext);
10171  if (!strlen(tablename)) {
10172  elog(NOTICE, "Table name must be provided");
10173  MemoryContextSwitchTo(oldcontext);
10174  SRF_RETURN_DONE(funcctx);
10175  }
10176  POSTGIS_RT_DEBUGF(3, "tablename = %s", tablename);
10177 
10178  /* column name is null, return null */
10179  if (PG_ARGISNULL(1)) {
10180  elog(NOTICE, "Column name must be provided");
10181  MemoryContextSwitchTo(oldcontext);
10182  SRF_RETURN_DONE(funcctx);
10183  }
10184  colnametext = PG_GETARG_TEXT_P(1);
10185  colname = text_to_cstring(colnametext);
10186  if (!strlen(colname)) {
10187  elog(NOTICE, "Column name must be provided");
10188  MemoryContextSwitchTo(oldcontext);
10189  SRF_RETURN_DONE(funcctx);
10190  }
10191  POSTGIS_RT_DEBUGF(3, "colname = %s", colname);
10192 
10193  /* band index is 1-based */
10194  if (!PG_ARGISNULL(2))
10195  bandindex = PG_GETARG_INT32(2);
10196 
10197  /* exclude_nodata_value flag */
10198  if (!PG_ARGISNULL(3))
10199  exclude_nodata_value = PG_GETARG_BOOL(3);
10200 
10201  /* search values */
10202  if (!PG_ARGISNULL(4)) {
10203  array = PG_GETARG_ARRAYTYPE_P(4);
10204  etype = ARR_ELEMTYPE(array);
10205  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
10206 
10207  switch (etype) {
10208  case FLOAT4OID:
10209  case FLOAT8OID:
10210  break;
10211  default:
10212  MemoryContextSwitchTo(oldcontext);
10213  elog(ERROR, "RASTER_valueCountCoverage: Invalid data type for values");
10214  SRF_RETURN_DONE(funcctx);
10215  break;
10216  }
10217 
10218  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
10219  &nulls, &n);
10220 
10221  search_values = palloc(sizeof(double) * n);
10222  for (i = 0, j = 0; i < n; i++) {
10223  if (nulls[i]) continue;
10224 
10225  switch (etype) {
10226  case FLOAT4OID:
10227  search_values[j] = (double) DatumGetFloat4(e[i]);
10228  break;
10229  case FLOAT8OID:
10230  search_values[j] = (double) DatumGetFloat8(e[i]);
10231  break;
10232  }
10233 
10234  POSTGIS_RT_DEBUGF(5, "search_values[%d] = %f", j, search_values[j]);
10235  j++;
10236  }
10237  search_values_count = j;
10238 
10239  if (j < 1) {
10240  pfree(search_values);
10241  search_values = NULL;
10242  }
10243  }
10244 
10245  /* roundto */
10246  if (!PG_ARGISNULL(5)) {
10247  roundto = PG_GETARG_FLOAT8(5);
10248  if (roundto < 0.) roundto = 0;
10249  }
10250 
10251  /* iterate through rasters of coverage */
10252  /* connect to database */
10253  spi_result = SPI_connect();
10254  if (spi_result != SPI_OK_CONNECT) {
10255 
10256  if (search_values_count) pfree(search_values);
10257 
10258  MemoryContextSwitchTo(oldcontext);
10259  elog(ERROR, "RASTER_valueCountCoverage: Could not connect to database using SPI");
10260  SRF_RETURN_DONE(funcctx);
10261  }
10262 
10263  /* create sql */
10264  len = sizeof(char) * (strlen("SELECT \"\" FROM \"\" WHERE \"\" IS NOT NULL") + (strlen(colname) * 2) + strlen(tablename) + 1);
10265  sql = (char *) palloc(len);
10266  if (NULL == sql) {
10267 
10268  if (SPI_tuptable) SPI_freetuptable(tuptable);
10269  SPI_finish();
10270 
10271  if (search_values_count) pfree(search_values);
10272 
10273  MemoryContextSwitchTo(oldcontext);
10274  elog(ERROR, "RASTER_valueCountCoverage: Could not allocate memory for sql");
10275  SRF_RETURN_DONE(funcctx);
10276  }
10277 
10278  /* get cursor */
10279  snprintf(sql, len, "SELECT \"%s\" FROM \"%s\" WHERE \"%s\" IS NOT NULL", colname, tablename, colname);
10280  POSTGIS_RT_DEBUGF(3, "RASTER_valueCountCoverage: %s", sql);
10281  portal = SPI_cursor_open_with_args(
10282  "coverage",
10283  sql,
10284  0, NULL,
10285  NULL, NULL,
10286  TRUE, 0
10287  );
10288  pfree(sql);
10289 
10290  /* process resultset */
10291  SPI_cursor_fetch(portal, TRUE, 1);
10292  while (SPI_processed == 1 && SPI_tuptable != NULL) {
10293  tupdesc = SPI_tuptable->tupdesc;
10294  tuptable = SPI_tuptable;
10295  tuple = tuptable->vals[0];
10296 
10297  datum = SPI_getbinval(tuple, tupdesc, 1, &isNull);
10298  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
10299 
10300  if (SPI_tuptable) SPI_freetuptable(tuptable);
10301  SPI_cursor_close(portal);
10302  SPI_finish();
10303 
10304  if (NULL != covvcnts) pfree(covvcnts);
10305  if (search_values_count) pfree(search_values);
10306 
10307  MemoryContextSwitchTo(oldcontext);
10308  elog(ERROR, "RASTER_valueCountCoverage: Could not get raster of coverage");
10309  SRF_RETURN_DONE(funcctx);
10310  }
10311  else if (isNull) {
10312  SPI_cursor_fetch(portal, TRUE, 1);
10313  continue;
10314  }
10315 
10316  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(datum);
10317 
10318  raster = rt_raster_deserialize(pgraster, FALSE);
10319  if (!raster) {
10320 
10321  if (SPI_tuptable) SPI_freetuptable(tuptable);
10322  SPI_cursor_close(portal);
10323  SPI_finish();
10324 
10325  if (NULL != covvcnts) pfree(covvcnts);
10326  if (search_values_count) pfree(search_values);
10327 
10328  MemoryContextSwitchTo(oldcontext);
10329  elog(ERROR, "RASTER_valueCountCoverage: Could not deserialize raster");
10330  SRF_RETURN_DONE(funcctx);
10331  }
10332 
10333  /* inspect number of bands*/
10334  num_bands = rt_raster_get_num_bands(raster);
10335  if (bandindex < 1 || bandindex > num_bands) {
10336  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
10337 
10338  rt_raster_destroy(raster);
10339 
10340  if (SPI_tuptable) SPI_freetuptable(tuptable);
10341  SPI_cursor_close(portal);
10342  SPI_finish();
10343 
10344  if (NULL != covvcnts) pfree(covvcnts);
10345  if (search_values_count) pfree(search_values);
10346 
10347  MemoryContextSwitchTo(oldcontext);
10348  SRF_RETURN_DONE(funcctx);
10349  }
10350 
10351  /* get band */
10352  band = rt_raster_get_band(raster, bandindex - 1);
10353  if (!band) {
10354  elog(NOTICE, "Could not find band at index %d. Returning NULL", bandindex);
10355 
10356  rt_raster_destroy(raster);
10357 
10358  if (SPI_tuptable) SPI_freetuptable(tuptable);
10359  SPI_cursor_close(portal);
10360  SPI_finish();
10361 
10362  if (NULL != covvcnts) pfree(covvcnts);
10363  if (search_values_count) pfree(search_values);
10364 
10365  MemoryContextSwitchTo(oldcontext);
10366  SRF_RETURN_DONE(funcctx);
10367  }
10368 
10369  /* get counts of values */
10370  vcnts = rt_band_get_value_count(band, (int) exclude_nodata_value, search_values, search_values_count, roundto, &total, &count);
10371  rt_band_destroy(band);
10372  rt_raster_destroy(raster);
10373  if (NULL == vcnts || !count) {
10374  elog(NOTICE, "Could not count the values for band at index %d", bandindex);
10375 
10376  if (SPI_tuptable) SPI_freetuptable(tuptable);
10377  SPI_cursor_close(portal);
10378  SPI_finish();
10379 
10380  if (NULL != covvcnts) free(covvcnts);
10381  if (search_values_count) pfree(search_values);
10382 
10383  MemoryContextSwitchTo(oldcontext);
10384  SRF_RETURN_DONE(funcctx);
10385  }
10386 
10387  POSTGIS_RT_DEBUGF(3, "%d value counts returned", count);
10388 
10389  if (NULL == covvcnts) {
10390  covvcnts = (rt_valuecount) SPI_palloc(sizeof(struct rt_valuecount_t) * count);
10391  if (NULL == covvcnts) {
10392 
10393  if (SPI_tuptable) SPI_freetuptable(tuptable);
10394  SPI_cursor_close(portal);
10395  SPI_finish();
10396 
10397  if (search_values_count) pfree(search_values);
10398 
10399  MemoryContextSwitchTo(oldcontext);
10400  elog(ERROR, "RASTER_valueCountCoverage: Could not allocate memory for value counts of coverage");
10401  SRF_RETURN_DONE(funcctx);
10402  }
10403 
10404  for (i = 0; i < count; i++) {
10405  covvcnts[i].value = vcnts[i].value;
10406  covvcnts[i].count = vcnts[i].count;
10407  covvcnts[i].percent = -1;
10408  }
10409 
10410  covcount = count;
10411  }
10412  else {
10413  for (i = 0; i < count; i++) {
10414  exists = 0;
10415 
10416  for (j = 0; j < covcount; j++) {
10417  if (FLT_EQ(vcnts[i].value, covvcnts[j].value)) {
10418  exists = 1;
10419  break;
10420  }
10421  }
10422 
10423  if (exists) {
10424  covvcnts[j].count += vcnts[i].count;
10425  }
10426  else {
10427  covcount++;
10428  covvcnts = SPI_repalloc(covvcnts, sizeof(struct rt_valuecount_t) * covcount);
10429  if (NULL == covvcnts) {
10430 
10431  if (SPI_tuptable) SPI_freetuptable(tuptable);
10432  SPI_cursor_close(portal);
10433  SPI_finish();
10434 
10435  if (search_values_count) pfree(search_values);
10436  if (NULL != covvcnts) free(covvcnts);
10437 
10438  MemoryContextSwitchTo(oldcontext);
10439  elog(ERROR, "RASTER_valueCountCoverage: Could not change allocated memory for value counts of coverage");
10440  SRF_RETURN_DONE(funcctx);
10441  }
10442 
10443  covvcnts[covcount - 1].value = vcnts[i].value;
10444  covvcnts[covcount - 1].count = vcnts[i].count;
10445  covvcnts[covcount - 1].percent = -1;
10446  }
10447  }
10448  }
10449 
10450  covtotal += total;
10451 
10452  pfree(vcnts);
10453 
10454  /* next record */
10455  SPI_cursor_fetch(portal, TRUE, 1);
10456  }
10457 
10458  if (SPI_tuptable) SPI_freetuptable(tuptable);
10459  SPI_cursor_close(portal);
10460  SPI_finish();
10461 
10462  if (search_values_count) pfree(search_values);
10463 
10464  /* compute percentages */
10465  for (i = 0; i < covcount; i++) {
10466  covvcnts[i].percent = (double) covvcnts[i].count / covtotal;
10467  }
10468 
10469  /* Store needed information */
10470  funcctx->user_fctx = covvcnts;
10471 
10472  /* total number of tuples to be returned */
10473  funcctx->max_calls = covcount;
10474 
10475  /* Build a tuple descriptor for our result type */
10476  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
10477  ereport(ERROR, (
10478  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
10479  errmsg(
10480  "function returning record called in context "
10481  "that cannot accept type record"
10482  )
10483  ));
10484  }
10485 
10486  BlessTupleDesc(tupdesc);
10487  funcctx->tuple_desc = tupdesc;
10488 
10489  MemoryContextSwitchTo(oldcontext);
10490  }
10491 
10492  /* stuff done on every call of the function */
10493  funcctx = SRF_PERCALL_SETUP();
10494 
10495  call_cntr = funcctx->call_cntr;
10496  max_calls = funcctx->max_calls;
10497  tupdesc = funcctx->tuple_desc;
10498  covvcnts2 = funcctx->user_fctx;
10499 
10500  /* do when there is more left to send */
10501  if (call_cntr < max_calls) {
10502  int values_length = 3;
10503  Datum values[values_length];
10504  bool nulls[values_length];
10505  HeapTuple tuple;
10506  Datum result;
10507 
10508  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
10509 
10510  memset(nulls, FALSE, sizeof(bool) * values_length);
10511 
10512  values[0] = Float8GetDatum(covvcnts2[call_cntr].value);
10513  values[1] = UInt32GetDatum(covvcnts2[call_cntr].count);
10514  values[2] = Float8GetDatum(covvcnts2[call_cntr].percent);
10515 
10516  /* build a tuple */
10517  tuple = heap_form_tuple(tupdesc, values, nulls);
10518 
10519  /* make the tuple into a datum */
10520  result = HeapTupleGetDatum(tuple);
10521 
10522  SRF_RETURN_NEXT(funcctx, result);
10523  }
10524  /* do when there is no more left */
10525  else {
10526  pfree(covvcnts2);
10527  SRF_RETURN_DONE(funcctx);
10528  }
10529 }
struct rt_valuecount_t * rt_valuecount
Definition: rt_api.h:140
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:2342
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
double value
Definition: rt_api.h:2341
char ** result
Definition: liblwgeom.h:218
#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
double percent
Definition: rt_api.h:2343
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
rt_valuecount rt_band_get_value_count(rt_band band, int exclude_nodata_value, double *search_values, uint32_t search_values_count, double roundto, uint32_t *rtn_total, uint32_t *rtn_count)
Count the number of times provided value(s) occur in the band.
Definition: rt_api.c:4750

Here is the call graph for this function: