PostGIS  2.5.7dev-r@@SVN_REVISION@@

◆ RASTER_quantileCoverage()

Datum RASTER_quantileCoverage ( PG_FUNCTION_ARGS  )

Definition at line 1969 of file rtpg_statistics.c.

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

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: