PostGIS  3.0.6dev-r@@SVN_REVISION@@

◆ RASTER_quantileCoverage()

Datum RASTER_quantileCoverage ( PG_FUNCTION_ARGS  )

Definition at line 1973 of file rtpg_statistics.c.

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

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, 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, text_to_cstring(), TRUE, rt_quantile_t::value, genraster::value, and VALUES_LENGTH.

Here is the call graph for this function: