PostGIS  3.4.0dev-r@@SVN_REVISION@@

◆ RASTER_quantile()

Datum RASTER_quantile ( PG_FUNCTION_ARGS  )

Definition at line 1196 of file rtpg_statistics.c.

1197 {
1198  FuncCallContext *funcctx;
1199  TupleDesc tupdesc;
1200 
1201  int i;
1202  rt_quantile quant;
1203  rt_quantile quant2;
1204  int call_cntr;
1205  int max_calls;
1206 
1207  /* first call of function */
1208  if (SRF_IS_FIRSTCALL()) {
1209  MemoryContext oldcontext;
1210 
1211  rt_pgraster *pgraster = NULL;
1212  rt_raster raster = NULL;
1213  rt_band band = NULL;
1214  int32_t bandindex = 0;
1215  int num_bands = 0;
1216  bool exclude_nodata_value = TRUE;
1217  double sample = 0;
1218  double *quantiles = NULL;
1219  uint32_t quantiles_count = 0;
1220  double quantile = 0;
1221  rt_bandstats stats = NULL;
1222  uint32_t count;
1223 
1224  int j;
1225  int n;
1226 
1227  ArrayType *array;
1228  Oid etype;
1229  Datum *e;
1230  bool *nulls;
1231  int16 typlen;
1232  bool typbyval;
1233  char typalign;
1234 
1235  /* create a function context for cross-call persistence */
1236  funcctx = SRF_FIRSTCALL_INIT();
1237 
1238  /* switch to memory context appropriate for multiple function calls */
1239  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
1240 
1241  /* pgraster is null, return nothing */
1242  if (PG_ARGISNULL(0)) {
1243  MemoryContextSwitchTo(oldcontext);
1244  SRF_RETURN_DONE(funcctx);
1245  }
1246  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1247 
1248  raster = rt_raster_deserialize(pgraster, FALSE);
1249  if (!raster) {
1250  PG_FREE_IF_COPY(pgraster, 0);
1251  MemoryContextSwitchTo(oldcontext);
1252  elog(ERROR, "RASTER_quantile: Cannot deserialize raster");
1253  SRF_RETURN_DONE(funcctx);
1254  }
1255 
1256  /* band index is 1-based */
1257  bandindex = PG_GETARG_INT32(1);
1258  num_bands = rt_raster_get_num_bands(raster);
1259  if (bandindex < 1 || bandindex > num_bands) {
1260  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
1261  rt_raster_destroy(raster);
1262  PG_FREE_IF_COPY(pgraster, 0);
1263  MemoryContextSwitchTo(oldcontext);
1264  SRF_RETURN_DONE(funcctx);
1265  }
1266 
1267  /* exclude_nodata_value flag */
1268  if (!PG_ARGISNULL(2))
1269  exclude_nodata_value = PG_GETARG_BOOL(2);
1270 
1271  /* sample % */
1272  if (!PG_ARGISNULL(3)) {
1273  sample = PG_GETARG_FLOAT8(3);
1274  if (sample < 0 || sample > 1) {
1275  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
1276  rt_raster_destroy(raster);
1277  PG_FREE_IF_COPY(pgraster, 0);
1278  MemoryContextSwitchTo(oldcontext);
1279  SRF_RETURN_DONE(funcctx);
1280  }
1281  else if (FLT_EQ(sample, 0.0))
1282  sample = 1;
1283  }
1284  else
1285  sample = 1;
1286 
1287  /* quantiles */
1288  if (!PG_ARGISNULL(4)) {
1289  array = PG_GETARG_ARRAYTYPE_P(4);
1290  etype = ARR_ELEMTYPE(array);
1291  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1292 
1293  switch (etype) {
1294  case FLOAT4OID:
1295  case FLOAT8OID:
1296  break;
1297  default:
1298  rt_raster_destroy(raster);
1299  PG_FREE_IF_COPY(pgraster, 0);
1300  MemoryContextSwitchTo(oldcontext);
1301  elog(ERROR, "RASTER_quantile: Invalid data type for quantiles");
1302  SRF_RETURN_DONE(funcctx);
1303  break;
1304  }
1305 
1306  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1307  &nulls, &n);
1308 
1309  quantiles = palloc(sizeof(double) * n);
1310  for (i = 0, j = 0; i < n; i++) {
1311  if (nulls[i]) continue;
1312 
1313  switch (etype) {
1314  case FLOAT4OID:
1315  quantile = (double) DatumGetFloat4(e[i]);
1316  break;
1317  case FLOAT8OID:
1318  quantile = (double) DatumGetFloat8(e[i]);
1319  break;
1320  }
1321 
1322  if (quantile < 0 || quantile > 1) {
1323  elog(NOTICE, "Invalid value for quantile (must be between 0 and 1). Returning NULL");
1324  pfree(quantiles);
1325  rt_raster_destroy(raster);
1326  PG_FREE_IF_COPY(pgraster, 0);
1327  MemoryContextSwitchTo(oldcontext);
1328  SRF_RETURN_DONE(funcctx);
1329  }
1330 
1331  quantiles[j] = quantile;
1332  POSTGIS_RT_DEBUGF(5, "quantiles[%d] = %f", j, quantiles[j]);
1333  j++;
1334  }
1335  quantiles_count = j;
1336 
1337  if (j < 1) {
1338  pfree(quantiles);
1339  quantiles = NULL;
1340  }
1341  }
1342 
1343  /* get band */
1344  band = rt_raster_get_band(raster, bandindex - 1);
1345  if (!band) {
1346  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
1347  rt_raster_destroy(raster);
1348  PG_FREE_IF_COPY(pgraster, 0);
1349  MemoryContextSwitchTo(oldcontext);
1350  SRF_RETURN_DONE(funcctx);
1351  }
1352 
1353  /* get stats */
1354  stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 1, NULL, NULL, NULL);
1355  rt_band_destroy(band);
1356  rt_raster_destroy(raster);
1357  PG_FREE_IF_COPY(pgraster, 0);
1358  if (NULL == stats || NULL == stats->values) {
1359  elog(NOTICE, "Cannot retrieve summary statistics for band at index %d", bandindex);
1360  MemoryContextSwitchTo(oldcontext);
1361  SRF_RETURN_DONE(funcctx);
1362  }
1363  else if (stats->count < 1) {
1364  elog(NOTICE, "Cannot compute quantiles for band at index %d as the band has no values", bandindex);
1365  MemoryContextSwitchTo(oldcontext);
1366  SRF_RETURN_DONE(funcctx);
1367  }
1368 
1369  /* get quantiles */
1370  quant = rt_band_get_quantiles(stats, quantiles, quantiles_count, &count);
1371  if (quantiles_count) pfree(quantiles);
1372  pfree(stats);
1373  if (NULL == quant || !count) {
1374  elog(NOTICE, "Cannot compute quantiles for band at index %d", bandindex);
1375  MemoryContextSwitchTo(oldcontext);
1376  SRF_RETURN_DONE(funcctx);
1377  }
1378 
1379  POSTGIS_RT_DEBUGF(3, "%d quantiles returned", count);
1380 
1381  /* Store needed information */
1382  funcctx->user_fctx = quant;
1383 
1384  /* total number of tuples to be returned */
1385  funcctx->max_calls = count;
1386 
1387  /* Build a tuple descriptor for our result type */
1388  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
1389  ereport(ERROR, (
1390  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1391  errmsg(
1392  "function returning record called in context "
1393  "that cannot accept type record"
1394  )
1395  ));
1396  }
1397 
1398  BlessTupleDesc(tupdesc);
1399  funcctx->tuple_desc = tupdesc;
1400 
1401  MemoryContextSwitchTo(oldcontext);
1402  }
1403 
1404  /* stuff done on every call of the function */
1405  funcctx = SRF_PERCALL_SETUP();
1406 
1407  call_cntr = funcctx->call_cntr;
1408  max_calls = funcctx->max_calls;
1409  tupdesc = funcctx->tuple_desc;
1410  quant2 = funcctx->user_fctx;
1411 
1412  /* do when there is more left to send */
1413  if (call_cntr < max_calls) {
1414  Datum values[VALUES_LENGTH];
1415  bool nulls[VALUES_LENGTH];
1416  HeapTuple tuple;
1417  Datum result;
1418 
1419  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
1420 
1421  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
1422 
1423  values[0] = Float8GetDatum(quant2[call_cntr].quantile);
1424  values[1] = Float8GetDatum(quant2[call_cntr].value);
1425 
1426  /* build a tuple */
1427  tuple = heap_form_tuple(tupdesc, values, nulls);
1428 
1429  /* make the tuple into a datum */
1430  result = HeapTupleGetDatum(tuple);
1431 
1432  SRF_RETURN_NEXT(funcctx, result);
1433  }
1434  /* do when there is no more left */
1435  else {
1436  pfree(quant2);
1437  SRF_RETURN_DONE(funcctx);
1438  }
1439 }
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_band_get_quantiles
rt_quantile rt_band_get_quantiles(rt_bandstats stats, double *quantiles, int quantiles_count, uint32_t *rtn_count)
Compute the default set of or requested quantiles for a set of data the quantile formula used is same...
Definition: rt_statistics.c:672
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_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
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_DEBUGF
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition: rtpostgis.h:69
rt_band_t
Definition: librtcore.h:2466
rt_bandstats_t::count
uint32_t count
Definition: librtcore.h:2513
rt_bandstats_t::values
double * values
Definition: librtcore.h:2521
rt_bandstats_t
Definition: librtcore.h:2511
rt_quantile_t
Definition: librtcore.h:2538
rt_raster_serialized_t
Struct definitions.
Definition: librtcore.h:2403
rt_raster_t
Definition: librtcore.h:2437

References ovdump::band, rt_bandstats_t::count, genraster::count, FALSE, FLT_EQ, POSTGIS_RT_DEBUGF, rtrowdump::raster, result, rt_band_destroy(), rt_band_get_quantiles(), rt_band_get_summary_stats(), rt_raster_deserialize(), rt_raster_destroy(), rt_raster_get_band(), rt_raster_get_num_bands(), TRUE, genraster::value, rt_bandstats_t::values, and VALUES_LENGTH.

Here is the call graph for this function: