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

◆ RASTER_quantileCoverage()

Datum RASTER_quantileCoverage ( PG_FUNCTION_ARGS  )

Definition at line 1973 of file rtpg_statistics.c.

1974{
1975 FuncCallContext *funcctx;
1976 TupleDesc tupdesc;
1977
1978 uint32_t i;
1979 rt_quantile covquant = NULL;
1980 rt_quantile covquant2;
1981 int call_cntr;
1982 int max_calls;
1983
1984 POSTGIS_RT_DEBUG(3, "RASTER_quantileCoverage: Starting");
1985
1986 /* first call of function */
1987 if (SRF_IS_FIRSTCALL()) {
1988 MemoryContext oldcontext;
1989
1990 text *tablenametext = NULL;
1991 char *tablename = NULL;
1992 text *colnametext = NULL;
1993 char *colname = NULL;
1994 int32_t bandindex = 1;
1995 bool exclude_nodata_value = TRUE;
1996 double sample = 0;
1997 double *quantiles = NULL;
1998 uint32_t quantiles_count = 0;
1999 double quantile = 0;
2000 uint32_t count;
2001
2002 int len = 0;
2003 char *sql = NULL;
2004 char *tmp = NULL;
2005 uint64_t cov_count = 0;
2006 int spi_result;
2007 Portal portal;
2008 SPITupleTable *tuptable = NULL;
2009 HeapTuple tuple;
2010 Datum datum;
2011 bool isNull = FALSE;
2012
2013 rt_pgraster *pgraster = NULL;
2014 rt_raster raster = NULL;
2015 rt_band band = NULL;
2016 int num_bands = 0;
2017 struct quantile_llist *qlls = NULL;
2018 uint32_t qlls_count;
2019
2020 int j;
2021 int n;
2022
2023 ArrayType *array;
2024 Oid etype;
2025 Datum *e;
2026 bool *nulls;
2027 int16 typlen;
2028 bool typbyval;
2029 char typalign;
2030
2031 POSTGIS_RT_DEBUG(3, "RASTER_quantileCoverage: first call of function");
2032
2033 /* create a function context for cross-call persistence */
2034 funcctx = SRF_FIRSTCALL_INIT();
2035
2036 /* switch to memory context appropriate for multiple function calls */
2037 oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
2038
2039 /* tablename is null, return null */
2040 if (PG_ARGISNULL(0)) {
2041 elog(NOTICE, "Table name must be provided");
2042 MemoryContextSwitchTo(oldcontext);
2043 SRF_RETURN_DONE(funcctx);
2044 }
2045 tablenametext = PG_GETARG_TEXT_P(0);
2046 tablename = text_to_cstring(tablenametext);
2047 if (!strlen(tablename)) {
2048 elog(NOTICE, "Table name must be provided");
2049 MemoryContextSwitchTo(oldcontext);
2050 SRF_RETURN_DONE(funcctx);
2051 }
2052 POSTGIS_RT_DEBUGF(3, "RASTER_quantileCoverage: tablename = %s", tablename);
2053
2054 /* column name is null, return null */
2055 if (PG_ARGISNULL(1)) {
2056 elog(NOTICE, "Column name must be provided");
2057 MemoryContextSwitchTo(oldcontext);
2058 SRF_RETURN_DONE(funcctx);
2059 }
2060 colnametext = PG_GETARG_TEXT_P(1);
2061 colname = text_to_cstring(colnametext);
2062 if (!strlen(colname)) {
2063 elog(NOTICE, "Column name must be provided");
2064 MemoryContextSwitchTo(oldcontext);
2065 SRF_RETURN_DONE(funcctx);
2066 }
2067 POSTGIS_RT_DEBUGF(3, "RASTER_quantileCoverage: colname = %s", colname);
2068
2069 /* band index is 1-based */
2070 if (!PG_ARGISNULL(2))
2071 bandindex = PG_GETARG_INT32(2);
2072
2073 /* exclude_nodata_value flag */
2074 if (!PG_ARGISNULL(3))
2075 exclude_nodata_value = PG_GETARG_BOOL(3);
2076
2077 /* sample % */
2078 if (!PG_ARGISNULL(4)) {
2079 sample = PG_GETARG_FLOAT8(4);
2080 if (sample < 0 || sample > 1) {
2081 elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
2082 MemoryContextSwitchTo(oldcontext);
2083 SRF_RETURN_DONE(funcctx);
2084 }
2085 else if (FLT_EQ(sample, 0.0))
2086 sample = 1;
2087 }
2088 else
2089 sample = 1;
2090
2091 /* quantiles */
2092 if (!PG_ARGISNULL(5)) {
2093 array = PG_GETARG_ARRAYTYPE_P(5);
2094 etype = ARR_ELEMTYPE(array);
2095 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
2096
2097 switch (etype) {
2098 case FLOAT4OID:
2099 case FLOAT8OID:
2100 break;
2101 default:
2102 MemoryContextSwitchTo(oldcontext);
2103 elog(ERROR, "RASTER_quantileCoverage: Invalid data type for quantiles");
2104 SRF_RETURN_DONE(funcctx);
2105 break;
2106 }
2107
2108 deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
2109 &nulls, &n);
2110
2111 quantiles = palloc(sizeof(double) * n);
2112 for (i = 0, j = 0; i < (uint32_t) n; i++) {
2113 if (nulls[i]) continue;
2114
2115 switch (etype) {
2116 case FLOAT4OID:
2117 quantile = (double) DatumGetFloat4(e[i]);
2118 break;
2119 case FLOAT8OID:
2120 quantile = (double) DatumGetFloat8(e[i]);
2121 break;
2122 }
2123
2124 if (quantile < 0 || quantile > 1) {
2125 elog(NOTICE, "Invalid value for quantile (must be between 0 and 1). Returning NULL");
2126 pfree(quantiles);
2127 MemoryContextSwitchTo(oldcontext);
2128 SRF_RETURN_DONE(funcctx);
2129 }
2130
2131 quantiles[j] = quantile;
2132 POSTGIS_RT_DEBUGF(5, "quantiles[%d] = %f", j, quantiles[j]);
2133 j++;
2134 }
2135 quantiles_count = j;
2136
2137 if (j < 1) {
2138 pfree(quantiles);
2139 quantiles = NULL;
2140 }
2141 }
2142
2143 /* coverage stats */
2144 /* connect to database */
2145 spi_result = SPI_connect();
2146 if (spi_result != SPI_OK_CONNECT) {
2147 MemoryContextSwitchTo(oldcontext);
2148 elog(ERROR, "RASTER_quantileCoverage: Cannot connect to database using SPI");
2149 SRF_RETURN_DONE(funcctx);
2150 }
2151
2152 len = sizeof(char) * (strlen("SELECT count FROM _st_summarystats('','',,::boolean,)") + strlen(tablename) + strlen(colname) + (MAX_INT_CHARLEN * 2) + MAX_DBL_CHARLEN + 1);
2153 sql = (char *) palloc(len);
2154 if (NULL == sql) {
2155
2156 if (SPI_tuptable) SPI_freetuptable(tuptable);
2157 SPI_finish();
2158
2159 MemoryContextSwitchTo(oldcontext);
2160 elog(ERROR, "RASTER_quantileCoverage: Cannot allocate memory for sql");
2161 SRF_RETURN_DONE(funcctx);
2162 }
2163
2164 /* get stats */
2165 snprintf(sql, len, "SELECT count FROM _st_summarystats('%s','%s',%d,%d::boolean,%f)", tablename, colname, bandindex, (exclude_nodata_value ? 1 : 0), sample);
2166 POSTGIS_RT_DEBUGF(3, "stats sql: %s", sql);
2167 spi_result = SPI_execute(sql, TRUE, 0);
2168 pfree(sql);
2169 if (spi_result != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
2170
2171 if (SPI_tuptable) SPI_freetuptable(tuptable);
2172 SPI_finish();
2173
2174 MemoryContextSwitchTo(oldcontext);
2175 elog(ERROR, "RASTER_quantileCoverage: Cannot get summary stats of coverage");
2176 SRF_RETURN_DONE(funcctx);
2177 }
2178
2179 tupdesc = SPI_tuptable->tupdesc;
2180 tuptable = SPI_tuptable;
2181 tuple = tuptable->vals[0];
2182
2183 tmp = SPI_getvalue(tuple, tupdesc, 1);
2184 if (NULL == tmp || !strlen(tmp)) {
2185
2186 if (SPI_tuptable) SPI_freetuptable(tuptable);
2187 SPI_finish();
2188
2189 MemoryContextSwitchTo(oldcontext);
2190 elog(ERROR, "RASTER_quantileCoverage: Cannot get summary stats of coverage");
2191 SRF_RETURN_DONE(funcctx);
2192 }
2193 cov_count = strtol(tmp, NULL, 10);
2194 POSTGIS_RT_DEBUGF(3, "covcount = %d", (int) cov_count);
2195 pfree(tmp);
2196
2197 /* iterate through rasters of coverage */
2198 /* create sql */
2199 len = sizeof(char) * (strlen("SELECT \"\" FROM \"\" WHERE \"\" IS NOT NULL") + (strlen(colname) * 2) + strlen(tablename) + 1);
2200 sql = (char *) palloc(len);
2201 if (NULL == sql) {
2202
2203 if (SPI_tuptable) SPI_freetuptable(tuptable);
2204 SPI_finish();
2205
2206 MemoryContextSwitchTo(oldcontext);
2207 elog(ERROR, "RASTER_quantileCoverage: Cannot allocate memory for sql");
2208 SRF_RETURN_DONE(funcctx);
2209 }
2210
2211 /* get cursor */
2212 snprintf(sql, len, "SELECT \"%s\" FROM \"%s\" WHERE \"%s\" IS NOT NULL", colname, tablename, colname);
2213 POSTGIS_RT_DEBUGF(3, "coverage sql: %s", sql);
2214 portal = SPI_cursor_open_with_args(
2215 "coverage",
2216 sql,
2217 0, NULL,
2218 NULL, NULL,
2219 TRUE, 0
2220 );
2221 pfree(sql);
2222
2223 /* process resultset */
2224 SPI_cursor_fetch(portal, TRUE, 1);
2225 while (SPI_processed == 1 && SPI_tuptable != NULL) {
2226 if (NULL != covquant) pfree(covquant);
2227
2228 tupdesc = SPI_tuptable->tupdesc;
2229 tuple = SPI_tuptable->vals[0];
2230
2231 datum = SPI_getbinval(tuple, tupdesc, 1, &isNull);
2232 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
2233 SPI_freetuptable(SPI_tuptable);
2234 SPI_cursor_close(portal);
2235 SPI_finish();
2236
2237 MemoryContextSwitchTo(oldcontext);
2238 elog(ERROR, "RASTER_quantileCoverage: Cannot get raster of coverage");
2239 SRF_RETURN_DONE(funcctx);
2240 }
2241 else if (isNull) {
2242 SPI_cursor_fetch(portal, TRUE, 1);
2243 continue;
2244 }
2245
2246 pgraster = (rt_pgraster *) PG_DETOAST_DATUM(datum);
2247
2248 raster = rt_raster_deserialize(pgraster, FALSE);
2249 if (!raster) {
2250
2251 SPI_freetuptable(SPI_tuptable);
2252 SPI_cursor_close(portal);
2253 SPI_finish();
2254
2255 MemoryContextSwitchTo(oldcontext);
2256 elog(ERROR, "RASTER_quantileCoverage: Cannot deserialize raster");
2257 SRF_RETURN_DONE(funcctx);
2258 }
2259
2260 /* inspect number of bands*/
2261 num_bands = rt_raster_get_num_bands(raster);
2262 if (bandindex < 1 || bandindex > num_bands) {
2263 elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
2264
2265 rt_raster_destroy(raster);
2266
2267 SPI_freetuptable(SPI_tuptable);
2268 SPI_cursor_close(portal);
2269 SPI_finish();
2270
2271 MemoryContextSwitchTo(oldcontext);
2272 SRF_RETURN_DONE(funcctx);
2273 }
2274
2275 /* get band */
2276 band = rt_raster_get_band(raster, bandindex - 1);
2277 if (!band) {
2278 elog(NOTICE, "Cannot find raster band of index %d. Returning NULL", bandindex);
2279
2280 rt_raster_destroy(raster);
2281
2282 SPI_freetuptable(SPI_tuptable);
2283 SPI_cursor_close(portal);
2284 SPI_finish();
2285
2286 MemoryContextSwitchTo(oldcontext);
2287 SRF_RETURN_DONE(funcctx);
2288 }
2289
2290 covquant = rt_band_get_quantiles_stream(
2291 band,
2292 exclude_nodata_value, sample, cov_count,
2293 &qlls, &qlls_count,
2294 quantiles, quantiles_count,
2295 &count
2296 );
2297
2298 rt_band_destroy(band);
2299 rt_raster_destroy(raster);
2300
2301 if (!covquant || !count) {
2302 elog(NOTICE, "Cannot compute quantiles for band at index %d", bandindex);
2303
2304 SPI_freetuptable(SPI_tuptable);
2305 SPI_cursor_close(portal);
2306 SPI_finish();
2307
2308 MemoryContextSwitchTo(oldcontext);
2309 SRF_RETURN_DONE(funcctx);
2310 }
2311
2312 /* next record */
2313 SPI_cursor_fetch(portal, TRUE, 1);
2314 }
2315
2316 covquant2 = SPI_palloc(sizeof(struct rt_quantile_t) * count);
2317 for (i = 0; i < count; i++) {
2318 covquant2[i].quantile = covquant[i].quantile;
2319 covquant2[i].has_value = covquant[i].has_value;
2320 if (covquant2[i].has_value)
2321 covquant2[i].value = covquant[i].value;
2322 }
2323
2324 pfree(covquant);
2325 quantile_llist_destroy(&qlls, qlls_count);
2326
2327 if (SPI_tuptable) SPI_freetuptable(SPI_tuptable);
2328 SPI_cursor_close(portal);
2329 SPI_finish();
2330
2331 if (quantiles_count) pfree(quantiles);
2332
2333 POSTGIS_RT_DEBUGF(3, "%d quantiles returned", count);
2334
2335 /* Store needed information */
2336 funcctx->user_fctx = covquant2;
2337
2338 /* total number of tuples to be returned */
2339 funcctx->max_calls = count;
2340
2341 /* Build a tuple descriptor for our result type */
2342 if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
2343 ereport(ERROR, (
2344 errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
2345 errmsg(
2346 "function returning record called in context "
2347 "that cannot accept type record"
2348 )
2349 ));
2350 }
2351
2352 BlessTupleDesc(tupdesc);
2353 funcctx->tuple_desc = tupdesc;
2354
2355 MemoryContextSwitchTo(oldcontext);
2356 }
2357
2358 /* stuff done on every call of the function */
2359 funcctx = SRF_PERCALL_SETUP();
2360
2361 call_cntr = funcctx->call_cntr;
2362 max_calls = funcctx->max_calls;
2363 tupdesc = funcctx->tuple_desc;
2364 covquant2 = funcctx->user_fctx;
2365
2366 /* do when there is more left to send */
2367 if (call_cntr < max_calls) {
2368 Datum values[VALUES_LENGTH];
2369 bool nulls[VALUES_LENGTH];
2370 HeapTuple tuple;
2371 Datum result;
2372
2373 POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
2374
2375 memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
2376
2377 values[0] = Float8GetDatum(covquant2[call_cntr].quantile);
2378 if (covquant2[call_cntr].has_value)
2379 values[1] = Float8GetDatum(covquant2[call_cntr].value);
2380 else
2381 nulls[1] = TRUE;
2382
2383 /* build a tuple */
2384 tuple = heap_form_tuple(tupdesc, values, nulls);
2385
2386 /* make the tuple into a datum */
2387 result = HeapTupleGetDatum(tuple);
2388
2389 SRF_RETURN_NEXT(funcctx, result);
2390 }
2391 /* do when there is no more left */
2392 else {
2393 POSTGIS_RT_DEBUG(3, "done");
2394 pfree(covquant2);
2395 SRF_RETURN_DONE(funcctx);
2396 }
2397}
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_get_quantiles_stream
rt_quantile rt_band_get_quantiles_stream(rt_band band, int exclude_nodata_value, double sample, uint64_t cov_count, struct quantile_llist **qlls, uint32_t *qlls_count, double *quantiles, uint32_t quantiles_count, uint32_t *rtn_count)
Compute the default set of or requested quantiles for a coverage.
Definition rt_statistics.c:1010
quantile_llist_destroy
int quantile_llist_destroy(struct quantile_llist **list, uint32_t list_count)
Definition rt_statistics.c:878
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_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
quantile_llist::count
uint32_t count
Definition librtcore.h:2400
quantile_llist::quantile
double quantile
Definition librtcore.h:2395
quantile_llist
Definition librtcore.h:2393
rt_band_t
Definition librtcore.h:2314
rt_quantile_t::value
double value
Definition librtcore.h:2388
rt_quantile_t::quantile
double quantile
Definition librtcore.h:2387
rt_quantile_t::has_value
uint32_t has_value
Definition librtcore.h:2389
rt_quantile_t
Definition librtcore.h:2386
rt_raster_serialized_t
Struct definitions.
Definition librtcore.h:2251
rt_raster_t
Definition librtcore.h:2285

References quantile_llist::count, FALSE, FLT_EQ, rt_quantile_t::has_value, MAX_DBL_CHARLEN, MAX_INT_CHARLEN, POSTGIS_RT_DEBUG, POSTGIS_RT_DEBUGF, rt_quantile_t::quantile, quantile_llist::quantile, quantile_llist_destroy(), rt_band_destroy(), rt_band_get_quantiles_stream(), rt_raster_deserialize(), rt_raster_destroy(), rt_raster_get_band(), rt_raster_get_num_bands(), text_to_cstring(), TRUE, rt_quantile_t::value, and VALUES_LENGTH.

Here is the call graph for this function: