PostGIS  2.3.8dev-r@@SVN_REVISION@@

◆ RASTER_quantileCoverage()

Datum RASTER_quantileCoverage ( PG_FUNCTION_ARGS  )

Definition at line 1975 of file rtpg_statistics.c.

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

Referenced by RASTER_quantile().

1976 {
1977  FuncCallContext *funcctx;
1978  TupleDesc tupdesc;
1979 
1980  int i;
1981  rt_quantile covquant = NULL;
1982  rt_quantile covquant2;
1983  int call_cntr;
1984  int max_calls;
1985 
1986  POSTGIS_RT_DEBUG(3, "RASTER_quantileCoverage: Starting");
1987 
1988  /* first call of function */
1989  if (SRF_IS_FIRSTCALL()) {
1990  MemoryContext oldcontext;
1991 
1992  text *tablenametext = NULL;
1993  char *tablename = NULL;
1994  text *colnametext = NULL;
1995  char *colname = NULL;
1996  int32_t bandindex = 1;
1997  bool exclude_nodata_value = TRUE;
1998  double sample = 0;
1999  double *quantiles = NULL;
2000  uint32_t quantiles_count = 0;
2001  double quantile = 0;
2002  uint32_t count;
2003 
2004  int len = 0;
2005  char *sql = NULL;
2006  char *tmp = NULL;
2007  uint64_t cov_count = 0;
2008  int spi_result;
2009  Portal portal;
2010  SPITupleTable *tuptable = NULL;
2011  HeapTuple tuple;
2012  Datum datum;
2013  bool isNull = FALSE;
2014 
2015  rt_pgraster *pgraster = NULL;
2016  rt_raster raster = NULL;
2017  rt_band band = NULL;
2018  int num_bands = 0;
2019  struct quantile_llist *qlls = NULL;
2020  uint32_t qlls_count;
2021 
2022  int j;
2023  int n;
2024 
2025  ArrayType *array;
2026  Oid etype;
2027  Datum *e;
2028  bool *nulls;
2029  int16 typlen;
2030  bool typbyval;
2031  char typalign;
2032 
2033  POSTGIS_RT_DEBUG(3, "RASTER_quantileCoverage: first call of function");
2034 
2035  /* create a function context for cross-call persistence */
2036  funcctx = SRF_FIRSTCALL_INIT();
2037 
2038  /* switch to memory context appropriate for multiple function calls */
2039  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
2040 
2041  /* tablename is null, return null */
2042  if (PG_ARGISNULL(0)) {
2043  elog(NOTICE, "Table name must be provided");
2044  MemoryContextSwitchTo(oldcontext);
2045  SRF_RETURN_DONE(funcctx);
2046  }
2047  tablenametext = PG_GETARG_TEXT_P(0);
2048  tablename = text_to_cstring(tablenametext);
2049  if (!strlen(tablename)) {
2050  elog(NOTICE, "Table name must be provided");
2051  MemoryContextSwitchTo(oldcontext);
2052  SRF_RETURN_DONE(funcctx);
2053  }
2054  POSTGIS_RT_DEBUGF(3, "RASTER_quantileCoverage: tablename = %s", tablename);
2055 
2056  /* column name is null, return null */
2057  if (PG_ARGISNULL(1)) {
2058  elog(NOTICE, "Column name must be provided");
2059  MemoryContextSwitchTo(oldcontext);
2060  SRF_RETURN_DONE(funcctx);
2061  }
2062  colnametext = PG_GETARG_TEXT_P(1);
2063  colname = text_to_cstring(colnametext);
2064  if (!strlen(colname)) {
2065  elog(NOTICE, "Column name must be provided");
2066  MemoryContextSwitchTo(oldcontext);
2067  SRF_RETURN_DONE(funcctx);
2068  }
2069  POSTGIS_RT_DEBUGF(3, "RASTER_quantileCoverage: colname = %s", colname);
2070 
2071  /* band index is 1-based */
2072  if (!PG_ARGISNULL(2))
2073  bandindex = PG_GETARG_INT32(2);
2074 
2075  /* exclude_nodata_value flag */
2076  if (!PG_ARGISNULL(3))
2077  exclude_nodata_value = PG_GETARG_BOOL(3);
2078 
2079  /* sample % */
2080  if (!PG_ARGISNULL(4)) {
2081  sample = PG_GETARG_FLOAT8(4);
2082  if (sample < 0 || sample > 1) {
2083  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
2084  MemoryContextSwitchTo(oldcontext);
2085  SRF_RETURN_DONE(funcctx);
2086  }
2087  else if (FLT_EQ(sample, 0.0))
2088  sample = 1;
2089  }
2090  else
2091  sample = 1;
2092 
2093  /* quantiles */
2094  if (!PG_ARGISNULL(5)) {
2095  array = PG_GETARG_ARRAYTYPE_P(5);
2096  etype = ARR_ELEMTYPE(array);
2097  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
2098 
2099  switch (etype) {
2100  case FLOAT4OID:
2101  case FLOAT8OID:
2102  break;
2103  default:
2104  MemoryContextSwitchTo(oldcontext);
2105  elog(ERROR, "RASTER_quantileCoverage: Invalid data type for quantiles");
2106  SRF_RETURN_DONE(funcctx);
2107  break;
2108  }
2109 
2110  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
2111  &nulls, &n);
2112 
2113  quantiles = palloc(sizeof(double) * n);
2114  for (i = 0, j = 0; i < n; i++) {
2115  if (nulls[i]) continue;
2116 
2117  switch (etype) {
2118  case FLOAT4OID:
2119  quantile = (double) DatumGetFloat4(e[i]);
2120  break;
2121  case FLOAT8OID:
2122  quantile = (double) DatumGetFloat8(e[i]);
2123  break;
2124  }
2125 
2126  if (quantile < 0 || quantile > 1) {
2127  elog(NOTICE, "Invalid value for quantile (must be between 0 and 1). Returning NULL");
2128  pfree(quantiles);
2129  MemoryContextSwitchTo(oldcontext);
2130  SRF_RETURN_DONE(funcctx);
2131  }
2132 
2133  quantiles[j] = quantile;
2134  POSTGIS_RT_DEBUGF(5, "quantiles[%d] = %f", j, quantiles[j]);
2135  j++;
2136  }
2137  quantiles_count = j;
2138 
2139  if (j < 1) {
2140  pfree(quantiles);
2141  quantiles = NULL;
2142  }
2143  }
2144 
2145  /* coverage stats */
2146  /* connect to database */
2147  spi_result = SPI_connect();
2148  if (spi_result != SPI_OK_CONNECT) {
2149  MemoryContextSwitchTo(oldcontext);
2150  elog(ERROR, "RASTER_quantileCoverage: Cannot connect to database using SPI");
2151  SRF_RETURN_DONE(funcctx);
2152  }
2153 
2154  len = sizeof(char) * (strlen("SELECT count FROM _st_summarystats('','',,::boolean,)") + strlen(tablename) + strlen(colname) + (MAX_INT_CHARLEN * 2) + MAX_DBL_CHARLEN + 1);
2155  sql = (char *) palloc(len);
2156  if (NULL == sql) {
2157 
2158  if (SPI_tuptable) SPI_freetuptable(tuptable);
2159  SPI_finish();
2160 
2161  MemoryContextSwitchTo(oldcontext);
2162  elog(ERROR, "RASTER_quantileCoverage: Cannot allocate memory for sql");
2163  SRF_RETURN_DONE(funcctx);
2164  }
2165 
2166  /* get stats */
2167  snprintf(sql, len, "SELECT count FROM _st_summarystats('%s','%s',%d,%d::boolean,%f)", tablename, colname, bandindex, (exclude_nodata_value ? 1 : 0), sample);
2168  POSTGIS_RT_DEBUGF(3, "stats sql: %s", sql);
2169  spi_result = SPI_execute(sql, TRUE, 0);
2170  pfree(sql);
2171  if (spi_result != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
2172 
2173  if (SPI_tuptable) SPI_freetuptable(tuptable);
2174  SPI_finish();
2175 
2176  MemoryContextSwitchTo(oldcontext);
2177  elog(ERROR, "RASTER_quantileCoverage: Cannot get summary stats of coverage");
2178  SRF_RETURN_DONE(funcctx);
2179  }
2180 
2181  tupdesc = SPI_tuptable->tupdesc;
2182  tuptable = SPI_tuptable;
2183  tuple = tuptable->vals[0];
2184 
2185  tmp = SPI_getvalue(tuple, tupdesc, 1);
2186  if (NULL == tmp || !strlen(tmp)) {
2187 
2188  if (SPI_tuptable) SPI_freetuptable(tuptable);
2189  SPI_finish();
2190 
2191  MemoryContextSwitchTo(oldcontext);
2192  elog(ERROR, "RASTER_quantileCoverage: Cannot get summary stats of coverage");
2193  SRF_RETURN_DONE(funcctx);
2194  }
2195  cov_count = strtol(tmp, NULL, 10);
2196  POSTGIS_RT_DEBUGF(3, "covcount = %d", (int) cov_count);
2197  pfree(tmp);
2198 
2199  /* iterate through rasters of coverage */
2200  /* create sql */
2201  len = sizeof(char) * (strlen("SELECT \"\" FROM \"\" WHERE \"\" IS NOT NULL") + (strlen(colname) * 2) + strlen(tablename) + 1);
2202  sql = (char *) palloc(len);
2203  if (NULL == sql) {
2204 
2205  if (SPI_tuptable) SPI_freetuptable(tuptable);
2206  SPI_finish();
2207 
2208  MemoryContextSwitchTo(oldcontext);
2209  elog(ERROR, "RASTER_quantileCoverage: Cannot allocate memory for sql");
2210  SRF_RETURN_DONE(funcctx);
2211  }
2212 
2213  /* get cursor */
2214  snprintf(sql, len, "SELECT \"%s\" FROM \"%s\" WHERE \"%s\" IS NOT NULL", colname, tablename, colname);
2215  POSTGIS_RT_DEBUGF(3, "coverage sql: %s", sql);
2216  portal = SPI_cursor_open_with_args(
2217  "coverage",
2218  sql,
2219  0, NULL,
2220  NULL, NULL,
2221  TRUE, 0
2222  );
2223  pfree(sql);
2224 
2225  /* process resultset */
2226  SPI_cursor_fetch(portal, TRUE, 1);
2227  while (SPI_processed == 1 && SPI_tuptable != NULL) {
2228  if (NULL != covquant) pfree(covquant);
2229 
2230  tupdesc = SPI_tuptable->tupdesc;
2231  tuptable = SPI_tuptable;
2232  tuple = tuptable->vals[0];
2233 
2234  datum = SPI_getbinval(tuple, tupdesc, 1, &isNull);
2235  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
2236 
2237  if (SPI_tuptable) SPI_freetuptable(tuptable);
2238  SPI_cursor_close(portal);
2239  SPI_finish();
2240 
2241  MemoryContextSwitchTo(oldcontext);
2242  elog(ERROR, "RASTER_quantileCoverage: Cannot get raster of coverage");
2243  SRF_RETURN_DONE(funcctx);
2244  }
2245  else if (isNull) {
2246  SPI_cursor_fetch(portal, TRUE, 1);
2247  continue;
2248  }
2249 
2250  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(datum);
2251 
2252  raster = rt_raster_deserialize(pgraster, FALSE);
2253  if (!raster) {
2254 
2255  if (SPI_tuptable) SPI_freetuptable(tuptable);
2256  SPI_cursor_close(portal);
2257  SPI_finish();
2258 
2259  MemoryContextSwitchTo(oldcontext);
2260  elog(ERROR, "RASTER_quantileCoverage: Cannot deserialize raster");
2261  SRF_RETURN_DONE(funcctx);
2262  }
2263 
2264  /* inspect number of bands*/
2265  num_bands = rt_raster_get_num_bands(raster);
2266  if (bandindex < 1 || bandindex > num_bands) {
2267  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
2268 
2269  rt_raster_destroy(raster);
2270 
2271  if (SPI_tuptable) SPI_freetuptable(tuptable);
2272  SPI_cursor_close(portal);
2273  SPI_finish();
2274 
2275  MemoryContextSwitchTo(oldcontext);
2276  SRF_RETURN_DONE(funcctx);
2277  }
2278 
2279  /* get band */
2280  band = rt_raster_get_band(raster, bandindex - 1);
2281  if (!band) {
2282  elog(NOTICE, "Cannot find raster band of index %d. Returning NULL", bandindex);
2283 
2284  rt_raster_destroy(raster);
2285 
2286  if (SPI_tuptable) SPI_freetuptable(tuptable);
2287  SPI_cursor_close(portal);
2288  SPI_finish();
2289 
2290  MemoryContextSwitchTo(oldcontext);
2291  SRF_RETURN_DONE(funcctx);
2292  }
2293 
2294  covquant = rt_band_get_quantiles_stream(
2295  band,
2296  exclude_nodata_value, sample, cov_count,
2297  &qlls, &qlls_count,
2298  quantiles, quantiles_count,
2299  &count
2300  );
2301 
2302  rt_band_destroy(band);
2303  rt_raster_destroy(raster);
2304 
2305  if (NULL == covquant || !count) {
2306  elog(NOTICE, "Cannot compute quantiles for band at index %d", bandindex);
2307 
2308  if (SPI_tuptable) SPI_freetuptable(tuptable);
2309  SPI_cursor_close(portal);
2310  SPI_finish();
2311 
2312  MemoryContextSwitchTo(oldcontext);
2313  SRF_RETURN_DONE(funcctx);
2314  }
2315 
2316  /* next record */
2317  SPI_cursor_fetch(portal, TRUE, 1);
2318  }
2319 
2320  covquant2 = SPI_palloc(sizeof(struct rt_quantile_t) * count);
2321  for (i = 0; i < count; i++) {
2322  covquant2[i].quantile = covquant[i].quantile;
2323  covquant2[i].has_value = covquant[i].has_value;
2324  if (covquant2[i].has_value)
2325  covquant2[i].value = covquant[i].value;
2326  }
2327 
2328  if (NULL != covquant) pfree(covquant);
2329  quantile_llist_destroy(&qlls, qlls_count);
2330 
2331  if (SPI_tuptable) SPI_freetuptable(tuptable);
2332  SPI_cursor_close(portal);
2333  SPI_finish();
2334 
2335  if (quantiles_count) pfree(quantiles);
2336 
2337  POSTGIS_RT_DEBUGF(3, "%d quantiles returned", count);
2338 
2339  /* Store needed information */
2340  funcctx->user_fctx = covquant2;
2341 
2342  /* total number of tuples to be returned */
2343  funcctx->max_calls = count;
2344 
2345  /* Build a tuple descriptor for our result type */
2346  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
2347  ereport(ERROR, (
2348  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
2349  errmsg(
2350  "function returning record called in context "
2351  "that cannot accept type record"
2352  )
2353  ));
2354  }
2355 
2356  BlessTupleDesc(tupdesc);
2357  funcctx->tuple_desc = tupdesc;
2358 
2359  MemoryContextSwitchTo(oldcontext);
2360  }
2361 
2362  /* stuff done on every call of the function */
2363  funcctx = SRF_PERCALL_SETUP();
2364 
2365  call_cntr = funcctx->call_cntr;
2366  max_calls = funcctx->max_calls;
2367  tupdesc = funcctx->tuple_desc;
2368  covquant2 = funcctx->user_fctx;
2369 
2370  /* do when there is more left to send */
2371  if (call_cntr < max_calls) {
2372  int values_length = 2;
2373  Datum values[values_length];
2374  bool nulls[values_length];
2375  HeapTuple tuple;
2376  Datum result;
2377 
2378  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
2379 
2380  memset(nulls, FALSE, sizeof(bool) * values_length);
2381 
2382  values[0] = Float8GetDatum(covquant2[call_cntr].quantile);
2383  if (covquant2[call_cntr].has_value)
2384  values[1] = Float8GetDatum(covquant2[call_cntr].value);
2385  else
2386  nulls[1] = TRUE;
2387 
2388  /* build a tuple */
2389  tuple = heap_form_tuple(tupdesc, values, nulls);
2390 
2391  /* make the tuple into a datum */
2392  result = HeapTupleGetDatum(tuple);
2393 
2394  SRF_RETURN_NEXT(funcctx, result);
2395  }
2396  /* do when there is no more left */
2397  else {
2398  POSTGIS_RT_DEBUG(3, "done");
2399  pfree(covquant2);
2400  SRF_RETURN_DONE(funcctx);
2401  }
2402 }
int quantile_llist_destroy(struct quantile_llist **list, uint32_t list_count)
int rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:372
double quantile
Definition: librtcore.h:2357
raster
Be careful!! Zeros function&#39;s input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
band
Definition: ovdump.py:57
double quantile
Definition: librtcore.h:2349
#define FLT_EQ(x, y)
Definition: librtcore.h:2197
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition: rt_band.c:242
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.
double value
Definition: librtcore.h:2350
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition: rtpostgis.h:57
#define MAX_INT_CHARLEN
Definition: rtpostgis.h:68
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:381
int count
Definition: genraster.py:56
#define MAX_DBL_CHARLEN
Definition: rtpostgis.h:67
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:82
#define FALSE
Definition: dbfopen.c:168
Struct definitions.
Definition: librtcore.h:2213
uint32_t has_value
Definition: librtcore.h:2351
int value
Definition: genraster.py:61
#define POSTGIS_RT_DEBUG(level, msg)
Definition: rtpostgis.h:53
#define TRUE
Definition: dbfopen.c:169
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition: rt_serialize.c:717
Here is the call graph for this function:
Here is the caller graph for this function: