From: Bborie Park Date: Mon, 16 May 2011 19:17:15 +0000 (+0000) Subject: Addition of ST_SummaryStats function. X-Git-Tag: 2.0.0alpha1~1663 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=bffadd2c7615e4bc875701077d87406cacc8dec1;p=postgis Addition of ST_SummaryStats function. - added function rt_band_get_summary_stats to rt_core/rt_api.c and rt_api.h - added test case to test/core/testapi.c - added function RASTER_summaryStats to rt_pg/rt_pg.c - added SQL functions for ST_SummaryStats to rt_pg/rtpostgis.sql.in.c - added regression tests in test/regress Associated ticket is #930. git-svn-id: http://svn.osgeo.org/postgis/trunk@7148 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c index dd5d58309..cb01a5430 100644 --- a/raster/rt_core/rt_api.c +++ b/raster/rt_core/rt_api.c @@ -33,6 +33,7 @@ #include #include /* for FLT_EPSILON and float type limits */ #include /* for integer type limits */ +#include /* for time */ #include "rt_api.h" #define POSTGIS_RASTER_WARN_ON_TRUNCATION @@ -1425,6 +1426,285 @@ rt_band_check_is_nodata(rt_band band) return TRUE; } +struct rt_bandstats_t { + double sample; + uint32_t count; + + double min; + double max; + double mean; + double stddev; + + double *values; + int sorted; /* flag indicating that values is sorted ascending by value */ +}; + +/** + * Compute summary statistics for a band + * + * @param band: the band to query for minimum and maximum pixel values + * @param hasnodata: if zero, ignore nodata value + * @param sample: percentage of pixels to sample + * @param inc_vals: flag to include values in return struct + * + * @return the summary statistics for a band + */ +rt_bandstats +rt_band_get_summary_stats(rt_band band, int hasnodata, double sample, + int inc_vals) { + rt_pixtype pixtype = PT_END; + uint8_t *data = NULL; + uint32_t x = 0; + uint32_t y = 0; + uint32_t z = 0; + uint32_t offset = 0; + uint32_t diff = 0; + int rtn; + int hasnodata_flag = FALSE; + double nodata = 0; + double *values = NULL; + double value; + rt_bandstats stats = NULL; + + uint32_t sample_size = 0; + int byY = 0; + uint32_t outer = 0; + uint32_t inner = 0; + uint32_t sample_per = 0; + uint32_t sample_int = 0; + uint32_t i = 0; + uint32_t j = 0; + double sum = 0; + uint32_t k = 0; + double M = 0; + double Q = 0; + +#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 != band); + + if (band->offline) { + rterror("rt_band_get_summary_stats not implemented yet for OFFDB bands"); + return NULL; + } + + data = rt_band_get_data(band); + pixtype = band->pixtype; + + hasnodata_flag = rt_band_get_hasnodata_flag(band); + if (hasnodata_flag != FALSE) { + nodata = rt_band_get_nodata(band); + } + else { + hasnodata = 1; + } + + RASTER_DEBUGF(3, "nodata = %f", nodata); + RASTER_DEBUGF(3, "hasnodata_flag = %d", hasnodata_flag); + RASTER_DEBUGF(3, "user hasnodata = %d", hasnodata); + + /* entire band is nodata */ + if (rt_band_get_isnodata_flag(band) != FALSE) { + if (!hasnodata) { + rtwarn("All pixels of band have the NODATA value"); + return NULL; + } + else { + stats = (rt_bandstats) rtalloc(sizeof(struct rt_bandstats_t)); + if (NULL == stats) { + rterror("rt_band_get_summary_stats: Unable to allocate memory for stats"); + return NULL; + } + + stats->sample = 1; + stats->count = band->width * band->height; + stats->min = stats->max = nodata; + stats->mean = nodata; + stats->stddev = 0; + stats->values = NULL; + stats->sorted = 0; + + return stats; + } + } + + if (band->height > band->width) { + byY = 1; + outer = band->height; + inner = band->width; + } + else { + byY = 0; + outer = band->width; + inner = band->height; + } + + /* clamp percentage */ + if (sample <= 0 || sample > 1) + sample = 1; + + /* sample all pixels */ + if (sample == 1) { + sample_size = band->width * band->height; + sample_per = inner; + } + /* + randomly sample a percentage of available pixels + sampling method is known as + "systematic random sample without replacement" + */ + else { + sample_size = round((band->width * band->height) * sample); + sample_per = round(sample_size / outer); + sample_int = round(inner / sample_per); + srand(time(NULL)); + } + + RASTER_DEBUGF(3, "sampling %d of %d available pixels w/ %d per set" + , sample_size, (band->width * band->height), sample_per); + + if (inc_vals) { + values = rtalloc(sizeof(double) * sample_size); + if (NULL == values) { + rtwarn("Unable to allocate memory for values"); + inc_vals = 0; + } + } + + for (x = 0, j = 0, k = 0; x < outer; x++) { + y = -1; + diff = 0; + + for (i = 0, z = 0; i < sample_per; i++) { + if (sample == 1) + y = i; + else { + offset = (rand() % sample_int) + 1; + y += diff + offset; + diff = sample_int - offset; + } + RASTER_DEBUGF(5, "(x, y, z) = (%d, %d, %d)", (byY ? y : x), (byY ? x : y), z); + if (y >= inner || z > sample_per) break; + + if (byY) + rtn = rt_band_get_pixel(band, y, x, &value); + else + rtn = rt_band_get_pixel(band, x, y, &value); + + j++; + if (rtn != -1) { + if ( + hasnodata || + (!hasnodata && (hasnodata_flag != FALSE) && (value != nodata)) + ) { + + /* inc_vals set, collect pixel values */ + if (inc_vals) { + values[k] = value; + } + + /* average */ + k++; + sum += value; + + /* + one-pass standard deviation + http://www.eecs.berkeley.edu/~mhoemmen/cs194/Tutorials/variance.pdf + */ + if (k == 1) { + Q = 0; + M = value; + } + else { + Q = Q + (((k - 1) * pow(value - M, 2)) / k); + M = M + ((value - M ) / k); + } + + /* min/max */ + if (NULL == stats) { + stats = (rt_bandstats) rtalloc(sizeof(struct rt_bandstats_t)); + if (NULL == stats) { + rterror("rt_band_get_summary_stats: Unable to allocate memory for stats"); + return NULL; + } + + stats->sample = sample; + stats->count = 0; + + stats->mean = 0; + stats->stddev = -1; + stats->min = stats->max = value; + + stats->values = NULL; + stats->sorted = 0; + + } + else { + if (value < stats->min) + stats->min = value; + if (value > stats->max) + stats->max = value; + } + + } + } + + z++; + } + } + + RASTER_DEBUG(3, "sampling complete"); + + if (k > 0) { + if (inc_vals) { + /* free unused memory */ + if (sample_size != k) { + rtrealloc(values, k * sizeof(double)); + } + + stats->values = values; + } + + stats->count = k; + stats->mean = sum / k; + + /* standard deviation */ + if (sample == 1) + stats->stddev = sqrt(Q / k); + /* sample deviation */ + else { + if (k < 2) + stats->stddev = -1; + else + stats->stddev = sqrt(Q / (k - 1)); + } + } + /* inc_vals thus values allocated but not used */ + else if (inc_vals) { + rtdealloc(values); + } + +#if POSTGIS_DEBUG_LEVEL > 0 + if (NULL != stats) { + stop = clock(); + elapsed = ((double) (stop - start)) / CLOCKS_PER_SEC; + RASTER_DEBUGF(3, "(time, count, mean, stddev, min, max) = (%0.4f, %d, %f, %f, %f, %f)", + elapsed, stats->count, stats->mean, stats->stddev, stats->min, stats->max); + } +#endif + + RASTER_DEBUG(3, "done"); + return stats; +} + /*- rt_raster --------------------------------------------------------*/ struct rt_raster_serialized_t { diff --git a/raster/rt_core/rt_api.h b/raster/rt_core/rt_api.h index 9308b999e..968294bca 100644 --- a/raster/rt_core/rt_api.h +++ b/raster/rt_core/rt_api.h @@ -82,7 +82,7 @@ #endif #endif -#include /* For size_t */ +#include /* For size_t, srand and rand */ #include /* For C99 int types */ #include "liblwgeom.h" @@ -443,6 +443,20 @@ int rt_band_is_nodata(rt_band band); int rt_band_check_is_nodata(rt_band band); +/** + * Compute summary statistics for a band + * + * @param band: the band to query for minimum and maximum pixel values + * @param hasnodata: if zero, ignore nodata value + * @param sample: percentage of pixels to sample + * @param inc_vals: flag to include values in return struct + * + * @return the summary statistics for a band + */ +typedef struct rt_bandstats_t* rt_bandstats; +rt_bandstats rt_band_get_summary_stats(rt_band band, int hasnodata, + double sample, int inc_vals); + /*- rt_raster --------------------------------------------------------*/ diff --git a/raster/rt_pg/rt_pg.c b/raster/rt_pg/rt_pg.c index 8eeea487a..4b169c8ed 100644 --- a/raster/rt_pg/rt_pg.c +++ b/raster/rt_pg/rt_pg.c @@ -187,6 +187,9 @@ Datum RASTER_mapAlgebra(PG_FUNCTION_ARGS); /* create new raster from existing raster's bands */ Datum RASTER_band(PG_FUNCTION_ARGS); +/* Get summary stats */ +Datum RASTER_summaryStats(PG_FUNCTION_ARGS); + /* Replace function taken from * http://ubuntuforums.org/showthread.php?s=aa6f015109fd7e4c7e30d2fd8b717497&t=141670&page=3 @@ -2745,6 +2748,219 @@ Datum RASTER_band(PG_FUNCTION_ARGS) PG_RETURN_POINTER(pgrast); } +struct rt_bandstats_t { + double sample; + uint32_t count; + + double min; + double max; + double mean; + double stddev; + + double *values; + int sorted; +}; + +/** + * Get summary stats of a band + */ +PG_FUNCTION_INFO_V1(RASTER_summaryStats); +Datum RASTER_summaryStats(PG_FUNCTION_ARGS) +{ + FuncCallContext *funcctx; + TupleDesc tupdesc; + AttInMetadata *attinmeta; + + rt_bandstats stats; + rt_bandstats stats2; + 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; + bool hasnodata = FALSE; + int num_bands = 0; + double sample = 0; + + /* 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_summaryStats: 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); + PG_RETURN_NULL(); + } + 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); + PG_RETURN_NULL(); + } + else if (sample == 0) + sample = 1; + } + else + sample = 1; + + /* 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); + PG_RETURN_NULL(); + } + + /* we don't need the raw values, hence the zero parameter */ + stats = rt_band_get_summary_stats(band, (int) hasnodata, sample, 0); + rt_band_destroy(band); + rt_raster_destroy(raster); + if (NULL == stats) { + elog(NOTICE, "Could not retrieve summary statistics of band of index %d. Returning NULL", bandindex); + PG_RETURN_NULL(); + } + + /* store needed data */ + funcctx->user_fctx = stats; + + /* total number of tuples to return */ + funcctx->max_calls = 1; + + /* 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; + stats2 = 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, "tuple #%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(5 * sizeof(char *)); + + values[0] = (char *) palloc( + (snprintf(NULL, 0, "%d", stats2[call_cntr].count) + 1) * sizeof(char) + ); + values[1] = (char *) palloc( + (snprintf(NULL, 0, "%f", stats2[call_cntr].mean) + 1) * sizeof(char) + ); + values[2] = (char *) palloc( + (snprintf(NULL, 0, "%f", stats2[call_cntr].stddev) + 1) * sizeof(char) + ); + values[3] = (char *) palloc( + (snprintf(NULL, 0, "%f", stats2[call_cntr].min) + 1) * sizeof(char) + ); + values[4] = (char *) palloc( + (snprintf(NULL, 0, "%f", stats2[call_cntr].max) + 1) * sizeof(char) + ); + + snprintf( + values[0], + (snprintf(NULL, 0, "%d", stats2[call_cntr].count) + 1) * sizeof(char), + "%d", + stats2[call_cntr].count + ); + snprintf( + values[1], + (snprintf(NULL, 0, "%f", stats2[call_cntr].mean) + 1) * sizeof(char), + "%f", + stats2[call_cntr].mean + ); + snprintf( + values[2], + (snprintf(NULL, 0, "%f", stats2[call_cntr].stddev) + 1) * sizeof(char), + "%f", + stats2[call_cntr].stddev + ); + snprintf( + values[3], + (snprintf(NULL, 0, "%f", stats2[call_cntr].min) + 1) * sizeof(char), + "%f", + stats2[call_cntr].min + ); + snprintf( + values[4], + (snprintf(NULL, 0, "%f", stats2[call_cntr].max) + 1) * sizeof(char), + "%f", + stats2[call_cntr].max + ); + + /* build a tuple */ + tuple = BuildTupleFromCStrings(attinmeta, values); + + /* make the tuple into a datum */ + result = HeapTupleGetDatum(tuple); + + /* clean up */ + pfree(values[1]); + pfree(values[0]); + pfree(values); + + SRF_RETURN_NEXT(funcctx, result); + } + else { + pfree(stats2); + 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 5b8fd3210..11937523b 100644 --- a/raster/rt_pg/rtpostgis.sql.in.c +++ b/raster/rt_pg/rtpostgis.sql.in.c @@ -256,6 +256,199 @@ CREATE OR REPLACE FUNCTION st_band(rast raster) AS $$ SELECT st_band($1, ARRAY[1]) $$ LANGUAGE 'SQL' IMMUTABLE STRICT; +----------------------------------------------------------------------- +-- ST_SummaryStats and ST_ApproxSummaryStats +----------------------------------------------------------------------- +CREATE TYPE summarystats AS ( + count integer, + mean double precision, + stddev double precision, + min double precision, + max double precision +); + +CREATE OR REPLACE FUNCTION _st_summarystats(rast raster, nband int, hasnodata boolean, sample_percent double precision) + RETURNS summarystats + AS 'MODULE_PATHNAME','RASTER_summaryStats' + LANGUAGE 'C' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_summarystats(rast raster, nband int, hasnodata boolean) + RETURNS summarystats + AS $$ SELECT count, mean, stddev, min, max FROM _st_summarystats($1, $2, $3, 1) $$ + LANGUAGE 'SQL' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_summarystats(rast raster, nband int) + RETURNS summarystats + AS $$ SELECT count, mean, stddev, min, max FROM _st_summarystats($1, $2, FALSE, 1) $$ + LANGUAGE 'SQL' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_summarystats(rast raster, hasnodata boolean) + RETURNS summarystats + AS $$ SELECT count, mean, stddev, min, max FROM _st_summarystats($1, 1, $2, 1) $$ + LANGUAGE 'SQL' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_summarystats(rast raster) + RETURNS summarystats + AS $$ SELECT count, mean, stddev, min, max FROM _st_summarystats($1, 1, FALSE, 1) $$ + LANGUAGE 'SQL' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxsummarystats(rast raster, nband int, hasnodata boolean, sample_percent double precision) + RETURNS summarystats + AS $$ SELECT count, mean, stddev, min, max FROM _st_summarystats($1, $2, $3, $4) $$ + LANGUAGE 'SQL' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxsummarystats(rast raster, nband int, sample_percent double precision) + RETURNS summarystats + AS $$ SELECT count, mean, stddev, min, max FROM _st_summarystats($1, $2, FALSE, $3) $$ + LANGUAGE 'SQL' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxsummarystats(rast raster, hasnodata boolean, sample_percent double precision) + RETURNS summarystats + AS $$ SELECT count, mean, stddev, min, max FROM _st_summarystats($1, 1, $2, $3) $$ + LANGUAGE 'SQL' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxsummarystats(rast raster, sample_percent double precision) + RETURNS summarystats + AS $$ SELECT count, mean, stddev, min, max FROM _st_summarystats($1, 1, FALSE, $2) $$ + LANGUAGE 'SQL' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxsummarystats(rast raster) + RETURNS summarystats + AS $$ SELECT count, mean, stddev, min, max FROM _st_summarystats($1, 1, FALSE, 0.1) $$ + LANGUAGE 'SQL' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION _st_summarystats(rastertable text, rastercolumn text, nband integer, hasnodata boolean, sample_percent double precision) + RETURNS summarystats + AS $$ + DECLARE + curs refcursor; + + ctable text; + ccolumn text; + rast raster; + stats summarystats; + rtn summarystats; + + mean double precision; + tmp double precision; + BEGIN + -- nband + IF nband < 1 THEN + RAISE WARNING 'Invalid band index (must use 1-based). Returning NULL'; + RETURN NULL; + END IF; + + -- sample percent + IF sample_percent < 0 OR sample_percent > 1 THEN + RAISE WARNING 'Invalid sample percentage (must be between 0 and 1). Returning NULL'; + RETURN NULL; + END IF; + + -- clean rastertable and rastercolumn + ctable := quote_ident(rastertable); + ccolumn := quote_ident(rastercolumn); + + BEGIN + OPEN curs FOR EXECUTE 'SELECT ' + || ccolumn + || ' FROM ' + || ctable + || ' WHERE ' + || ccolumn + || ' IS NOT NULL'; + EXCEPTION + WHEN OTHERS THEN + RAISE WARNING 'Invalid table or column name. Returning NULL'; + RETURN NULL; + END; + + mean := 0; + tmp := 0; + LOOP + FETCH curs INTO rast; + EXIT WHEN NOT FOUND; + + SELECT * FROM _st_summarystats(rast, nband, hasnodata, sample_percent) INTO stats; + IF stats IS NULL OR stats.count < 1 THEN + CONTINUE; + END IF; + --RAISE NOTICE 'stats = %', stats; + + mean := mean + (stats.count * stats.mean); -- not the mean + -- sum of "sum of squares" + IF sample_percent > 0 AND sample_percent < 1 THEN + tmp := tmp + (power(stats.stddev, 2) * (stats.count - 1)); + ELSE + tmp := tmp + (power(stats.stddev, 2) * stats.count); + END IF; + + IF rtn.count IS NULL THEN + rtn := stats; + ELSE + rtn.count = rtn.count + stats.count; + + IF stats.min < rtn.min THEN + rtn.min := stats.min; + END IF; + IF stats.max > rtn.max THEN + rtn.max := stats.max; + END IF; + END IF; + + END LOOP; + + CLOSE curs; + + -- weighted mean and cumulative stddev + IF rtn.count IS NOT NULL AND rtn.count > 0 THEN + rtn.mean := mean / rtn.count; + rtn.stddev := sqrt(tmp / rtn.count); + END IF; + + RETURN rtn; + END; + $$ LANGUAGE 'plpgsql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_summarystats(rastertable text, rastercolumn text, nband integer, hasnodata boolean) + RETURNS summarystats + AS $$ SELECT count, mean, stddev, min, max FROM _st_summarystats($1, $2, $3, $4, 1) $$ + LANGUAGE 'SQL' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_summarystats(rastertable text, rastercolumn text, nband integer) + RETURNS summarystats + AS $$ SELECT count, mean, stddev, min, max FROM _st_summarystats($1, $2, $3, FALSE, 1) $$ + LANGUAGE 'SQL' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_summarystats(rastertable text, rastercolumn text) + RETURNS summarystats + AS $$ SELECT count, mean, stddev, min, max FROM _st_summarystats($1, $2, 1, FALSE, 1) $$ + LANGUAGE 'SQL' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxsummarystats(rastertable text, rastercolumn text, nband integer, hasnodata boolean, sample_percent double precision) + RETURNS summarystats + AS $$ SELECT count, mean, stddev, min, max FROM _st_summarystats($1, $2, $3, $4, $5) $$ + LANGUAGE 'SQL' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxsummarystats(rastertable text, rastercolumn text, nband integer, sample_percent double precision) + RETURNS summarystats + AS $$ SELECT count, mean, stddev, min, max FROM _st_summarystats($1, $2, $3, FALSE, $4) $$ + LANGUAGE 'SQL' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxsummarystats(rastertable text, rastercolumn text, nband integer) + RETURNS summarystats + AS $$ SELECT count, mean, stddev, min, max FROM _st_summarystats($1, $2, $3, FALSE, 0.1) $$ + LANGUAGE 'SQL' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxsummarystats(rastertable text, rastercolumn text, sample_percent double precision) + RETURNS summarystats + AS $$ SELECT count, mean, stddev, min, max FROM _st_summarystats($1, $2, 1, FALSE, $3) $$ + LANGUAGE 'SQL' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxsummarystats(rastertable text, rastercolumn text) + RETURNS summarystats + AS $$ SELECT count, mean, stddev, min, max FROM _st_summarystats($1, $2, 1, FALSE, 0.1) $$ + LANGUAGE 'SQL' IMMUTABLE STRICT; + ----------------------------------------------------------------------- -- MapAlgebra ----------------------------------------------------------------------- diff --git a/raster/test/core/testapi.c b/raster/test/core/testapi.c index 352dceb52..7e2af187a 100644 --- a/raster/test/core/testapi.c +++ b/raster/test/core/testapi.c @@ -987,6 +987,90 @@ static void testRasterFromBand(rt_raster raster) { rt_raster_destroy(rast); } +struct rt_bandstats_t { + double sample; + uint32_t count; + + double min; + double max; + double mean; + double stddev; + + double *values; + int sorted; +}; +static void testBandStats() { + rt_bandstats stats = NULL; + + rt_raster raster; + rt_band band; + uint32_t x; + uint32_t xmax = 10000; + uint32_t y; + uint32_t ymax = 10000; + double nodata; + int rtn; + + raster = rt_raster_new(xmax, ymax); + assert(raster); /* or we're out of virtual memory */ + band = addBand(raster, PT_32BUI, 0, 0); + CHECK(band); + rt_band_set_nodata(band, 0); + + for (x = 0; x < xmax; x++) { + for (y = 0; y < ymax; y++) { + rtn = rt_band_set_pixel(band, x, y, x + y); + CHECK((rtn != -1)); + } + } + + nodata = rt_band_get_nodata(band); + CHECK_EQUALS(nodata, 0); + + stats = (rt_bandstats) rt_band_get_summary_stats(band, 0, 0, 1); + CHECK(stats); + CHECK_EQUALS(stats->min, 1); + CHECK_EQUALS(stats->max, 19998); + + rtdealloc(stats->values); + rtdealloc(stats); + + stats = (rt_bandstats) rt_band_get_summary_stats(band, 0, 0.1, 1); + CHECK(stats); + + rtdealloc(stats->values); + rtdealloc(stats); + + stats = (rt_bandstats) rt_band_get_summary_stats(band, 0, 0.15, 0); + CHECK(stats); + rtdealloc(stats); + + stats = (rt_bandstats) rt_band_get_summary_stats(band, 0, 0.2, 0); + CHECK(stats); + rtdealloc(stats); + + stats = (rt_bandstats) rt_band_get_summary_stats(band, 0, 0.25, 0); + CHECK(stats); + rtdealloc(stats); + + stats = (rt_bandstats) rt_band_get_summary_stats(band, 1, 0, 1); + CHECK(stats); + CHECK_EQUALS(stats->min, 0); + CHECK_EQUALS(stats->max, 19998); + + rtdealloc(stats->values); + rtdealloc(stats); + + stats = (rt_bandstats) rt_band_get_summary_stats(band, 1, 0.1, 1); + CHECK(stats); + + rtdealloc(stats->values); + rtdealloc(stats); + + rt_band_destroy(band); + rt_raster_destroy(raster); +} + int main() { @@ -1297,6 +1381,10 @@ main() testRasterFromBand(raster); printf("Successfully tested rt_raster_from_band\n"); + printf("Testing band stats\n"); + testBandStats(); + printf("Successfully tested band stats\n"); + deepRelease(raster); diff --git a/raster/test/regress/Makefile.in b/raster/test/regress/Makefile.in index 3cb356767..96c8679b8 100644 --- a/raster/test/regress/Makefile.in +++ b/raster/test/regress/Makefile.in @@ -61,6 +61,7 @@ TEST_BANDPROPS = \ create_rt_band_properties_test.sql \ rt_band_properties.sql \ rt_set_band_properties.sql \ + rt_summarystats.sql \ $(NULL) TEST_PIXEL = \ diff --git a/raster/test/regress/rt_summarystats.sql b/raster/test/regress/rt_summarystats.sql new file mode 100644 index 000000000..9b09b10ed --- /dev/null +++ b/raster/test/regress/rt_summarystats.sql @@ -0,0 +1,125 @@ +SELECT * FROM ST_SummaryStats( + 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 + ) + , FALSE +); +SELECT ST_SummaryStats( + 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 + ) + , FALSE +); +SELECT count FROM ST_SummaryStats( + 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 + ) + , FALSE +); +SELECT count FROM ST_SummaryStats( + 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 +); +SELECT mean, stddev FROM ST_SummaryStats( + 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 + ) + , FALSE +); +SELECT mean, stddev FROM ST_SummaryStats( + 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 +); +BEGIN; +CREATE TEMP TABLE test + ON COMMIT DROP AS + SELECT + rast.rast + FROM ( + SELECT 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 + ) AS rast + ) AS rast + FULL JOIN ( + SELECT generate_series(1, 10) AS id + ) AS id + ON 1 = 1; +SELECT * FROM ST_SummaryStats('test', 'rast', 1, FALSE); +SELECT * FROM ST_SummaryStats('test', 'rast', 1, TRUE); +SELECT * FROM ST_SummaryStats('test', 'rast', 1); +SELECT * FROM ST_SummaryStats('test', 'rast'); +ROLLBACK; diff --git a/raster/test/regress/rt_summarystats_expected b/raster/test/regress/rt_summarystats_expected new file mode 100644 index 000000000..3fdf85f33 --- /dev/null +++ b/raster/test/regress/rt_summarystats_expected @@ -0,0 +1,13 @@ +2|-3.429205|6.570795|-10|3.14159 +(2,-3.429205,6.570795,-10,3.14159) +2 +100 +-3.429205|6.570795 +-0.068584|1.045941 +BEGIN +SELECT 10 +20|-3.429205|6.570795|-10|3.14159 +1000|-0.068584|1.045941|-10|3.14159 +20|-3.429205|6.570795|-10|3.14159 +20|-3.429205|6.570795|-10|3.14159 +COMMIT