PostGIS  2.5.7dev-r@@SVN_REVISION@@

◆ RASTER_quantile()

Datum RASTER_quantile ( PG_FUNCTION_ARGS  )

Definition at line 1720 of file rtpg_statistics.c.

1721 {
1722  FuncCallContext *funcctx;
1723  TupleDesc tupdesc;
1724 
1725  int i;
1726  rt_quantile quant;
1727  rt_quantile quant2;
1728  int call_cntr;
1729  int max_calls;
1730 
1731  /* first call of function */
1732  if (SRF_IS_FIRSTCALL()) {
1733  MemoryContext oldcontext;
1734 
1735  rt_pgraster *pgraster = NULL;
1736  rt_raster raster = NULL;
1737  rt_band band = NULL;
1738  int32_t bandindex = 0;
1739  int num_bands = 0;
1740  bool exclude_nodata_value = TRUE;
1741  double sample = 0;
1742  double *quantiles = NULL;
1743  uint32_t quantiles_count = 0;
1744  double quantile = 0;
1745  rt_bandstats stats = NULL;
1746  uint32_t count;
1747 
1748  int j;
1749  int n;
1750 
1751  ArrayType *array;
1752  Oid etype;
1753  Datum *e;
1754  bool *nulls;
1755  int16 typlen;
1756  bool typbyval;
1757  char typalign;
1758 
1759  /* create a function context for cross-call persistence */
1760  funcctx = SRF_FIRSTCALL_INIT();
1761 
1762  /* switch to memory context appropriate for multiple function calls */
1763  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
1764 
1765  /* pgraster is null, return nothing */
1766  if (PG_ARGISNULL(0)) {
1767  MemoryContextSwitchTo(oldcontext);
1768  SRF_RETURN_DONE(funcctx);
1769  }
1770  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1771 
1772  raster = rt_raster_deserialize(pgraster, FALSE);
1773  if (!raster) {
1774  PG_FREE_IF_COPY(pgraster, 0);
1775  MemoryContextSwitchTo(oldcontext);
1776  elog(ERROR, "RASTER_quantile: Cannot deserialize raster");
1777  SRF_RETURN_DONE(funcctx);
1778  }
1779 
1780  /* band index is 1-based */
1781  bandindex = PG_GETARG_INT32(1);
1782  num_bands = rt_raster_get_num_bands(raster);
1783  if (bandindex < 1 || bandindex > num_bands) {
1784  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
1785  rt_raster_destroy(raster);
1786  PG_FREE_IF_COPY(pgraster, 0);
1787  MemoryContextSwitchTo(oldcontext);
1788  SRF_RETURN_DONE(funcctx);
1789  }
1790 
1791  /* exclude_nodata_value flag */
1792  if (!PG_ARGISNULL(2))
1793  exclude_nodata_value = PG_GETARG_BOOL(2);
1794 
1795  /* sample % */
1796  if (!PG_ARGISNULL(3)) {
1797  sample = PG_GETARG_FLOAT8(3);
1798  if (sample < 0 || sample > 1) {
1799  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
1800  rt_raster_destroy(raster);
1801  PG_FREE_IF_COPY(pgraster, 0);
1802  MemoryContextSwitchTo(oldcontext);
1803  SRF_RETURN_DONE(funcctx);
1804  }
1805  else if (FLT_EQ(sample, 0.0))
1806  sample = 1;
1807  }
1808  else
1809  sample = 1;
1810 
1811  /* quantiles */
1812  if (!PG_ARGISNULL(4)) {
1813  array = PG_GETARG_ARRAYTYPE_P(4);
1814  etype = ARR_ELEMTYPE(array);
1815  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1816 
1817  switch (etype) {
1818  case FLOAT4OID:
1819  case FLOAT8OID:
1820  break;
1821  default:
1822  rt_raster_destroy(raster);
1823  PG_FREE_IF_COPY(pgraster, 0);
1824  MemoryContextSwitchTo(oldcontext);
1825  elog(ERROR, "RASTER_quantile: Invalid data type for quantiles");
1826  SRF_RETURN_DONE(funcctx);
1827  break;
1828  }
1829 
1830  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1831  &nulls, &n);
1832 
1833  quantiles = palloc(sizeof(double) * n);
1834  for (i = 0, j = 0; i < n; i++) {
1835  if (nulls[i]) continue;
1836 
1837  switch (etype) {
1838  case FLOAT4OID:
1839  quantile = (double) DatumGetFloat4(e[i]);
1840  break;
1841  case FLOAT8OID:
1842  quantile = (double) DatumGetFloat8(e[i]);
1843  break;
1844  }
1845 
1846  if (quantile < 0 || quantile > 1) {
1847  elog(NOTICE, "Invalid value for quantile (must be between 0 and 1). Returning NULL");
1848  pfree(quantiles);
1849  rt_raster_destroy(raster);
1850  PG_FREE_IF_COPY(pgraster, 0);
1851  MemoryContextSwitchTo(oldcontext);
1852  SRF_RETURN_DONE(funcctx);
1853  }
1854 
1855  quantiles[j] = quantile;
1856  POSTGIS_RT_DEBUGF(5, "quantiles[%d] = %f", j, quantiles[j]);
1857  j++;
1858  }
1859  quantiles_count = j;
1860 
1861  if (j < 1) {
1862  pfree(quantiles);
1863  quantiles = NULL;
1864  }
1865  }
1866 
1867  /* get band */
1868  band = rt_raster_get_band(raster, bandindex - 1);
1869  if (!band) {
1870  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
1871  rt_raster_destroy(raster);
1872  PG_FREE_IF_COPY(pgraster, 0);
1873  MemoryContextSwitchTo(oldcontext);
1874  SRF_RETURN_DONE(funcctx);
1875  }
1876 
1877  /* get stats */
1878  stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 1, NULL, NULL, NULL);
1879  rt_band_destroy(band);
1880  rt_raster_destroy(raster);
1881  PG_FREE_IF_COPY(pgraster, 0);
1882  if (NULL == stats || NULL == stats->values) {
1883  elog(NOTICE, "Cannot retrieve summary statistics for band at index %d", bandindex);
1884  MemoryContextSwitchTo(oldcontext);
1885  SRF_RETURN_DONE(funcctx);
1886  }
1887  else if (stats->count < 1) {
1888  elog(NOTICE, "Cannot compute quantiles for band at index %d as the band has no values", bandindex);
1889  MemoryContextSwitchTo(oldcontext);
1890  SRF_RETURN_DONE(funcctx);
1891  }
1892 
1893  /* get quantiles */
1894  quant = rt_band_get_quantiles(stats, quantiles, quantiles_count, &count);
1895  if (quantiles_count) pfree(quantiles);
1896  pfree(stats);
1897  if (NULL == quant || !count) {
1898  elog(NOTICE, "Cannot compute quantiles for band at index %d", bandindex);
1899  MemoryContextSwitchTo(oldcontext);
1900  SRF_RETURN_DONE(funcctx);
1901  }
1902 
1903  POSTGIS_RT_DEBUGF(3, "%d quantiles returned", count);
1904 
1905  /* Store needed information */
1906  funcctx->user_fctx = quant;
1907 
1908  /* total number of tuples to be returned */
1909  funcctx->max_calls = count;
1910 
1911  /* Build a tuple descriptor for our result type */
1912  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
1913  ereport(ERROR, (
1914  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1915  errmsg(
1916  "function returning record called in context "
1917  "that cannot accept type record"
1918  )
1919  ));
1920  }
1921 
1922  BlessTupleDesc(tupdesc);
1923  funcctx->tuple_desc = tupdesc;
1924 
1925  MemoryContextSwitchTo(oldcontext);
1926  }
1927 
1928  /* stuff done on every call of the function */
1929  funcctx = SRF_PERCALL_SETUP();
1930 
1931  call_cntr = funcctx->call_cntr;
1932  max_calls = funcctx->max_calls;
1933  tupdesc = funcctx->tuple_desc;
1934  quant2 = funcctx->user_fctx;
1935 
1936  /* do when there is more left to send */
1937  if (call_cntr < max_calls) {
1938  Datum values[VALUES_LENGTH];
1939  bool nulls[VALUES_LENGTH];
1940  HeapTuple tuple;
1941  Datum result;
1942 
1943  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
1944 
1945  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
1946 
1947  values[0] = Float8GetDatum(quant2[call_cntr].quantile);
1948  values[1] = Float8GetDatum(quant2[call_cntr].value);
1949 
1950  /* build a tuple */
1951  tuple = heap_form_tuple(tupdesc, values, nulls);
1952 
1953  /* make the tuple into a datum */
1954  result = HeapTupleGetDatum(tuple);
1955 
1956  SRF_RETURN_NEXT(funcctx, result);
1957  }
1958  /* do when there is no more left */
1959  else {
1960  pfree(quant2);
1961  SRF_RETURN_DONE(funcctx);
1962  }
1963 }
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_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_band_get_quantiles
rt_quantile rt_band_get_quantiles(rt_bandstats stats, double *quantiles, int quantiles_count, uint32_t *rtn_count)
Compute the default set of or requested quantiles for a set of data the quantile formula used is same...
Definition: rt_statistics.c:672
rt_band_get_summary_stats
rt_bandstats rt_band_get_summary_stats(rt_band band, int exclude_nodata_value, double sample, int inc_vals, uint64_t *cK, double *cM, double *cQ)
Compute summary statistics for a band.
Definition: rt_statistics.c:110
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
VALUES_LENGTH
#define VALUES_LENGTH
Definition: rtpg_statistics.c:2396
POSTGIS_RT_DEBUGF
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition: rtpostgis.h:65
rt_band_t
Definition: librtcore.h:2313
rt_bandstats_t::count
uint32_t count
Definition: librtcore.h:2360
rt_bandstats_t::values
double * values
Definition: librtcore.h:2368
rt_bandstats_t
Definition: librtcore.h:2358
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, rt_bandstats_t::count, genraster::count, FALSE, FLT_EQ, POSTGIS_RT_DEBUGF, rtrowdump::raster, rt_band_destroy(), rt_band_get_quantiles(), rt_band_get_summary_stats(), rt_raster_deserialize(), rt_raster_destroy(), rt_raster_get_band(), rt_raster_get_num_bands(), TRUE, genraster::value, rt_bandstats_t::values, and VALUES_LENGTH.

Here is the call graph for this function: