PostGIS 3.0.6dev-r@@SVN_REVISION@@
Loading...
Searching...
No Matches

◆ RASTER_quantile()

Datum RASTER_quantile ( PG_FUNCTION_ARGS  )

Definition at line 1724 of file rtpg_statistics.c.

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

References rt_bandstats_t::count, FALSE, FLT_EQ, POSTGIS_RT_DEBUGF, 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, rt_bandstats_t::values, and VALUES_LENGTH.

Here is the call graph for this function: