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

◆ RASTER_valueCountCoverage()

Datum RASTER_valueCountCoverage ( PG_FUNCTION_ARGS  )

Definition at line 1656 of file rtpg_statistics.c.

1656 {
1657 FuncCallContext *funcctx;
1658 TupleDesc tupdesc;
1659
1660 uint32_t i;
1661 uint64_t covcount = 0;
1662 uint64_t covtotal = 0;
1663 rt_valuecount covvcnts = NULL;
1664 rt_valuecount covvcnts2;
1665 int call_cntr;
1666 int max_calls;
1667
1668 POSTGIS_RT_DEBUG(3, "RASTER_valueCountCoverage: Starting");
1669
1670 /* first call of function */
1671 if (SRF_IS_FIRSTCALL()) {
1672 MemoryContext oldcontext;
1673
1674 text *tablenametext = NULL;
1675 char *tablename = NULL;
1676 text *colnametext = NULL;
1677 char *colname = NULL;
1678 int32_t bandindex = 1;
1679 bool exclude_nodata_value = TRUE;
1680 double *search_values = NULL;
1681 uint32_t search_values_count = 0;
1682 double roundto = 0;
1683
1684 int len = 0;
1685 char *sql = NULL;
1686 int spi_result;
1687 Portal portal;
1688 HeapTuple tuple;
1689 Datum datum;
1690 bool isNull = FALSE;
1691 rt_pgraster *pgraster = NULL;
1692 rt_raster raster = NULL;
1693 rt_band band = NULL;
1694 int num_bands = 0;
1695 uint32_t count;
1696 uint32_t total;
1697 rt_valuecount vcnts = NULL;
1698 int exists = 0;
1699
1700 uint32_t j;
1701 int n;
1702
1703 ArrayType *array;
1704 Oid etype;
1705 Datum *e;
1706 bool *nulls;
1707 int16 typlen;
1708 bool typbyval;
1709 char typalign;
1710
1711 /* create a function context for cross-call persistence */
1712 funcctx = SRF_FIRSTCALL_INIT();
1713
1714 /* switch to memory context appropriate for multiple function calls */
1715 oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
1716
1717 /* tablename is null, return null */
1718 if (PG_ARGISNULL(0)) {
1719 elog(NOTICE, "Table name must be provided");
1720 MemoryContextSwitchTo(oldcontext);
1721 SRF_RETURN_DONE(funcctx);
1722 }
1723 tablenametext = PG_GETARG_TEXT_P(0);
1724 tablename = text_to_cstring(tablenametext);
1725 if (!strlen(tablename)) {
1726 elog(NOTICE, "Table name must be provided");
1727 MemoryContextSwitchTo(oldcontext);
1728 SRF_RETURN_DONE(funcctx);
1729 }
1730 POSTGIS_RT_DEBUGF(3, "tablename = %s", tablename);
1731
1732 /* column name is null, return null */
1733 if (PG_ARGISNULL(1)) {
1734 elog(NOTICE, "Column name must be provided");
1735 MemoryContextSwitchTo(oldcontext);
1736 SRF_RETURN_DONE(funcctx);
1737 }
1738 colnametext = PG_GETARG_TEXT_P(1);
1739 colname = text_to_cstring(colnametext);
1740 if (!strlen(colname)) {
1741 elog(NOTICE, "Column name must be provided");
1742 MemoryContextSwitchTo(oldcontext);
1743 SRF_RETURN_DONE(funcctx);
1744 }
1745 POSTGIS_RT_DEBUGF(3, "colname = %s", colname);
1746
1747 /* band index is 1-based */
1748 if (!PG_ARGISNULL(2))
1749 bandindex = PG_GETARG_INT32(2);
1750
1751 /* exclude_nodata_value flag */
1752 if (!PG_ARGISNULL(3))
1753 exclude_nodata_value = PG_GETARG_BOOL(3);
1754
1755 /* search values */
1756 if (!PG_ARGISNULL(4)) {
1757 array = PG_GETARG_ARRAYTYPE_P(4);
1758 etype = ARR_ELEMTYPE(array);
1759 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1760
1761 switch (etype) {
1762 case FLOAT4OID:
1763 case FLOAT8OID:
1764 break;
1765 default:
1766 MemoryContextSwitchTo(oldcontext);
1767 elog(ERROR, "RASTER_valueCountCoverage: Invalid data type for values");
1768 SRF_RETURN_DONE(funcctx);
1769 break;
1770 }
1771
1772 deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1773 &nulls, &n);
1774
1775 search_values = palloc(sizeof(double) * n);
1776 for (i = 0, j = 0; i < (uint32_t) n; i++) {
1777 if (nulls[i]) continue;
1778
1779 switch (etype) {
1780 case FLOAT4OID:
1781 search_values[j] = (double) DatumGetFloat4(e[i]);
1782 break;
1783 case FLOAT8OID:
1784 search_values[j] = (double) DatumGetFloat8(e[i]);
1785 break;
1786 }
1787
1788 POSTGIS_RT_DEBUGF(5, "search_values[%d] = %f", j, search_values[j]);
1789 j++;
1790 }
1791 search_values_count = j;
1792
1793 if (j < 1) {
1794 pfree(search_values);
1795 search_values = NULL;
1796 }
1797 }
1798
1799 /* roundto */
1800 if (!PG_ARGISNULL(5)) {
1801 roundto = PG_GETARG_FLOAT8(5);
1802 if (roundto < 0.) roundto = 0;
1803 }
1804
1805 /* iterate through rasters of coverage */
1806 /* connect to database */
1807 spi_result = SPI_connect();
1808 if (spi_result != SPI_OK_CONNECT) {
1809
1810 if (search_values_count) pfree(search_values);
1811
1812 MemoryContextSwitchTo(oldcontext);
1813 elog(ERROR, "RASTER_valueCountCoverage: Cannot connect to database using SPI");
1814 SRF_RETURN_DONE(funcctx);
1815 }
1816
1817 /* create sql */
1818 len = sizeof(char) * (strlen("SELECT \"\" FROM \"\" WHERE \"\" IS NOT NULL") + (strlen(colname) * 2) + strlen(tablename) + 1);
1819 sql = (char *) palloc(len);
1820 if (NULL == sql) {
1821
1822 if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
1823 SPI_finish();
1824
1825 if (search_values_count) pfree(search_values);
1826
1827 MemoryContextSwitchTo(oldcontext);
1828 elog(ERROR, "RASTER_valueCountCoverage: Cannot allocate memory for sql");
1829 SRF_RETURN_DONE(funcctx);
1830 }
1831
1832 /* get cursor */
1833 snprintf(sql, len, "SELECT \"%s\" FROM \"%s\" WHERE \"%s\" IS NOT NULL", colname, tablename, colname);
1834 POSTGIS_RT_DEBUGF(3, "RASTER_valueCountCoverage: %s", sql);
1835 portal = SPI_cursor_open_with_args(
1836 "coverage",
1837 sql,
1838 0, NULL,
1839 NULL, NULL,
1840 TRUE, 0
1841 );
1842 pfree(sql);
1843
1844 /* process resultset */
1845 SPI_cursor_fetch(portal, TRUE, 1);
1846 while (SPI_processed == 1 && SPI_tuptable != NULL) {
1847 tupdesc = SPI_tuptable->tupdesc;
1848 tuple = SPI_tuptable->vals[0];
1849
1850 datum = SPI_getbinval(tuple, tupdesc, 1, &isNull);
1851 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1852
1853 SPI_freetuptable(SPI_tuptable);
1854 SPI_cursor_close(portal);
1855 SPI_finish();
1856
1857 if (NULL != covvcnts) pfree(covvcnts);
1858 if (search_values_count) pfree(search_values);
1859
1860 MemoryContextSwitchTo(oldcontext);
1861 elog(ERROR, "RASTER_valueCountCoverage: Cannot get raster of coverage");
1862 SRF_RETURN_DONE(funcctx);
1863 }
1864 else if (isNull) {
1865 SPI_cursor_fetch(portal, TRUE, 1);
1866 continue;
1867 }
1868
1869 pgraster = (rt_pgraster *) PG_DETOAST_DATUM(datum);
1870
1871 raster = rt_raster_deserialize(pgraster, FALSE);
1872 if (!raster) {
1873
1874 SPI_freetuptable(SPI_tuptable);
1875 SPI_cursor_close(portal);
1876 SPI_finish();
1877
1878 if (NULL != covvcnts) pfree(covvcnts);
1879 if (search_values_count) pfree(search_values);
1880
1881 MemoryContextSwitchTo(oldcontext);
1882 elog(ERROR, "RASTER_valueCountCoverage: Cannot deserialize raster");
1883 SRF_RETURN_DONE(funcctx);
1884 }
1885
1886 /* inspect number of bands*/
1887 num_bands = rt_raster_get_num_bands(raster);
1888 if (bandindex < 1 || bandindex > num_bands) {
1889 elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
1890
1891 rt_raster_destroy(raster);
1892
1893 SPI_freetuptable(SPI_tuptable);
1894 SPI_cursor_close(portal);
1895 SPI_finish();
1896
1897 if (NULL != covvcnts) pfree(covvcnts);
1898 if (search_values_count) pfree(search_values);
1899
1900 MemoryContextSwitchTo(oldcontext);
1901 SRF_RETURN_DONE(funcctx);
1902 }
1903
1904 /* get band */
1905 band = rt_raster_get_band(raster, bandindex - 1);
1906 if (!band) {
1907 elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
1908
1909 rt_raster_destroy(raster);
1910
1911 SPI_freetuptable(SPI_tuptable);
1912 SPI_cursor_close(portal);
1913 SPI_finish();
1914
1915 if (NULL != covvcnts) pfree(covvcnts);
1916 if (search_values_count) pfree(search_values);
1917
1918 MemoryContextSwitchTo(oldcontext);
1919 SRF_RETURN_DONE(funcctx);
1920 }
1921
1922 /* get counts of values */
1923 vcnts = rt_band_get_value_count(band, (int) exclude_nodata_value, search_values, search_values_count, roundto, &total, &count);
1924 rt_band_destroy(band);
1925 rt_raster_destroy(raster);
1926 if (NULL == vcnts || !count) {
1927 elog(NOTICE, "Cannot count the values for band at index %d", bandindex);
1928
1929 SPI_freetuptable(SPI_tuptable);
1930 SPI_cursor_close(portal);
1931 SPI_finish();
1932
1933 if (NULL != covvcnts) free(covvcnts);
1934 if (search_values_count) pfree(search_values);
1935
1936 MemoryContextSwitchTo(oldcontext);
1937 SRF_RETURN_DONE(funcctx);
1938 }
1939
1940 POSTGIS_RT_DEBUGF(3, "%d value counts returned", count);
1941
1942 if (NULL == covvcnts) {
1943 covvcnts = (rt_valuecount) SPI_palloc(sizeof(struct rt_valuecount_t) * count);
1944 if (NULL == covvcnts) {
1945
1946 SPI_freetuptable(SPI_tuptable);
1947 SPI_cursor_close(portal);
1948 SPI_finish();
1949
1950 if (search_values_count) pfree(search_values);
1951
1952 MemoryContextSwitchTo(oldcontext);
1953 elog(ERROR, "RASTER_valueCountCoverage: Cannot allocate memory for value counts of coverage");
1954 SRF_RETURN_DONE(funcctx);
1955 }
1956
1957 for (i = 0; i < count; i++) {
1958 covvcnts[i].value = vcnts[i].value;
1959 covvcnts[i].count = vcnts[i].count;
1960 covvcnts[i].percent = -1;
1961 }
1962
1963 covcount = count;
1964 }
1965 else {
1966 for (i = 0; i < count; i++) {
1967 exists = 0;
1968
1969 for (j = 0; j < covcount; j++) {
1970 if (FLT_EQ(vcnts[i].value, covvcnts[j].value)) {
1971 exists = 1;
1972 break;
1973 }
1974 }
1975
1976 if (exists) {
1977 covvcnts[j].count += vcnts[i].count;
1978 }
1979 else {
1980 covcount++;
1981 covvcnts = SPI_repalloc(covvcnts, sizeof(struct rt_valuecount_t) * covcount);
1982 if (!covvcnts) {
1983 SPI_freetuptable(SPI_tuptable);
1984 SPI_cursor_close(portal);
1985 SPI_finish();
1986
1987 if (search_values_count) pfree(search_values);
1988
1989 MemoryContextSwitchTo(oldcontext);
1990 elog(ERROR, "RASTER_valueCountCoverage: Cannot change allocated memory for value counts of coverage");
1991 SRF_RETURN_DONE(funcctx);
1992 }
1993
1994 covvcnts[covcount - 1].value = vcnts[i].value;
1995 covvcnts[covcount - 1].count = vcnts[i].count;
1996 covvcnts[covcount - 1].percent = -1;
1997 }
1998 }
1999 }
2000
2001 covtotal += total;
2002
2003 pfree(vcnts);
2004
2005 /* next record */
2006 SPI_cursor_fetch(portal, TRUE, 1);
2007 }
2008
2009 if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
2010 SPI_cursor_close(portal);
2011 SPI_finish();
2012
2013 if (search_values_count) pfree(search_values);
2014
2015 /* compute percentages */
2016 for (i = 0; i < covcount; i++) {
2017 covvcnts[i].percent = (double) covvcnts[i].count / covtotal;
2018 }
2019
2020 /* Store needed information */
2021 funcctx->user_fctx = covvcnts;
2022
2023 /* total number of tuples to be returned */
2024 funcctx->max_calls = covcount;
2025
2026 /* Build a tuple descriptor for our result type */
2027 if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
2028 ereport(ERROR, (
2029 errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
2030 errmsg(
2031 "function returning record called in context "
2032 "that cannot accept type record"
2033 )
2034 ));
2035 }
2036
2037 BlessTupleDesc(tupdesc);
2038 funcctx->tuple_desc = tupdesc;
2039
2040 MemoryContextSwitchTo(oldcontext);
2041 }
2042
2043 /* stuff done on every call of the function */
2044 funcctx = SRF_PERCALL_SETUP();
2045
2046 call_cntr = funcctx->call_cntr;
2047 max_calls = funcctx->max_calls;
2048 tupdesc = funcctx->tuple_desc;
2049 covvcnts2 = funcctx->user_fctx;
2050
2051 /* do when there is more left to send */
2052 if (call_cntr < max_calls) {
2053 Datum values[VALUES_LENGTH];
2054 bool nulls[VALUES_LENGTH];
2055 HeapTuple tuple;
2056 Datum result;
2057
2058 POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
2059
2060 memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
2061
2062 values[0] = Float8GetDatum(covvcnts2[call_cntr].value);
2063 values[1] = UInt32GetDatum(covvcnts2[call_cntr].count);
2064 values[2] = Float8GetDatum(covvcnts2[call_cntr].percent);
2065
2066 /* build a tuple */
2067 tuple = heap_form_tuple(tupdesc, values, nulls);
2068
2069 /* make the tuple into a datum */
2070 result = HeapTupleGetDatum(tuple);
2071
2072 SRF_RETURN_NEXT(funcctx, result);
2073 }
2074 /* do when there is no more left */
2075 else {
2076 if (covvcnts2) pfree(covvcnts2);
2077 SRF_RETURN_DONE(funcctx);
2078 }
2079}
result
char result[OUT_DOUBLE_BUFFER_SIZE]
Definition cu_print.c:267
TRUE
#define TRUE
Definition dbfopen.c:73
FALSE
#define FALSE
Definition dbfopen.c:72
rt_raster_destroy
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition rt_raster.c:86
FLT_EQ
#define FLT_EQ(x, y)
Definition librtcore.h:2436
rt_band_destroy
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition rt_band.c:499
rt_raster_get_num_bands
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition rt_raster.c:376
rt_valuecount
struct rt_valuecount_t * rt_valuecount
Definition librtcore.h:154
rt_band_get_value_count
rt_valuecount rt_band_get_value_count(rt_band band, int exclude_nodata_value, double *search_values, uint32_t search_values_count, double roundto, uint32_t *rtn_total, uint32_t *rtn_count)
Count the number of times provided value(s) occur in the band.
Definition rt_statistics.c:1609
rt_raster_deserialize
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition rt_serialize.c:732
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:385
free
void free(void *)
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:125
VALUES_LENGTH
#define VALUES_LENGTH
Definition rtpg_statistics.c:65
POSTGIS_RT_DEBUG
#define POSTGIS_RT_DEBUG(level, msg)
Definition rtpostgis.h:65
POSTGIS_RT_DEBUGF
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition rtpostgis.h:69
rt_band_t
Definition librtcore.h:2515
rt_raster_serialized_t
Struct definitions.
Definition librtcore.h:2452
rt_raster_t
Definition librtcore.h:2486
rt_valuecount_t::count
uint32_t count
Definition librtcore.h:2627
rt_valuecount_t::percent
double percent
Definition librtcore.h:2628
rt_valuecount_t::value
double value
Definition librtcore.h:2626
rt_valuecount_t
Definition librtcore.h:2625

References rt_valuecount_t::count, FALSE, FLT_EQ, free(), rt_valuecount_t::percent, POSTGIS_RT_DEBUG, POSTGIS_RT_DEBUGF, result, rt_band_destroy(), rt_band_get_value_count(), rt_raster_deserialize(), rt_raster_destroy(), rt_raster_get_band(), rt_raster_get_num_bands(), TRUE, rt_valuecount_t::value, and VALUES_LENGTH.

Here is the call graph for this function: