PostGIS  2.5.0beta2dev-r@@SVN_REVISION@@

◆ RASTER_quantileCoverage()

Datum RASTER_quantileCoverage ( PG_FUNCTION_ARGS  )

Definition at line 1967 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, text_to_cstring(), TRUE, genraster::value, and rt_quantile_t::value.

Referenced by RASTER_quantile().

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