PostGIS  2.4.9dev-r@@SVN_REVISION@@

◆ RASTER_quantileCoverage()

Datum RASTER_quantileCoverage ( PG_FUNCTION_ARGS  )

Definition at line 1977 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, POSTGIS_RT_DEBUG, POSTGIS_RT_DEBUGF, rt_quantile_t::quantile, quantile_llist::quantile, quantile_llist_destroy(), rtrowdump::raster, 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, rt_quantile_t::value, and VALUES_LENGTH.

Referenced by RASTER_quantile().

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