PostGIS  2.4.9dev-r@@SVN_REVISION@@

◆ RASTER_quantile()

Datum RASTER_quantile ( PG_FUNCTION_ARGS  )

Definition at line 1728 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, rt_bandstats_t::values, and VALUES_LENGTH.

Referenced by RASTER_histogramCoverage().

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