From: Bborie Park Date: Mon, 16 May 2011 19:39:41 +0000 (+0000) Subject: Added ST_Quantile functions X-Git-Tag: 2.0.0alpha1~1658 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=581c2dd9c3ac1bd12e746748b2a8abb2109fc1ed;p=postgis Added ST_Quantile functions - added function rt_raster_get_quantiles to rt_core/rt_api.c and rt_api.h - added test case to test/core/testapi.c - added function RASTER_quantile to rt_pg/rt_pg.c - added SQL functions for ST_Quantile - added regression tests Associated ticket is #935 git-svn-id: http://svn.osgeo.org/postgis/trunk@7153 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c index ac126ab3e..fdbbe8f11 100644 --- a/raster/rt_core/rt_api.c +++ b/raster/rt_core/rt_api.c @@ -127,6 +127,61 @@ rt_util_clamp_to_32F(double value) { return (float)fmin(fmax((value), -FLT_MAX), FLT_MAX); } +/* quicksort */ +#define swap(x, y) { double t; t = x; x = y; y = t; } +#define order(x, y) if (x > y) swap(x, y) + +static double pivot(double *left, double *right) { + double l, m, r, *p; + + l = *left; + m = *(left + (right - left) / 2); + r = *right; + + /* order */ + order(l, m); + order(l, r); + order(m, r); + + /* pivot is higher of two values */ + if (l < m) return m; + if (m < r) return r; + + /* find pivot that isn't left */ + for (p = left + 1; p <= right; ++p) { + if (*p != *left) + return (*p < *left) ? *left : *p; + } + + /* all values are same */ + return -1; +} + +static double *partition(double *left, double *right, double pivot) { + while (left <= right) { + while (*left < pivot) ++left; + while (*right >= pivot) --right; + + if (left < right) { + swap(*left, *right); + ++left; + --right; + } + } + + return left; +} + +static void quicksort(double *left, double *right) { + double p = pivot(left, right); + double *pos; + + if (p != -1) { + pos = partition(left, right, p); + quicksort(left, pos - 1); + quicksort(pos, right); + } +} /*- rt_context -------------------------------------------------------*/ @@ -1948,6 +2003,122 @@ rt_band_get_histogram(rt_bandstats stats, return bins; } +struct rt_quantile_t { + double quantile; + double value; +}; + +/** + * Compute the default set of or requested quantiles for a set of data + * the quantile formula used is same as Excel and R default method + * + * @param stats: a populated stats struct for processing + * @param quantiles: the quantiles to be computed + * @param quantiles_count: the number of quantiles to be computed + * @param rtn_count: set to the number of quantiles being returned + * + * @return the default set of or requested quantiles for a band + */ +rt_quantile +rt_band_get_quantiles(rt_bandstats stats, + double *quantiles, int quantiles_count, int *rtn_count) { + rt_quantile rtn; + int init_quantiles = 0; + int i = 0; + double h; + int hl; + +#if POSTGIS_DEBUG_LEVEL > 0 + clock_t start, stop; + double elapsed = 0; +#endif + + RASTER_DEBUG(3, "starting"); +#if POSTGIS_DEBUG_LEVEL > 0 + start = clock(); +#endif + + assert(NULL != stats); + + if (stats->count < 1 || NULL == stats->values) { + rterror("rt_util_get_quantile: rt_bandstats object has no value"); + return NULL; + } + + /* quantiles not provided */ + if (NULL == quantiles) { + /* quantile count not specified, default to quartiles */ + if (quantiles_count < 2) + quantiles_count = 5; + + quantiles = rtalloc(sizeof(double) * quantiles_count); + init_quantiles = 1; + if (NULL == quantiles) { + rterror("rt_util_get_quantile: Unable to allocate memory for quantile input"); + return NULL; + } + + quantiles_count--; + for (i = 0; i <= quantiles_count; i++) + quantiles[i] = ((double) i) / quantiles_count; + quantiles_count++; + } + + /* check quantiles */ + for (i = 0; i < quantiles_count; i++) { + if (quantiles[i] < 0. || quantiles[i] > 1.) { + rterror("rt_util_get_quantile: Quantile value not between 0 and 1"); + if (init_quantiles) rtdealloc(quantiles); + return NULL; + } + } + quicksort(quantiles, quantiles + quantiles_count - 1); + + /* initialize rt_quantile */ + rtn = rtalloc(sizeof(struct rt_quantile_t) * quantiles_count); + if (NULL == rtn) { + rterror("rt_util_get_quantile: Unable to allocate memory for quantile output"); + if (init_quantiles) rtdealloc(quantiles); + return NULL; + } + + /* sort values */ + if (!stats->sorted) { + quicksort(stats->values, stats->values + stats->count - 1); + stats->sorted = 1; + } + + /* + make quantiles + + formula is that used in R (method 7) and Excel from + http://en.wikipedia.org/wiki/Quantile + */ + for (i = 0; i < quantiles_count; i++) { + rtn[i].quantile = quantiles[i]; + + h = ((stats->count - 1.) * quantiles[i]) + 1.; + hl = floor(h); + + /* h greater than hl, do full equation */ + if (h > hl) + rtn[i].value = stats->values[hl - 1] + ((h - hl) * (stats->values[hl] - stats->values[hl - 1])); + /* shortcut as second part of equation is zero */ + else + rtn[i].value = stats->values[hl - 1]; + } + +#if POSTGIS_DEBUG_LEVEL > 0 + stop = clock(); + elapsed = ((double) (stop - start)) / CLOCKS_PER_SEC; + RASTER_DEBUGF(3, "elapsed time = %0.4f", elapsed); +#endif + + *rtn_count = quantiles_count; + RASTER_DEBUG(3, "done"); + return rtn; +} + /*- rt_raster --------------------------------------------------------*/ struct rt_raster_serialized_t { diff --git a/raster/rt_core/rt_api.h b/raster/rt_core/rt_api.h index b975ba373..c97ad173d 100644 --- a/raster/rt_core/rt_api.h +++ b/raster/rt_core/rt_api.h @@ -474,6 +474,21 @@ rt_histogram rt_band_get_histogram(rt_bandstats stats, int bin_count, double *bin_widths, int bin_widths_count, int right, int *rtn_count); +/** + * Compute the default set of or requested quantiles for a set of data + * the quantile formula used is same as Excel and R default method + * + * @param stats: a populated stats struct for processing + * @param quantiles: the quantiles to be computed + * @param quantiles_count: the number of quantiles to be computed + * @param rtn_count: the number of quantiles being returned + * + * @return the default set of or requested quantiles for a band + */ +typedef struct rt_quantile_t* rt_quantile; +rt_quantile rt_band_get_quantiles(rt_bandstats stats, + double *quantiles, int quantiles_count, int *rtn_count); + /*- rt_raster --------------------------------------------------------*/ diff --git a/raster/rt_pg/rt_pg.c b/raster/rt_pg/rt_pg.c index 3e2adb8dc..f3baf0a7b 100644 --- a/raster/rt_pg/rt_pg.c +++ b/raster/rt_pg/rt_pg.c @@ -193,6 +193,9 @@ Datum RASTER_summaryStats(PG_FUNCTION_ARGS); /* get histogram */ Datum RASTER_histogram(PG_FUNCTION_ARGS); +/* get quantiles */ +Datum RASTER_quantile(PG_FUNCTION_ARGS); + /* Replace function taken from * http://ubuntuforums.org/showthread.php?s=aa6f015109fd7e4c7e30d2fd8b717497&t=141670&page=3 @@ -3268,6 +3271,272 @@ Datum RASTER_histogram(PG_FUNCTION_ARGS) } } +/* get quantiles */ +struct rt_quantile_t { + double quantile; + double value; +}; + +/** + * Returns quantiles for a band + */ +PG_FUNCTION_INFO_V1(RASTER_quantile); +Datum RASTER_quantile(PG_FUNCTION_ARGS) +{ + FuncCallContext *funcctx; + TupleDesc tupdesc; + AttInMetadata *attinmeta; + + int count; + rt_quantile quant; + rt_quantile quant2; + int call_cntr; + int max_calls; + + /* first call of function */ + if (SRF_IS_FIRSTCALL()) { + MemoryContext oldcontext; + + rt_pgraster *pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + rt_raster raster = NULL; + rt_band band = NULL; + int32_t bandindex = 0; + int num_bands = 0; + bool hasnodata = FALSE; + double sample = 0; + double *quantiles = NULL; + int quantiles_count = 0; + double quantile; + rt_bandstats stats = NULL; + + int i; + int j; + int n; + + ArrayType *array; + Oid etype; + Datum *e; + bool *nulls; + int16 typlen; + bool typbyval; + char typalign; + int ndims = 1; + int *dims; + int *lbs; + + /* create a function context for cross-call persistence */ + funcctx = SRF_FIRSTCALL_INIT(); + + /* switch to memory context appropriate for multiple function calls */ + oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx); + + raster = rt_raster_deserialize(pgraster); + if (!raster) { + elog(ERROR, "RASTER_quantile: Could not deserialize raster"); + PG_RETURN_NULL(); + } + + /* band index is 1-based */ + bandindex = PG_GETARG_INT32(1); + num_bands = rt_raster_get_num_bands(raster); + if (bandindex < 1 || bandindex > num_bands) { + elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL"); + rt_raster_destroy(raster); + SRF_RETURN_DONE(funcctx); + } + assert(0 <= (bandindex - 1)); + + /* hasnodata flag */ + if (!PG_ARGISNULL(2)) + hasnodata = PG_GETARG_BOOL(2); + + /* sample % */ + if (!PG_ARGISNULL(3)) { + sample = PG_GETARG_FLOAT8(3); + if (sample < 0 || sample > 1) { + elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL"); + rt_raster_destroy(raster); + SRF_RETURN_DONE(funcctx); + } + else if (sample == 0) + sample = 1; + } + else + sample = 1; + + /* quantiles */ + if (!PG_ARGISNULL(4)) { + array = PG_GETARG_ARRAYTYPE_P(4); + etype = ARR_ELEMTYPE(array); + get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign); + + switch (etype) { + case FLOAT4OID: + case FLOAT8OID: + break; + default: + elog(ERROR, "RASTER_quantile: Invalid data type for quanitiles"); + rt_raster_destroy(raster); + PG_RETURN_NULL(); + break; + } + + ndims = ARR_NDIM(array); + dims = ARR_DIMS(array); + lbs = ARR_LBOUND(array); + + deconstruct_array(array, etype, typlen, typbyval, typalign, &e, + &nulls, &n); + + quantiles = palloc(sizeof(double) * n); + for (i = 0, j = 0; i < n; i++) { + if (nulls[i]) continue; + + switch (etype) { + case FLOAT4OID: + quantile = (double) DatumGetFloat4(e[i]); + break; + case FLOAT8OID: + quantile = (double) DatumGetFloat8(e[i]); + break; + } + + if (quantile < 0 || quantile > 1) { + elog(NOTICE, "Invalid value for quantile (must be between 0 and 1). Returning NULL"); + pfree(quantiles); + rt_raster_destroy(raster); + SRF_RETURN_DONE(funcctx); + } + + quantiles[j] = quantile; + POSTGIS_RT_DEBUGF(5, "quantiles[%d] = %f", j, quantiles[j]); + j++; + } + quantiles_count = j; + + if (j < 1) { + pfree(quantiles); + quantiles = NULL; + } + } + + /* get band */ + band = rt_raster_get_band(raster, bandindex - 1); + if (!band) { + elog(NOTICE, "Could not find raster band of index %d. Returning NULL", bandindex); + rt_raster_destroy(raster); + SRF_RETURN_DONE(funcctx); + } + + /* get stats */ + stats = rt_band_get_summary_stats(band, (int) hasnodata, sample, 1); + rt_band_destroy(band); + rt_raster_destroy(raster); + if (NULL == stats || NULL == stats->values) { + elog(NOTICE, "Could not retrieve summary statistics of raster band of index %d", bandindex); + SRF_RETURN_DONE(funcctx); + } + + /* get quantiles */ + quant = rt_band_get_quantiles(stats, quantiles, quantiles_count, &count); + if (quantiles_count) pfree(quantiles); + pfree(stats); + if (NULL == quant || !count) { + elog(NOTICE, "Could not retrieve quantiles of raster band of index %d", bandindex); + SRF_RETURN_DONE(funcctx); + } + + POSTGIS_RT_DEBUGF(3, "%d quantiles returned", count); + + /* Store needed information */ + funcctx->user_fctx = quant; + + /* total number of tuples to be returned */ + funcctx->max_calls = count; + + /* Build a tuple descriptor for our result type */ + if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) { + ereport(ERROR, ( + errcode(ERRCODE_FEATURE_NOT_SUPPORTED), + errmsg( + "function returning record called in context " + "that cannot accept type record" + ) + )); + } + + /* + * generate attribute metadata needed later to produce tuples from raw + * C strings + */ + attinmeta = TupleDescGetAttInMetadata(tupdesc); + funcctx->attinmeta = attinmeta; + MemoryContextSwitchTo(oldcontext); + } + + /* stuff done on every call of the function */ + funcctx = SRF_PERCALL_SETUP(); + + call_cntr = funcctx->call_cntr; + max_calls = funcctx->max_calls; + attinmeta = funcctx->attinmeta; + quant2 = funcctx->user_fctx; + + /* do when there is more left to send */ + if (call_cntr < max_calls) { + char **values; + HeapTuple tuple; + Datum result; + + POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr); + + /* + * Prepare a values array for building the returned tuple. + * This should be an array of C strings which will + * be processed later by the type input functions. + */ + values = (char **) palloc(2 * sizeof(char *)); + + values[0] = (char *) palloc( + (snprintf(NULL, 0, "%f", quant2[call_cntr].quantile) + 1) * sizeof(char) + ); + values[1] = (char *) palloc( + (snprintf(NULL, 0, "%f", quant2[call_cntr].value) + 1) * sizeof(char) + ); + + snprintf( + values[0], + (snprintf(NULL, 0, "%f", quant2[call_cntr].quantile) + 1) * sizeof(char), + "%f", + quant2[call_cntr].quantile + ); + snprintf( + values[1], + (snprintf(NULL, 0, "%f", quant2[call_cntr].value) + 1) * sizeof(char), + "%f", + quant2[call_cntr].value + ); + + /* build a tuple */ + tuple = BuildTupleFromCStrings(attinmeta, values); + + /* make the tuple into a datum */ + result = HeapTupleGetDatum(tuple); + + /* clean up (this is not really necessary) */ + pfree(values[1]); + pfree(values[0]); + pfree(values); + + SRF_RETURN_NEXT(funcctx, result); + } + /* do when there is no more left */ + else { + pfree(quant2); + SRF_RETURN_DONE(funcctx); + } +} + /* ---------------------------------------------------------------- */ /* Memory allocation / error reporting hooks */ /* ---------------------------------------------------------------- */ diff --git a/raster/rt_pg/rtpostgis.sql.in.c b/raster/rt_pg/rtpostgis.sql.in.c index e72cee456..503e2ed0d 100644 --- a/raster/rt_pg/rtpostgis.sql.in.c +++ b/raster/rt_pg/rtpostgis.sql.in.c @@ -849,6 +849,130 @@ CREATE OR REPLACE FUNCTION st_approxhistogram(rast raster, nband int, sample_per AS $$ SELECT min, max, count, proportion FROM _st_histogram($1, $2, FALSE, $3, 0, NULL, FALSE) $$ LANGUAGE 'sql' IMMUTABLE STRICT; +----------------------------------------------------------------------- +-- ST_Quantile and ST_ApproxQuantile +----------------------------------------------------------------------- +CREATE TYPE quantile AS ( + quantile double precision, + value double precision +); + +-- Cannot be strict as "quantile" can be NULL +CREATE OR REPLACE FUNCTION _st_quantile(rast raster, nband int, hasnodata boolean, sample_percent double precision, quantiles double precision[]) + RETURNS SETOF quantile + AS 'MODULE_PATHNAME','RASTER_quantile' + LANGUAGE 'C' IMMUTABLE; + +CREATE OR REPLACE FUNCTION st_quantile(rast raster, nband int, hasnodata boolean, quantiles double precision[]) + RETURNS SETOF quantile + AS $$ SELECT quantile, value FROM _st_quantile($1, $2, $3, 1, $4) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_quantile(rast raster, nband int, quantiles double precision[]) + RETURNS SETOF quantile + AS $$ SELECT quantile, value FROM _st_quantile($1, $2, FALSE, 1, $3) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_quantile(rast raster, nband int, hasnodata boolean) + RETURNS SETOF quantile + AS $$ SELECT quantile, value FROM _st_quantile($1, $2, $3, 1, NULL) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_quantile(rast raster, nband int) + RETURNS SETOF quantile + AS $$ SELECT quantile, value FROM _st_quantile($1, $2, FALSE, 1, NULL) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_quantile(rast raster, quantiles double precision[]) + RETURNS SETOF quantile + AS $$ SELECT quantile, value FROM _st_quantile($1, 1, FALSE, 1, $2) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_quantile(rast raster) + RETURNS SETOF quantile + AS $$ SELECT quantile, value FROM _st_quantile($1, 1, FALSE, 1, NULL) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_quantile(rast raster, nband int, hasnodata boolean, quantile double precision) + RETURNS quantile + AS $$ SELECT quantile, value FROM _st_quantile($1, $2, $3, 1, ARRAY[$4]::double precision[]) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_quantile(rast raster, nband int, quantile double precision) + RETURNS quantile + AS $$ SELECT quantile, value FROM _st_quantile($1, $2, FALSE, 1, ARRAY[$3]::double precision[]) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_quantile(rast raster, hasnodata boolean, quantile double precision) + RETURNS quantile + AS $$ SELECT quantile, value FROM _st_quantile($1, 1, $2, 1, ARRAY[$3]::double precision[]) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_quantile(rast raster, quantile double precision) + RETURNS quantile + AS $$ SELECT quantile, value FROM _st_quantile($1, 1, FALSE, 1, ARRAY[$2]::double precision[]) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxquantile(rast raster, nband int, hasnodata boolean, sample_percent double precision, quantiles double precision[]) + RETURNS SETOF quantile + AS $$ SELECT quantile, value FROM _st_quantile($1, $2, $3, $4, $5) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxquantile(rast raster, nband int, sample_percent double precision, quantiles double precision[]) + RETURNS SETOF quantile + AS $$ SELECT quantile, value FROM _st_quantile($1, $2, FALSE, $3, $4) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxquantile(rast raster, nband int, sample_percent double precision) + RETURNS SETOF quantile + AS $$ SELECT quantile, value FROM _st_quantile($1, $2, FALSE, $3, NULL) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxquantile(rast raster, sample_percent double precision, quantiles double precision[]) + RETURNS SETOF quantile + AS $$ SELECT quantile, value FROM _st_quantile($1, 1, FALSE, $2, $3) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxquantile(rast raster, sample_percent double precision) + RETURNS SETOF quantile + AS $$ SELECT quantile, value FROM _st_quantile($1, 1, FALSE, $2, NULL) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxquantile(rast raster, quantiles double precision[]) + RETURNS SETOF quantile + AS $$ SELECT quantile, value FROM _st_quantile($1, 1, FALSE, 0.1, $2) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxquantile(rast raster, nband int) + RETURNS SETOF quantile + AS $$ SELECT quantile, value FROM _st_quantile($1, $2, FALSE, 0.1, NULL) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxquantile(rast raster) + RETURNS SETOF quantile + AS $$ SELECT quantile, value FROM _st_quantile($1, 1, FALSE, 0.1, NULL) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxquantile(rast raster, nband int, hasnodata boolean, sample_percent double precision, quantile double precision) + RETURNS quantile + AS $$ SELECT quantile, value FROM _st_quantile($1, $2, $3, $4, ARRAY[$5]::double precision[]) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxquantile(rast raster, nband int, sample_percent double precision, quantile double precision) + RETURNS quantile + AS $$ SELECT quantile, value FROM _st_quantile($1, $2, FALSE, $3, ARRAY[$4]::double precision[]) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxquantile(rast raster, sample_percent double precision, quantile double precision) + RETURNS quantile + AS $$ SELECT quantile, value FROM _st_quantile($1, 1, FALSE, $2, ARRAY[$3]::double precision[]) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxquantile(rast raster, hasnodata boolean, quantile double precision) + RETURNS quantile + AS $$ SELECT quantile, value FROM _st_quantile($1, 1, $2, 1, ARRAY[$3]::double precision[]) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + ----------------------------------------------------------------------- -- MapAlgebra ----------------------------------------------------------------------- diff --git a/raster/test/core/testapi.c b/raster/test/core/testapi.c index 320702856..c878ec8b5 100644 --- a/raster/test/core/testapi.c +++ b/raster/test/core/testapi.c @@ -1009,10 +1009,16 @@ struct rt_histogram_t { int inc_min; int inc_max; }; +struct rt_quantile_t { + double quantile; + double value; +}; static void testBandStats() { rt_bandstats stats = NULL; rt_histogram histogram = NULL; double bin_width[] = {10000}; + double quantiles[] = {0.1, 0.3, 0.5, 0.7, 0.9}; + rt_quantile quantile = NULL; int count = 0; rt_raster raster; @@ -1045,6 +1051,10 @@ static void testBandStats() { CHECK_EQUALS(stats->min, 1); CHECK_EQUALS(stats->max, 19998); + quantile = (rt_quantile) rt_band_get_quantiles(stats, NULL, 0, &count); + CHECK(quantile); + rtdealloc(quantile); + histogram = (rt_histogram) rt_band_get_histogram(stats, 0, NULL, 0, 0, &count); CHECK(histogram); rtdealloc(histogram); @@ -1063,6 +1073,15 @@ static void testBandStats() { stats = (rt_bandstats) rt_band_get_summary_stats(band, 0, 0.1, 1); CHECK(stats); + quantile = (rt_quantile) rt_band_get_quantiles(stats, NULL, 0, &count); + CHECK(quantile); + rtdealloc(quantile); + + quantile = (rt_quantile) rt_band_get_quantiles(stats, quantiles, 5, &count); + CHECK(quantile); + CHECK((count == 5)); + rtdealloc(quantile); + histogram = (rt_histogram) rt_band_get_histogram(stats, 0, NULL, 0, 0, &count); CHECK(histogram); rtdealloc(histogram); @@ -1087,12 +1106,20 @@ static void testBandStats() { CHECK_EQUALS(stats->min, 0); CHECK_EQUALS(stats->max, 19998); + quantile = (rt_quantile) rt_band_get_quantiles(stats, NULL, 0, &count); + CHECK(quantile); + rtdealloc(quantile); + rtdealloc(stats->values); rtdealloc(stats); stats = (rt_bandstats) rt_band_get_summary_stats(band, 1, 0.1, 1); CHECK(stats); + quantile = (rt_quantile) rt_band_get_quantiles(stats, NULL, 0, &count); + CHECK(quantile); + rtdealloc(quantile); + rtdealloc(stats->values); rtdealloc(stats); diff --git a/raster/test/regress/Makefile.in b/raster/test/regress/Makefile.in index ac80ebef8..c9ed69ea5 100644 --- a/raster/test/regress/Makefile.in +++ b/raster/test/regress/Makefile.in @@ -66,6 +66,7 @@ TEST_BANDPROPS = \ rt_stddev.sql \ rt_minmax.sql \ rt_histogram.sql \ + rt_quantile.sql \ $(NULL) TEST_PIXEL = \ diff --git a/raster/test/regress/rt_quantile.sql b/raster/test/regress/rt_quantile.sql new file mode 100644 index 000000000..b0204f0c3 --- /dev/null +++ b/raster/test/regress/rt_quantile.sql @@ -0,0 +1,159 @@ +SELECT * FROM ST_Quantile( + ST_SetValue( + ST_SetValue( + ST_SetValue( + ST_AddBand( + ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1) + , 1, '64BF', 0, 0 + ) + , 1, 1, 1, -10 + ) + , 1, 5, 4, 0 + ) + , 1, 5, 5, 3.14159 + ) + , 1, TRUE, ARRAY[0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1]::double precision[] +); +SELECT * FROM ST_Quantile( + ST_SetValue( + ST_SetValue( + ST_SetValue( + ST_AddBand( + ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1) + , 1, '64BF', 0, 0 + ) + , 1, 1, 1, -10 + ) + , 1, 5, 4, 0 + ) + , 1, 5, 5, 3.14159 + ) + , 1, ARRAY[0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1]::double precision[] +); +SELECT * FROM ST_Quantile( + ST_SetValue( + ST_SetValue( + ST_SetValue( + ST_AddBand( + ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1) + , 1, '64BF', 0, 0 + ) + , 1, 1, 1, -10 + ) + , 1, 5, 4, 0 + ) + , 1, 5, 5, 3.14159 + ), + 1, TRUE +); +SELECT * FROM ST_Quantile( + ST_SetValue( + ST_SetValue( + ST_SetValue( + ST_AddBand( + ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1) + , 1, '64BF', 0, 0 + ) + , 1, 1, 1, -10 + ) + , 1, 5, 4, 0 + ) + , 1, 5, 5, 3.14159 + ), + 1 +); +SELECT * FROM ST_Quantile( + ST_SetValue( + ST_SetValue( + ST_SetValue( + ST_AddBand( + ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1) + , 1, '64BF', 0, 0 + ) + , 1, 1, 1, -10 + ) + , 1, 5, 4, 0 + ) + , 1, 5, 5, 3.14159 + ), + ARRAY[0.05, 0.95]::double precision[] +); +SELECT * FROM ST_Quantile( + ST_SetValue( + ST_SetValue( + ST_SetValue( + ST_AddBand( + ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1) + , 1, '64BF', 0, 0 + ) + , 1, 1, 1, -10 + ) + , 1, 5, 4, 0 + ) + , 1, 5, 5, 3.14159 + ) +); +SELECT * FROM ST_Quantile( + ST_SetValue( + ST_SetValue( + ST_SetValue( + ST_AddBand( + ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1) + , 1, '64BF', 0, 0 + ) + , 1, 1, 1, -10 + ) + , 1, 5, 4, 0 + ) + , 1, 5, 5, 3.14159 + ), + 1, TRUE, 0.05 +); +SELECT * FROM ST_Quantile( + ST_SetValue( + ST_SetValue( + ST_SetValue( + ST_AddBand( + ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1) + , 1, '64BF', 0, 0 + ) + , 1, 1, 1, -10 + ) + , 1, 5, 4, 0 + ) + , 1, 5, 5, 3.14159 + ), + 1, 0.95 +); +SELECT * FROM ST_Quantile( + ST_SetValue( + ST_SetValue( + ST_SetValue( + ST_AddBand( + ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1) + , 1, '64BF', 0, 0 + ) + , 1, 1, 1, -10 + ) + , 1, 5, 4, 0 + ) + , 1, 5, 5, 3.14159 + ), + TRUE, 0.7 +); +SELECT * FROM ST_Quantile( + ST_SetValue( + ST_SetValue( + ST_SetValue( + ST_AddBand( + ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1) + , 1, '64BF', 0, 0 + ) + , 1, 1, 1, -10 + ) + , 1, 5, 4, 0 + ) + , 1, 5, 5, 3.14159 + ), + 0.45 +); diff --git a/raster/test/regress/rt_quantile_expected b/raster/test/regress/rt_quantile_expected new file mode 100644 index 000000000..d612998b2 --- /dev/null +++ b/raster/test/regress/rt_quantile_expected @@ -0,0 +1,43 @@ +0|-10 +0.1|0 +0.2|0 +0.3|0 +0.4|0 +0.5|0 +0.6|0 +0.7|0 +0.8|0 +0.9|0 +1|3.14159 +0|-10 +0.1|-8.685841 +0.2|-7.371682 +0.3|-6.057523 +0.4|-4.743364 +0.5|-3.429205 +0.6|-2.115046 +0.7|-0.800887 +0.8|0.513272 +0.9|1.827431 +1|3.14159 +0|-10 +0.25|0 +0.5|0 +0.75|0 +1|3.14159 +0|-10 +0.25|-6.714602 +0.5|-3.429205 +0.75|-0.143807 +1|3.14159 +0.05|-9.34292 +0.95|2.484511 +0|-10 +0.25|-6.714602 +0.5|-3.429205 +0.75|-0.143807 +1|3.14159 +0.05|0 +0.95|2.484511 +0.7|0 +0.45|-4.086285