PostGIS  3.4.0dev-r@@SVN_REVISION@@

◆ 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  pfree(covvcnts2);
2077  SRF_RETURN_DONE(funcctx);
2078  }
2079 }
result
char result[OUT_DOUBLE_BUFFER_SIZE]
Definition: cu_print.c:262
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:2387
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:376
rt_valuecount
struct rt_valuecount_t * rt_valuecount
Definition: librtcore.h:155
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:1608
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:385
free
void free(void *)
genraster.value
int value
Definition: genraster.py:62
genraster.count
int count
Definition: genraster.py:57
ovdump.band
band
Definition: ovdump.py:58
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:1442
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:2466
rt_raster_serialized_t
Struct definitions.
Definition: librtcore.h:2403
rt_raster_t
Definition: librtcore.h:2437
rt_valuecount_t::count
uint32_t count
Definition: librtcore.h:2578
rt_valuecount_t::percent
double percent
Definition: librtcore.h:2579
rt_valuecount_t::value
double value
Definition: librtcore.h:2577
rt_valuecount_t
Definition: librtcore.h:2576

References ovdump::band, rt_valuecount_t::count, genraster::count, FALSE, FLT_EQ, free(), rt_valuecount_t::percent, POSTGIS_RT_DEBUG, POSTGIS_RT_DEBUGF, rtrowdump::raster, result, rt_band_destroy(), rt_band_get_value_count(), rt_raster_deserialize(), rt_raster_destroy(), rt_raster_get_band(), rt_raster_get_num_bands(), rtgdalraster::sql, TRUE, rt_valuecount_t::value, genraster::value, and VALUES_LENGTH.

Here is the call graph for this function: