PostGIS  2.4.9dev-r@@SVN_REVISION@@

◆ RASTER_histogramCoverage()

Datum RASTER_histogramCoverage ( PG_FUNCTION_ARGS  )

Definition at line 1195 of file rtpg_statistics.c.

References ovdump::band, genraster::count, rt_bandstats_t::count, rt_histogram_t::count, rtpg_summarystats_arg_t::exclude_nodata_value, FALSE, FLT_EQ, rt_histogram_t::max, MAX_DBL_CHARLEN, MAX_INT_CHARLEN, rt_histogram_t::min, rt_histogram_t::percent, PG_FUNCTION_INFO_V1(), POSTGIS_RT_DEBUG, POSTGIS_RT_DEBUGF, rtrowdump::raster, RASTER_quantile(), rt_band_destroy(), rt_band_get_histogram(), 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, rtgdalraster::sql, rtpg_summarystats_arg_t::stats, TRUE, and VALUES_LENGTH.

Referenced by RASTER_histogram().

1196 {
1197  FuncCallContext *funcctx;
1198  TupleDesc tupdesc;
1199 
1200  int i;
1201  rt_histogram covhist = NULL;
1202  rt_histogram covhist2;
1203  int call_cntr;
1204  int max_calls;
1205 
1206  POSTGIS_RT_DEBUG(3, "RASTER_histogramCoverage: Starting");
1207 
1208  /* first call of function */
1209  if (SRF_IS_FIRSTCALL()) {
1210  MemoryContext oldcontext;
1211 
1212  text *tablenametext = NULL;
1213  char *tablename = NULL;
1214  text *colnametext = NULL;
1215  char *colname = NULL;
1216  int32_t bandindex = 1;
1217  bool exclude_nodata_value = TRUE;
1218  double sample = 0;
1219  uint32_t bin_count = 0;
1220  double *bin_width = NULL;
1221  uint32_t bin_width_count = 0;
1222  double width = 0;
1223  bool right = FALSE;
1224  uint32_t count;
1225 
1226  int len = 0;
1227  char *sql = NULL;
1228  char *tmp = NULL;
1229  double min = 0;
1230  double max = 0;
1231  int spi_result;
1232  Portal portal;
1233  SPITupleTable *tuptable = NULL;
1234  HeapTuple tuple;
1235  Datum datum;
1236  bool isNull = FALSE;
1237 
1238  rt_pgraster *pgraster = NULL;
1239  rt_raster raster = NULL;
1240  rt_band band = NULL;
1241  int num_bands = 0;
1242  rt_bandstats stats = NULL;
1243  rt_histogram hist;
1244  uint64_t sum = 0;
1245 
1246  int j;
1247  int n;
1248 
1249  ArrayType *array;
1250  Oid etype;
1251  Datum *e;
1252  bool *nulls;
1253  int16 typlen;
1254  bool typbyval;
1255  char typalign;
1256 
1257  POSTGIS_RT_DEBUG(3, "RASTER_histogramCoverage: first call of function");
1258 
1259  /* create a function context for cross-call persistence */
1260  funcctx = SRF_FIRSTCALL_INIT();
1261 
1262  /* switch to memory context appropriate for multiple function calls */
1263  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
1264 
1265  /* tablename is null, return null */
1266  if (PG_ARGISNULL(0)) {
1267  elog(NOTICE, "Table name must be provided");
1268  MemoryContextSwitchTo(oldcontext);
1269  SRF_RETURN_DONE(funcctx);
1270  }
1271  tablenametext = PG_GETARG_TEXT_P(0);
1272  tablename = text_to_cstring(tablenametext);
1273  if (!strlen(tablename)) {
1274  elog(NOTICE, "Table name must be provided");
1275  MemoryContextSwitchTo(oldcontext);
1276  SRF_RETURN_DONE(funcctx);
1277  }
1278  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: tablename = %s", tablename);
1279 
1280  /* column name is null, return null */
1281  if (PG_ARGISNULL(1)) {
1282  elog(NOTICE, "Column name must be provided");
1283  MemoryContextSwitchTo(oldcontext);
1284  SRF_RETURN_DONE(funcctx);
1285  }
1286  colnametext = PG_GETARG_TEXT_P(1);
1287  colname = text_to_cstring(colnametext);
1288  if (!strlen(colname)) {
1289  elog(NOTICE, "Column name must be provided");
1290  MemoryContextSwitchTo(oldcontext);
1291  SRF_RETURN_DONE(funcctx);
1292  }
1293  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: colname = %s", colname);
1294 
1295  /* band index is 1-based */
1296  if (!PG_ARGISNULL(2))
1297  bandindex = PG_GETARG_INT32(2);
1298 
1299  /* exclude_nodata_value flag */
1300  if (!PG_ARGISNULL(3))
1301  exclude_nodata_value = PG_GETARG_BOOL(3);
1302 
1303  /* sample % */
1304  if (!PG_ARGISNULL(4)) {
1305  sample = PG_GETARG_FLOAT8(4);
1306  if (sample < 0 || sample > 1) {
1307  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
1308  MemoryContextSwitchTo(oldcontext);
1309  SRF_RETURN_DONE(funcctx);
1310  }
1311  else if (FLT_EQ(sample, 0.0))
1312  sample = 1;
1313  }
1314  else
1315  sample = 1;
1316 
1317  /* bin_count */
1318  if (!PG_ARGISNULL(5)) {
1319  bin_count = PG_GETARG_INT32(5);
1320  if (bin_count < 1) bin_count = 0;
1321  }
1322 
1323  /* bin_width */
1324  if (!PG_ARGISNULL(6)) {
1325  array = PG_GETARG_ARRAYTYPE_P(6);
1326  etype = ARR_ELEMTYPE(array);
1327  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1328 
1329  switch (etype) {
1330  case FLOAT4OID:
1331  case FLOAT8OID:
1332  break;
1333  default:
1334  MemoryContextSwitchTo(oldcontext);
1335  elog(ERROR, "RASTER_histogramCoverage: Invalid data type for width");
1336  SRF_RETURN_DONE(funcctx);
1337  break;
1338  }
1339 
1340  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1341  &nulls, &n);
1342 
1343  bin_width = palloc(sizeof(double) * n);
1344  for (i = 0, j = 0; i < n; i++) {
1345  if (nulls[i]) continue;
1346 
1347  switch (etype) {
1348  case FLOAT4OID:
1349  width = (double) DatumGetFloat4(e[i]);
1350  break;
1351  case FLOAT8OID:
1352  width = (double) DatumGetFloat8(e[i]);
1353  break;
1354  }
1355 
1356  if (width < 0 || FLT_EQ(width, 0.0)) {
1357  elog(NOTICE, "Invalid value for width (must be greater than 0). Returning NULL");
1358  pfree(bin_width);
1359  MemoryContextSwitchTo(oldcontext);
1360  SRF_RETURN_DONE(funcctx);
1361  }
1362 
1363  bin_width[j] = width;
1364  POSTGIS_RT_DEBUGF(5, "bin_width[%d] = %f", j, bin_width[j]);
1365  j++;
1366  }
1367  bin_width_count = j;
1368 
1369  if (j < 1) {
1370  pfree(bin_width);
1371  bin_width = NULL;
1372  }
1373  }
1374 
1375  /* right */
1376  if (!PG_ARGISNULL(7))
1377  right = PG_GETARG_BOOL(7);
1378 
1379  /* connect to database */
1380  spi_result = SPI_connect();
1381  if (spi_result != SPI_OK_CONNECT) {
1382 
1383  if (bin_width_count) pfree(bin_width);
1384 
1385  MemoryContextSwitchTo(oldcontext);
1386  elog(ERROR, "RASTER_histogramCoverage: Cannot connect to database using SPI");
1387  SRF_RETURN_DONE(funcctx);
1388  }
1389 
1390  /* coverage stats */
1391  len = sizeof(char) * (strlen("SELECT min, max FROM _st_summarystats('','',,::boolean,)") + strlen(tablename) + strlen(colname) + (MAX_INT_CHARLEN * 2) + MAX_DBL_CHARLEN + 1);
1392  sql = (char *) palloc(len);
1393  if (NULL == sql) {
1394 
1395  if (SPI_tuptable) SPI_freetuptable(tuptable);
1396  SPI_finish();
1397 
1398  if (bin_width_count) pfree(bin_width);
1399 
1400  MemoryContextSwitchTo(oldcontext);
1401  elog(ERROR, "RASTER_histogramCoverage: Cannot allocate memory for sql");
1402  SRF_RETURN_DONE(funcctx);
1403  }
1404 
1405  /* get stats */
1406  snprintf(sql, len, "SELECT min, max FROM _st_summarystats('%s','%s',%d,%d::boolean,%f)", tablename, colname, bandindex, (exclude_nodata_value ? 1 : 0), sample);
1407  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: %s", sql);
1408  spi_result = SPI_execute(sql, TRUE, 0);
1409  pfree(sql);
1410  if (spi_result != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
1411 
1412  if (SPI_tuptable) SPI_freetuptable(tuptable);
1413  SPI_finish();
1414 
1415  if (bin_width_count) pfree(bin_width);
1416 
1417  MemoryContextSwitchTo(oldcontext);
1418  elog(ERROR, "RASTER_histogramCoverage: Cannot get summary stats of coverage");
1419  SRF_RETURN_DONE(funcctx);
1420  }
1421 
1422  tupdesc = SPI_tuptable->tupdesc;
1423  tuptable = SPI_tuptable;
1424  tuple = tuptable->vals[0];
1425 
1426  tmp = SPI_getvalue(tuple, tupdesc, 1);
1427  if (NULL == tmp || !strlen(tmp)) {
1428 
1429  if (SPI_tuptable) SPI_freetuptable(tuptable);
1430  SPI_finish();
1431 
1432  if (bin_width_count) pfree(bin_width);
1433 
1434  MemoryContextSwitchTo(oldcontext);
1435  elog(ERROR, "RASTER_histogramCoverage: Cannot get summary stats of coverage");
1436  SRF_RETURN_DONE(funcctx);
1437  }
1438  min = strtod(tmp, NULL);
1439  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: min = %f", min);
1440  pfree(tmp);
1441 
1442  tmp = SPI_getvalue(tuple, tupdesc, 2);
1443  if (NULL == tmp || !strlen(tmp)) {
1444 
1445  if (SPI_tuptable) SPI_freetuptable(tuptable);
1446  SPI_finish();
1447 
1448  if (bin_width_count) pfree(bin_width);
1449 
1450  MemoryContextSwitchTo(oldcontext);
1451  elog(ERROR, "RASTER_histogramCoverage: Cannot get summary stats of coverage");
1452  SRF_RETURN_DONE(funcctx);
1453  }
1454  max = strtod(tmp, NULL);
1455  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: max = %f", max);
1456  pfree(tmp);
1457 
1458  /* iterate through rasters of coverage */
1459  /* create sql */
1460  len = sizeof(char) * (strlen("SELECT \"\" FROM \"\" WHERE \"\" IS NOT NULL") + (strlen(colname) * 2) + strlen(tablename) + 1);
1461  sql = (char *) palloc(len);
1462  if (NULL == sql) {
1463 
1464  if (SPI_tuptable) SPI_freetuptable(tuptable);
1465  SPI_finish();
1466 
1467  if (bin_width_count) pfree(bin_width);
1468 
1469  MemoryContextSwitchTo(oldcontext);
1470  elog(ERROR, "RASTER_histogramCoverage: Cannot allocate memory for sql");
1471  SRF_RETURN_DONE(funcctx);
1472  }
1473 
1474  /* get cursor */
1475  snprintf(sql, len, "SELECT \"%s\" FROM \"%s\" WHERE \"%s\" IS NOT NULL", colname, tablename, colname);
1476  POSTGIS_RT_DEBUGF(3, "RASTER_histogramCoverage: %s", sql);
1477  portal = SPI_cursor_open_with_args(
1478  "coverage",
1479  sql,
1480  0, NULL,
1481  NULL, NULL,
1482  TRUE, 0
1483  );
1484  pfree(sql);
1485 
1486  /* process resultset */
1487  SPI_cursor_fetch(portal, TRUE, 1);
1488  while (SPI_processed == 1 && SPI_tuptable != NULL) {
1489  tupdesc = SPI_tuptable->tupdesc;
1490  tuptable = SPI_tuptable;
1491  tuple = tuptable->vals[0];
1492 
1493  datum = SPI_getbinval(tuple, tupdesc, 1, &isNull);
1494  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1495 
1496  if (SPI_tuptable) SPI_freetuptable(tuptable);
1497  SPI_cursor_close(portal);
1498  SPI_finish();
1499 
1500  if (NULL != covhist) pfree(covhist);
1501  if (bin_width_count) pfree(bin_width);
1502 
1503  MemoryContextSwitchTo(oldcontext);
1504  elog(ERROR, "RASTER_histogramCoverage: Cannot get raster of coverage");
1505  SRF_RETURN_DONE(funcctx);
1506  }
1507  else if (isNull) {
1508  SPI_cursor_fetch(portal, TRUE, 1);
1509  continue;
1510  }
1511 
1512  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(datum);
1513 
1514  raster = rt_raster_deserialize(pgraster, FALSE);
1515  if (!raster) {
1516 
1517  if (SPI_tuptable) SPI_freetuptable(tuptable);
1518  SPI_cursor_close(portal);
1519  SPI_finish();
1520 
1521  if (NULL != covhist) pfree(covhist);
1522  if (bin_width_count) pfree(bin_width);
1523 
1524  MemoryContextSwitchTo(oldcontext);
1525  elog(ERROR, "RASTER_histogramCoverage: Cannot deserialize raster");
1526  SRF_RETURN_DONE(funcctx);
1527  }
1528 
1529  /* inspect number of bands*/
1530  num_bands = rt_raster_get_num_bands(raster);
1531  if (bandindex < 1 || bandindex > num_bands) {
1532  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
1533 
1534  rt_raster_destroy(raster);
1535 
1536  if (SPI_tuptable) SPI_freetuptable(tuptable);
1537  SPI_cursor_close(portal);
1538  SPI_finish();
1539 
1540  if (NULL != covhist) pfree(covhist);
1541  if (bin_width_count) pfree(bin_width);
1542 
1543  MemoryContextSwitchTo(oldcontext);
1544  SRF_RETURN_DONE(funcctx);
1545  }
1546 
1547  /* get band */
1548  band = rt_raster_get_band(raster, bandindex - 1);
1549  if (!band) {
1550  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
1551 
1552  rt_raster_destroy(raster);
1553 
1554  if (SPI_tuptable) SPI_freetuptable(tuptable);
1555  SPI_cursor_close(portal);
1556  SPI_finish();
1557 
1558  if (NULL != covhist) pfree(covhist);
1559  if (bin_width_count) pfree(bin_width);
1560 
1561  MemoryContextSwitchTo(oldcontext);
1562  SRF_RETURN_DONE(funcctx);
1563  }
1564 
1565  /* we need the raw values, hence the non-zero parameter */
1566  stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 1, NULL, NULL, NULL);
1567 
1568  rt_band_destroy(band);
1569  rt_raster_destroy(raster);
1570 
1571  if (NULL == stats) {
1572  elog(NOTICE, "Cannot compute summary statistics for band at index %d. Returning NULL", bandindex);
1573 
1574  if (SPI_tuptable) SPI_freetuptable(tuptable);
1575  SPI_cursor_close(portal);
1576  SPI_finish();
1577 
1578  if (NULL != covhist) pfree(covhist);
1579  if (bin_width_count) pfree(bin_width);
1580 
1581  MemoryContextSwitchTo(oldcontext);
1582  SRF_RETURN_DONE(funcctx);
1583  }
1584 
1585  /* get histogram */
1586  if (stats->count > 0) {
1587  hist = rt_band_get_histogram(stats, bin_count, bin_width, bin_width_count, right, min, max, &count);
1588  pfree(stats);
1589  if (NULL == hist || !count) {
1590  elog(NOTICE, "Cannot compute histogram for band at index %d", bandindex);
1591 
1592  if (SPI_tuptable) SPI_freetuptable(tuptable);
1593  SPI_cursor_close(portal);
1594  SPI_finish();
1595 
1596  if (NULL != covhist) pfree(covhist);
1597  if (bin_width_count) pfree(bin_width);
1598 
1599  MemoryContextSwitchTo(oldcontext);
1600  SRF_RETURN_DONE(funcctx);
1601  }
1602 
1603  POSTGIS_RT_DEBUGF(3, "%d bins returned", count);
1604 
1605  /* coverage histogram */
1606  if (NULL == covhist) {
1607  covhist = (rt_histogram) SPI_palloc(sizeof(struct rt_histogram_t) * count);
1608  if (NULL == covhist) {
1609 
1610  pfree(hist);
1611  if (SPI_tuptable) SPI_freetuptable(tuptable);
1612  SPI_cursor_close(portal);
1613  SPI_finish();
1614 
1615  if (bin_width_count) pfree(bin_width);
1616 
1617  MemoryContextSwitchTo(oldcontext);
1618  elog(ERROR, "RASTER_histogramCoverage: Cannot allocate memory for histogram of coverage");
1619  SRF_RETURN_DONE(funcctx);
1620  }
1621 
1622  for (i = 0; i < count; i++) {
1623  sum += hist[i].count;
1624  covhist[i].count = hist[i].count;
1625  covhist[i].percent = 0;
1626  covhist[i].min = hist[i].min;
1627  covhist[i].max = hist[i].max;
1628  }
1629  }
1630  else {
1631  for (i = 0; i < count; i++) {
1632  sum += hist[i].count;
1633  covhist[i].count += hist[i].count;
1634  }
1635  }
1636 
1637  pfree(hist);
1638 
1639  /* assuming bin_count wasn't set, force consistency */
1640  if (bin_count <= 0) bin_count = count;
1641  }
1642 
1643  /* next record */
1644  SPI_cursor_fetch(portal, TRUE, 1);
1645  }
1646 
1647  if (SPI_tuptable) SPI_freetuptable(tuptable);
1648  SPI_cursor_close(portal);
1649  SPI_finish();
1650 
1651  if (bin_width_count) pfree(bin_width);
1652 
1653  /* finish percent of histogram */
1654  if (sum > 0) {
1655  for (i = 0; i < count; i++)
1656  covhist[i].percent = covhist[i].count / (double) sum;
1657  }
1658 
1659  /* Store needed information */
1660  funcctx->user_fctx = covhist;
1661 
1662  /* total number of tuples to be returned */
1663  funcctx->max_calls = count;
1664 
1665  /* Build a tuple descriptor for our result type */
1666  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
1667  ereport(ERROR, (
1668  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1669  errmsg(
1670  "function returning record called in context "
1671  "that cannot accept type record"
1672  )
1673  ));
1674  }
1675 
1676  BlessTupleDesc(tupdesc);
1677  funcctx->tuple_desc = tupdesc;
1678 
1679  MemoryContextSwitchTo(oldcontext);
1680  }
1681 
1682  /* stuff done on every call of the function */
1683  funcctx = SRF_PERCALL_SETUP();
1684 
1685  call_cntr = funcctx->call_cntr;
1686  max_calls = funcctx->max_calls;
1687  tupdesc = funcctx->tuple_desc;
1688  covhist2 = funcctx->user_fctx;
1689 
1690  /* do when there is more left to send */
1691  if (call_cntr < max_calls) {
1692  Datum values[VALUES_LENGTH];
1693  bool nulls[VALUES_LENGTH];
1694  HeapTuple tuple;
1695  Datum result;
1696 
1697  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
1698 
1699  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
1700 
1701  values[0] = Float8GetDatum(covhist2[call_cntr].min);
1702  values[1] = Float8GetDatum(covhist2[call_cntr].max);
1703  values[2] = Int64GetDatum(covhist2[call_cntr].count);
1704  values[3] = Float8GetDatum(covhist2[call_cntr].percent);
1705 
1706  /* build a tuple */
1707  tuple = heap_form_tuple(tupdesc, values, nulls);
1708 
1709  /* make the tuple into a datum */
1710  result = HeapTupleGetDatum(tuple);
1711 
1712  SRF_RETURN_NEXT(funcctx, result);
1713  }
1714  /* do when there is no more left */
1715  else {
1716  pfree(covhist2);
1717  SRF_RETURN_DONE(funcctx);
1718  }
1719 }
#define VALUES_LENGTH
int rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:372
uint32_t count
Definition: librtcore.h:2311
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:2185
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition: rt_band.c:242
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition: rtpostgis.h:65
unsigned int uint32_t
Definition: uthash.h:78
#define MAX_INT_CHARLEN
Definition: rtpostgis.h:76
double percent
Definition: librtcore.h:2326
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
#define MAX_DBL_CHARLEN
Definition: rtpostgis.h:75
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.
struct rt_histogram_t * rt_histogram
Definition: librtcore.h:151
#define FALSE
Definition: dbfopen.c:168
Struct definitions.
Definition: librtcore.h:2201
rt_histogram rt_band_get_histogram(rt_bandstats stats, int bin_count, double *bin_widths, int bin_widths_count, int right, double min, double max, uint32_t *rtn_count)
Count the distribution of data.
#define POSTGIS_RT_DEBUG(level, msg)
Definition: rtpostgis.h: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
uint32_t count
Definition: librtcore.h:2325
Here is the call graph for this function:
Here is the caller graph for this function: