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

◆ 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
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:67
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 rt_bandstats_t::count, rt_histogram_t::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, 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(), text_to_cstring(), TRUE, and VALUES_LENGTH.

Here is the call graph for this function: