#include <assert.h>
#include <float.h> /* for FLT_EPSILON and float type limits */
#include <limits.h> /* for integer type limits */
+#include <time.h> /* for time */
#include "rt_api.h"
#define POSTGIS_RASTER_WARN_ON_TRUNCATION
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 {
#endif
#endif
-#include <stdlib.h> /* For size_t */
+#include <stdlib.h> /* For size_t, srand and rand */
#include <stdint.h> /* For C99 int types */
#include "liblwgeom.h"
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 --------------------------------------------------------*/
/* 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
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 */
/* ---------------------------------------------------------------- */
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
-----------------------------------------------------------------------
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()
{
testRasterFromBand(raster);
printf("Successfully tested rt_raster_from_band\n");
+ printf("Testing band stats\n");
+ testBandStats();
+ printf("Successfully tested band stats\n");
+
deepRelease(raster);
create_rt_band_properties_test.sql \
rt_band_properties.sql \
rt_set_band_properties.sql \
+ rt_summarystats.sql \
$(NULL)
TEST_PIXEL = \
--- /dev/null
+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;
--- /dev/null
+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
projects / postgis / commitdiff