PostGIS  2.5.0beta2dev-r@@SVN_REVISION@@

◆ RASTER_quantile()

Datum RASTER_quantile ( PG_FUNCTION_ARGS  )

Definition at line 1717 of file rtpg_statistics.c.

References ovdump::band, genraster::count, rt_bandstats_t::count, rtpg_summarystats_arg_t::exclude_nodata_value, FALSE, FLT_EQ, PG_FUNCTION_INFO_V1(), POSTGIS_RT_DEBUGF, rtrowdump::raster, RASTER_quantileCoverage(), 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(), rtpg_summarystats_arg_t::sample, rtpg_summarystats_arg_t::stats, TRUE, genraster::value, and rt_bandstats_t::values.

Referenced by RASTER_histogramCoverage().

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