PostGIS  2.5.7dev-r@@SVN_REVISION@@

◆ RASTER_histogramCoverage()

Datum RASTER_histogramCoverage ( PG_FUNCTION_ARGS  )

Definition at line 1191 of file rtpg_statistics.c.

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

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: