PostGIS  2.2.7dev-r@@SVN_REVISION@@
Datum RASTER_quantile ( PG_FUNCTION_ARGS  )

Definition at line 1725 of file rtpg_statistics.c.

References ovdump::band, genraster::count, rt_bandstats_t::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, and rt_bandstats_t::values.

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