PostGIS  3.0.6dev-r@@SVN_REVISION@@

◆ RASTER_histogramCoverage()

Datum RASTER_histogramCoverage ( PG_FUNCTION_ARGS  )

Definition at line 1195 of file rtpg_statistics.c.

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

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

Here is the call graph for this function: