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

◆ 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:267
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:2436
rt_band_destroy
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition rt_band.c:499
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:673
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:111
rt_raster_deserialize
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition rt_serialize.c:732
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.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:125
VALUES_LENGTH
#define VALUES_LENGTH
Definition rtpg_statistics.c:65
POSTGIS_RT_DEBUGF
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition rtpostgis.h:69
rt_band_t
Definition librtcore.h:2515
rt_bandstats_t::count
uint32_t count
Definition librtcore.h:2562
rt_bandstats_t::values
double * values
Definition librtcore.h:2570
rt_bandstats_t
Definition librtcore.h:2560
rt_quantile_t
Definition librtcore.h:2587
rt_raster_serialized_t
Struct definitions.
Definition librtcore.h:2452
rt_raster_t
Definition librtcore.h:2486

References rt_bandstats_t::count, FALSE, FLT_EQ, POSTGIS_RT_DEBUGF, 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, rt_bandstats_t::values, and VALUES_LENGTH.

Here is the call graph for this function: