}
/* quicksort */
-#define swap(x, y) { double t; t = x; x = y; y = t; }
-#define order(x, y) if (x > y) swap(x, y)
+#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;
r = *right;
/* order */
- order(l, m);
- order(l, r);
- order(m, r);
+ ORDER(l, m);
+ ORDER(l, r);
+ ORDER(m, r);
/* pivot is higher of two values */
if (l < m) return m;
while (*right >= pivot) --right;
if (left < right) {
- swap(*left, *right);
+ SWAP(*left, *right);
++left;
--right;
}
/**
* Compute summary statistics for a band
*
- * @param band: the band to query for minimum and maximum pixel values
+ * @param band: the band to query for summary stats
* @param hasnodata: if non-zero, ignore nodata values
* @param sample: percentage of pixels to sample
* @param inc_vals: flag to include values in return struct
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;
}
data = rt_band_get_data(band);
- pixtype = band->pixtype;
hasnodata_flag = rt_band_get_hasnodata_flag(band);
- if (hasnodata_flag != FALSE) {
+ if (hasnodata_flag != FALSE)
nodata = rt_band_get_nodata(band);
- }
- else {
+ else
hasnodata = 0;
- }
RASTER_DEBUGF(3, "nodata = %f", nodata);
RASTER_DEBUGF(3, "hasnodata_flag = %d", hasnodata_flag);
return rtn;
}
+/* symmetrical rounding */
+#define ROUND(x, y) (((x > 0.0) ? floor((x * pow(10, y) + 0.5)) : ceil((x * pow(10, y) - 0.5))) / pow(10, y));
+
+struct rt_valuecount_t {
+ double value;
+ uint32_t count;
+ double proportion;
+};
+
+/**
+ * Count the number of times provided value(s) occur in
+ * the band
+ *
+ * @param band: the band to query for minimum and maximum pixel values
+ * @param hasnodata: if non-zero, ignore nodata values
+ * @param search_values: array of values to count
+ * @param search_values_count: the number of search values
+ * @param roundto: the decimal place to round the values to
+ * @param rtn_count: the number of value counts being returned
+ *
+ * @return the default set of or requested quantiles for a band
+ */
+rt_valuecount
+rt_band_get_value_count(rt_band band, int hasnodata,
+ double *search_values, uint32_t search_values_count, double roundto,
+ int *rtn_count) {
+ rt_valuecount vcnts = NULL;
+ rt_pixtype pixtype = PT_END;
+ uint8_t *data = NULL;
+ int hasnodata_flag = FALSE;
+ double nodata = 0;
+
+ int scale = 0;
+ int doround = 0;
+ double tmpd = 0;
+ int i = 0;
+
+ uint32_t x = 0;
+ uint32_t y = 0;
+ int rtn;
+ double pxlval;
+ double rpxlval;
+ uint32_t total = 0;
+ int vcnts_count = 0;
+ int new_valuecount = 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_count_of_values 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 = 0;
+
+ RASTER_DEBUGF(3, "nodata = %f", nodata);
+ RASTER_DEBUGF(3, "hasnodata_flag = %d", hasnodata_flag);
+ RASTER_DEBUGF(3, "user hasnodata = %d", hasnodata);
+
+ /* process roundto */
+ if (roundto < 0 || (fabs(roundto - 0.0) < FLT_EPSILON)) {
+ roundto = 0;
+ scale = 0;
+ }
+ /* tenths, hundredths, thousandths, etc */
+ else if (roundto < 1) {
+ switch (pixtype) {
+ /* integer band types don't have digits after the decimal place */
+ case PT_1BB:
+ case PT_2BUI:
+ case PT_4BUI:
+ case PT_8BSI:
+ case PT_8BUI:
+ case PT_16BSI:
+ case PT_16BUI:
+ case PT_32BSI:
+ case PT_32BUI:
+ roundto = 0;
+ break;
+ /* floating points, check the rounding */
+ case PT_32BF:
+ case PT_64BF:
+ for (scale = 0; scale <= 20; scale++) {
+ tmpd = roundto * pow(10, scale);
+ if (fabs((tmpd - ((int) tmpd)) - 0.0) < FLT_EPSILON) break;
+ }
+ break;
+ case PT_END:
+ break;
+ }
+ }
+ /* ones, tens, hundreds, etc */
+ else {
+ for (scale = 0; scale >= -20; scale--) {
+ tmpd = roundto * pow(10, scale);
+ if (tmpd < 1 || fabs(tmpd - 1.0) < FLT_EPSILON) {
+ if (scale == 0) doround = 1;
+ break;
+ }
+ }
+ }
+
+ if (scale != 0 || doround)
+ doround = 1;
+ else
+ doround = 0;
+ RASTER_DEBUGF(3, "scale = %d", scale);
+ RASTER_DEBUGF(3, "doround = %d", doround);
+
+ /* process search_values */
+ if (search_values_count > 0 && NULL != search_values) {
+ vcnts = (rt_valuecount) rtalloc(sizeof(struct rt_valuecount_t) * search_values_count);
+ if (NULL == vcnts) {
+ rterror("rt_band_get_count_of_values: Unable to allocate memory for value counts");
+ *rtn_count = 0;
+ return NULL;
+ }
+
+ for (i = 0; i < search_values_count; i++) {
+ vcnts[i].count = 0;
+ vcnts[i].proportion = 0;
+ if (!doround)
+ vcnts[i].value = search_values[i];
+ else
+ vcnts[i].value = ROUND(search_values[i], scale);
+ }
+ vcnts_count = i;
+ }
+ else
+ search_values_count = 0;
+ RASTER_DEBUGF(3, "search_values_count = %d", search_values_count);
+
+ /* 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 {
+ if (search_values_count > 0) {
+ /* check for nodata match */
+ for (i = 0; i < search_values_count; i++) {
+ if (!doround)
+ tmpd = nodata;
+ else
+ tmpd = ROUND(nodata, scale);
+
+ if (fabs(tmpd - vcnts[i].value) > FLT_EPSILON)
+ continue;
+
+ vcnts[i].count = band->width * band->height;
+ vcnts->proportion = 1.0;
+ }
+
+ *rtn_count = vcnts_count;
+ }
+ /* no defined search values */
+ else {
+ vcnts = (rt_valuecount) rtalloc(sizeof(struct rt_valuecount_t));
+ if (NULL == vcnts) {
+ rterror("rt_band_get_count_of_values: Unable to allocate memory for value counts");
+ *rtn_count = 0;
+ return NULL;
+ }
+
+ vcnts->value = nodata;
+ vcnts->count = band->width * band->height;
+ vcnts->proportion = 1.0;
+
+ *rtn_count = 1;
+ }
+
+ return vcnts;
+ }
+ }
+
+ for (x = 0; x < band->width; x++) {
+ for (y = 0; y < band->height; y++) {
+ rtn = rt_band_get_pixel(band, x, y, &pxlval);
+
+ /* error getting value, continue */
+ if (rtn == -1) continue;
+
+ if (
+ !hasnodata || (
+ hasnodata &&
+ (hasnodata_flag != FALSE) &&
+ (fabs(pxlval - nodata) > FLT_EPSILON)
+ )
+ ) {
+ total++;
+ if (doround) {
+ rpxlval = ROUND(pxlval, scale);
+ }
+ else
+ rpxlval = pxlval;
+ RASTER_DEBUGF(5, "(pxlval, rpxlval) => (%0.6f, %0.6f)", pxlval, rpxlval);
+
+ new_valuecount = 1;
+ /* search for match in existing valuecounts */
+ for (i = 0; i < vcnts_count; i++) {
+ /* match found */
+ if (fabs(vcnts[i].value - rpxlval) < FLT_EPSILON) {
+ vcnts[i].count++;
+ new_valuecount = 0;
+ RASTER_DEBUGF(5, "(value, count) => (%0.6f, %d)", vcnts[i].value, vcnts[i].count);
+ break;
+ }
+ }
+
+ /*
+ don't add new valuecount either because
+ - no need for new one
+ - user-defined search values
+ */
+ if (!new_valuecount || search_values_count > 0) continue;
+
+ /* add new valuecount */
+ vcnts = rtrealloc(vcnts, sizeof(struct rt_valuecount_t) * (vcnts_count + 1));
+ if (NULL == vcnts) {
+ rterror("rt_band_get_count_of_values: Unable to allocate memory for value counts");
+ *rtn_count = 0;
+ return NULL;
+ }
+
+ vcnts[vcnts_count].value = rpxlval;
+ vcnts[vcnts_count].count = 1;
+ vcnts[vcnts_count].proportion = 0;
+ RASTER_DEBUGF(5, "(value, count) => (%0.6f, %d)", vcnts[vcnts_count].value, vcnts[vcnts_count].count);
+ vcnts_count++;
+ }
+ }
+ }
+
+#if POSTGIS_DEBUG_LEVEL > 0
+ stop = clock();
+ elapsed = ((double) (stop - start)) / CLOCKS_PER_SEC;
+ RASTER_DEBUGF(3, "elapsed time = %0.4f", elapsed);
+#endif
+
+ for (i = 0; i < vcnts_count; i++) {
+ vcnts[i].proportion = (double) vcnts[i].count / total;
+ RASTER_DEBUGF(5, "(value, count) => (%0.6f, %d)", vcnts[i].value, vcnts[i].count);
+ }
+
+ RASTER_DEBUG(3, "done");
+ *rtn_count = vcnts_count;
+ return vcnts;
+}
+
struct rt_reclassexpr_t {
struct rt_reclassrange {
double min;
uint32_t hasnodata, double nodataval, rt_reclassexpr *exprset,
int exprcount) {
rt_band band = NULL;
- int width = 0;
- int height = 0;
+ uint32_t width = 0;
+ uint32_t height = 0;
int numval = 0;
int memsize = 0;
void *mem = NULL;
double src_nodataval = 0.0;
int rtn;
- int x;
- int y;
+ uint32_t x;
+ uint32_t y;
int i;
double or = 0;
double ov = 0;
/**
* Compute summary statistics for a band
*
- * @param band: the band to query for minimum and maximum pixel values
+ * @param band: the band to query for summary stats
* @param hasnodata: if non-zero, ignore nodata values
* @param sample: percentage of pixels to sample
* @param inc_vals: flag to include values in return struct
rt_quantile rt_band_get_quantiles(rt_bandstats stats,
double *quantiles, int quantiles_count, int *rtn_count);
+/**
+ * Count the number of times provided value(s) occur in
+ * the band
+ *
+ * @param band: the band to query for minimum and maximum pixel values
+ * @param hasnodata: if non-zero, ignore nodata values
+ * @param search_values: array of values to count
+ * @param search_values_count: the number of search values
+ * @param roundto: the decimal place to round the values to
+ * @param rtn_count: the number of value counts being returned
+ *
+ * @return the default set of or requested quantiles for a band
+ */
+typedef struct rt_valuecount_t* rt_valuecount;
+rt_valuecount rt_band_get_value_count(rt_band band, int hasnodata,
+ double *search_values, uint32_t search_values_count,
+ double roundto, int *rtn_count);
+
/**
* Returns new band with values reclassified
*
/* get quantiles */
Datum RASTER_quantile(PG_FUNCTION_ARGS);
+/* get counts of values */
+Datum RASTER_valueCount(PG_FUNCTION_ARGS);
+
/* reclassify specified bands of a raster */
Datum RASTER_reclass(PG_FUNCTION_ARGS);
}
}
+struct rt_valuecount_t {
+ double value;
+ uint32_t count;
+ double proportion;
+};
+
+/* get counts of values */
+PG_FUNCTION_INFO_V1(RASTER_valueCount);
+Datum RASTER_valueCount(PG_FUNCTION_ARGS) {
+ FuncCallContext *funcctx;
+ TupleDesc tupdesc;
+ AttInMetadata *attinmeta;
+
+ int count;
+ rt_valuecount vcnts;
+ rt_valuecount vcnts2;
+ int call_cntr;
+ int max_calls;
+
+ /* first call of function */
+ if (SRF_IS_FIRSTCALL()) {
+ MemoryContext oldcontext;
+
+ rt_pgraster *pgraster = NULL;
+ rt_raster raster = NULL;
+ rt_band band = NULL;
+ int32_t bandindex = 0;
+ int num_bands = 0;
+ bool hasnodata = TRUE;
+ double *search_values = NULL;
+ int search_values_count = 0;
+ double roundto = 0;
+
+ 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);
+
+ /* pgraster is null, return nothing */
+ if (PG_ARGISNULL(0)) SRF_RETURN_DONE(funcctx);
+ pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+
+ raster = rt_raster_deserialize(pgraster);
+ if (!raster) {
+ elog(ERROR, "RASTER_valueCount: 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);
+
+ /* search values */
+ if (!PG_ARGISNULL(3)) {
+ array = PG_GETARG_ARRAYTYPE_P(3);
+ etype = ARR_ELEMTYPE(array);
+ get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
+
+ switch (etype) {
+ case FLOAT4OID:
+ case FLOAT8OID:
+ break;
+ default:
+ elog(ERROR, "RASTER_valueCount: Invalid data type for values");
+ 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);
+
+ search_values = palloc(sizeof(double) * n);
+ for (i = 0, j = 0; i < n; i++) {
+ if (nulls[i]) continue;
+
+ switch (etype) {
+ case FLOAT4OID:
+ search_values[j] = (double) DatumGetFloat4(e[i]);
+ break;
+ case FLOAT8OID:
+ search_values[j] = (double) DatumGetFloat8(e[i]);
+ break;
+ }
+
+ POSTGIS_RT_DEBUGF(5, "search_values[%d] = %f", j, search_values[j]);
+ j++;
+ }
+ search_values_count = j;
+
+ if (j < 1) {
+ pfree(search_values);
+ search_values = NULL;
+ }
+ }
+
+ /* roundto */
+ if (!PG_ARGISNULL(4)) {
+ roundto = PG_GETARG_FLOAT8(4);
+ if (roundto < 0.) roundto = 0;
+ }
+
+ /* 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 counts of values */
+ vcnts = rt_band_get_value_count(band, (int) hasnodata, search_values, search_values_count, roundto, &count);
+ rt_band_destroy(band);
+ rt_raster_destroy(raster);
+ if (NULL == vcnts || !count) {
+ elog(NOTICE, "Could not count the values of raster band of index %d", bandindex);
+ SRF_RETURN_DONE(funcctx);
+ }
+
+ POSTGIS_RT_DEBUGF(3, "%d value counts returned", count);
+
+ /* Store needed information */
+ funcctx->user_fctx = vcnts;
+
+ /* 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;
+ vcnts2 = 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(3 * sizeof(char *));
+
+ values[0] = (char *) palloc(sizeof(char) * (MAX_DBL_CHARLEN + 1));
+ values[1] = (char *) palloc(sizeof(char) * (MAX_INT_CHARLEN + 1));
+ values[2] = (char *) palloc(sizeof(char) * (MAX_DBL_CHARLEN + 1));
+
+ snprintf(
+ values[0],
+ sizeof(char) * (MAX_DBL_CHARLEN + 1),
+ "%f",
+ vcnts2[call_cntr].value
+ );
+ snprintf(
+ values[1],
+ sizeof(char) * (MAX_INT_CHARLEN + 1),
+ "%d",
+ vcnts2[call_cntr].count
+ );
+ snprintf(
+ values[2],
+ sizeof(char) * (MAX_DBL_CHARLEN + 1),
+ "%f",
+ vcnts2[call_cntr].proportion
+ );
+
+ /* 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[2]);
+ pfree(values[1]);
+ pfree(values[0]);
+ pfree(values);
+
+ SRF_RETURN_NEXT(funcctx, result);
+ }
+ /* do when there is no more left */
+ else {
+ pfree(vcnts2);
+ SRF_RETURN_DONE(funcctx);
+ }
+}
+
struct rt_reclassexpr_t {
struct rt_reclassrange {
double min;
AS $$ SELECT quantile, value FROM _st_quantile($1, 1, $2, 1, ARRAY[$3]::double precision[]) $$
LANGUAGE 'sql' IMMUTABLE STRICT;
+-----------------------------------------------------------------------
+-- ST_ValueCount
+-----------------------------------------------------------------------
+CREATE TYPE valuecount AS (
+ value double precision,
+ count integer,
+ proportion double precision
+);
+
+-- None of the functions can be strict as "searchvalues" and "roundto" can be NULL
+
+CREATE OR REPLACE FUNCTION _st_valuecount(rast raster, nband integer, hasnodata boolean, searchvalues double precision[], roundto double precision)
+ RETURNS SETOF valuecount
+ AS 'MODULE_PATHNAME', 'RASTER_valueCount'
+ LANGUAGE 'C' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rast raster, nband integer, hasnodata boolean, searchvalues double precision[], roundto double precision, OUT value double precision, OUT count integer)
+ RETURNS SETOF record
+ AS $$ SELECT value, count FROM _st_valuecount($1, $2, $3, $4, $5) $$
+ LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rast raster, nband integer, searchvalues double precision[], roundto double precision, OUT value double precision, OUT count integer)
+ RETURNS SETOF record
+ AS $$ SELECT value, count FROM _st_valuecount($1, $2, TRUE, $3, $4) $$
+ LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rast raster, nband integer, searchvalues double precision[], OUT value double precision, OUT count integer)
+ RETURNS SETOF record
+ AS $$ SELECT value, count FROM _st_valuecount($1, $2, TRUE, $3, 0) $$
+ LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rast raster, searchvalues double precision[], roundto double precision, OUT value double precision, OUT count integer)
+ RETURNS SETOF record
+ AS $$ SELECT value, count FROM _st_valuecount($1, 1, TRUE, $2, $3) $$
+ LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rast raster, searchvalues double precision[], OUT value double precision, OUT count integer)
+ RETURNS SETOF record
+ AS $$ SELECT value, count FROM _st_valuecount($1, 1, TRUE, $2, 0) $$
+ LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rast raster, nband integer, hasnodata boolean, searchvalue double precision, roundto double precision, OUT value double precision, OUT count integer)
+ RETURNS record
+ AS $$ SELECT value, count FROM _st_valuecount($1, $2, $3, ARRAY[$4]::double precision[], $5) $$
+ LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rast raster, nband integer, searchvalue double precision, roundto double precision, OUT value double precision, OUT count integer)
+ RETURNS record
+ AS $$ SELECT value, count FROM _st_valuecount($1, $2, TRUE, ARRAY[$3]::double precision[], $4) $$
+ LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rast raster, nband integer, searchvalue double precision, OUT value double precision, OUT count integer)
+ RETURNS record
+ AS $$ SELECT value, count FROM _st_valuecount($1, $2, TRUE, ARRAY[$3]::double precision[], 0) $$
+ LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rast raster, searchvalue double precision, roundto double precision, OUT value double precision, OUT count integer)
+ RETURNS record
+ AS $$ SELECT value, count FROM _st_valuecount($1, 1, TRUE, ARRAY[$2]::double precision[], $3) $$
+ LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rast raster, searchvalue double precision, OUT value double precision, OUT count integer)
+ RETURNS record
+ AS $$ SELECT value, count FROM _st_valuecount($1, 1, TRUE, ARRAY[$2]::double precision[], 0) $$
+ LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION _st_valuecount(rastertable text, rastercolumn text, nband integer, hasnodata boolean, searchvalues double precision[], roundto double precision)
+ RETURNS SETOF valuecount
+ AS $$
+ DECLARE
+ curs refcursor;
+
+ ctable text;
+ ccolumn text;
+ rast raster;
+
+ vcnts valuecount;
+ BEGIN
+ -- nband
+ IF nband < 1 THEN
+ RAISE WARNING 'Invalid band index (must use 1-based). Returning NULL';
+ RETURN;
+ END IF;
+
+ -- rastertable and rastercolumn
+ IF rastertable IS NULL THEN
+ RAISE WARNING 'rastertable cannot be NULL. Returning NULL';
+ RETURN;
+ ELSEIF rastercolumn IS NULL THEN
+ RAISE WARNING 'rastercolumn cannot be NULL. Returning NULL';
+ RETURN;
+ 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;
+ END;
+
+ LOOP
+ FETCH curs INTO rast;
+ EXIT WHEN NOT FOUND;
+
+ FOR vcnts IN SELECT * FROM _st_valuecount(rast, nband, hasnodata, searchvalues, roundto) LOOP
+ IF vcnts IS NULL THEN
+ CONTINUE;
+ END IF;
+ --RAISE NOTICE 'vcnts = %', vcnts;
+
+ RETURN NEXT vcnts;
+ END LOOP;
+
+ END LOOP;
+
+ CLOSE curs;
+
+ RETURN;
+ END;
+ $$ LANGUAGE 'plpgsql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rastertable text, rastercolumn text, nband integer, hasnodata boolean, searchvalues double precision[], roundto double precision, OUT value double precision, OUT count bigint)
+ RETURNS SETOF record
+ AS $$ SELECT value, sum(count) AS count FROM _st_valuecount($1, $2, $3, $4, $5, $6) GROUP BY 1 ORDER BY 1 $$
+ LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rastertable text, rastercolumn text, nband integer, searchvalues double precision[], roundto double precision, OUT value double precision, OUT count bigint)
+ RETURNS SETOF record
+ AS $$ SELECT value, count FROM st_valuecount($1, $2, $3, TRUE, $4, $5) $$
+ LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rastertable text, rastercolumn text, nband integer, searchvalues double precision[], OUT value double precision, OUT count bigint)
+ RETURNS SETOF record
+ AS $$ SELECT value, count FROM st_valuecount($1, $2, $3, TRUE, $4, 0) $$
+ LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rastertable text, rastercolumn text, searchvalues double precision[], roundto double precision, OUT value double precision, OUT count bigint)
+ RETURNS SETOF record
+ AS $$ SELECT value, count FROM st_valuecount($1, $2, 1, TRUE, $3, $4) $$
+ LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rastertable text, rastercolumn text, searchvalues double precision[], OUT value double precision, OUT count bigint)
+ RETURNS SETOF record
+ AS $$ SELECT value, count FROM st_valuecount($1, $2, 1, TRUE, $3, 0) $$
+ LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rastertable text, rastercolumn text, nband integer, hasnodata boolean, searchvalue double precision, roundto double precision, OUT value double precision, OUT count bigint)
+ RETURNS record
+ AS $$ SELECT value, count FROM st_valuecount($1, $2, $3, $4, ARRAY[$5]::double precision[], $6) $$
+ LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rastertable text, rastercolumn text, nband integer, searchvalue double precision, roundto double precision, OUT value double precision, OUT count bigint)
+ RETURNS record
+ AS $$ SELECT value, count FROM st_valuecount($1, $2, $3, TRUE, ARRAY[$4]::double precision[], $5) $$
+ LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rastertable text, rastercolumn text, nband integer, searchvalue double precision, OUT value double precision, OUT count bigint)
+ RETURNS record
+ AS $$ SELECT value, count FROM st_valuecount($1, $2, $3, TRUE, ARRAY[$4]::double precision[], 0) $$
+ LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rastertable text, rastercolumn text, searchvalue double precision, roundto double precision, OUT value double precision, OUT count bigint)
+ RETURNS record
+ AS $$ SELECT value, count FROM st_valuecount($1, $2, 1, TRUE, ARRAY[$3]::double precision[], $4) $$
+ LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rastertable text, rastercolumn text, searchvalue double precision, OUT value double precision, OUT count bigint)
+ RETURNS record
+ AS $$ SELECT value, count FROM st_valuecount($1, $2, 1, TRUE, ARRAY[$3]::double precision[], 0) $$
+ LANGUAGE 'sql' IMMUTABLE;
+
-----------------------------------------------------------------------
-- ST_Reclass
-----------------------------------------------------------------------
rtdealloc(stats->values);
rtdealloc(stats);
- rt_band_destroy(band);
- rt_raster_destroy(raster);
+ deepRelease(raster);
}
static void testRasterReplaceBand() {
for (i = cnt - 1; i >= 0; i--)
free(exprset[i]);
free(exprset);
- rt_band_destroy(band);
rt_band_destroy(newband);
- rt_raster_destroy(raster);
+ deepRelease(raster);
}
struct rt_gdaldriver_t {
rt_raster_destroy(rast);
}
+struct rt_valuecount_t {
+ double value;
+ uint32_t count;
+ double proportion;
+};
+static void testValueCount() {
+ rt_valuecount vcnts = NULL;
+
+ rt_raster raster;
+ rt_band band;
+ uint32_t x;
+ uint32_t xmax = 100;
+ uint32_t y;
+ uint32_t ymax = 100;
+ int rtn = 0;
+
+ double count[] = {3, 4, 5};
+
+ raster = rt_raster_new(xmax, ymax);
+ assert(raster); /* or we're out of virtual memory */
+ band = addBand(raster, PT_32BF, 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, (((double) x * y) + (x + y) + (x + y * x)) / (x + y + 1));
+ CHECK((rtn != -1));
+ }
+ }
+ vcnts = rt_band_get_value_count(band, 1, NULL, 0, 0, &rtn);
+ CHECK(vcnts);
+ CHECK((rtn > 0));
+ rtdealloc(vcnts);
+
+ vcnts = rt_band_get_value_count(band, 1, NULL, 0, 0.01, &rtn);
+ CHECK(vcnts);
+ CHECK((rtn > 0));
+ rtdealloc(vcnts);
+
+ vcnts = rt_band_get_value_count(band, 1, NULL, 0, 0.1, &rtn);
+ CHECK(vcnts);
+ CHECK((rtn > 0));
+ rtdealloc(vcnts);
+
+ vcnts = rt_band_get_value_count(band, 1, NULL, 0, 1, &rtn);
+ CHECK(vcnts);
+ CHECK((rtn > 0));
+ rtdealloc(vcnts);
+
+ vcnts = rt_band_get_value_count(band, 1, NULL, 0, 10, &rtn);
+ CHECK(vcnts);
+ CHECK((rtn > 0));
+ rtdealloc(vcnts);
+
+ vcnts = rt_band_get_value_count(band, 1, count, 3, 1, &rtn);
+ CHECK(vcnts);
+ CHECK((rtn > 0));
+ rtdealloc(vcnts);
+
+ deepRelease(raster);
+}
+
int
main()
{
testGDALDrivers();
printf("Successfully tested rt_raster_gdal_drivers\n");
+ printf("Testing rt_band_get_value_count\n");
+ testValueCount();
+ printf("Successfully tested rt_band_get_value_count\n");
+
deepRelease(raster);
return EXIT_SUCCESS;
rt_minmax.sql \
rt_histogram.sql \
rt_quantile.sql \
+ rt_valuecount.sql \
$(NULL)
TEST_PIXEL = \
--- /dev/null
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+ 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[]::double precision[], 0);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+ 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, NULL::double precision[], 0);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+ 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, FALSE, NULL::double precision[], 0);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+ 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, NULL::double precision[], 0.1);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+ 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[3.1], 0.1);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+ 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[-10]);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+ 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[-10], 0);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+ 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[-10, 3]);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+ 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, 3.14, 0.19);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+ 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, FALSE, 3.14, 0.01);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+ 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, -10, 0.1);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+ 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, -10);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+ 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
+ )
+, -10., 10);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+ 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
+ )
+, 3.14159);
+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 round(value::numeric, 3), count FROM st_valuecount('test', 'rast', 1, FALSE, NULL::double precision[], 0);
+SELECT round(value::numeric, 3), count FROM st_valuecount('test', 'rast', 1, ARRAY[-10]::double precision[], 0);
+SELECT round(value::numeric, 3), count FROM st_valuecount('test', 'rast', 1, ARRAY[1]::double precision[]);
+SELECT round(value::numeric, 3), count FROM st_valuecount('test', 'rast', NULL::double precision[], 0.1);
+SELECT round(value::numeric, 3), count FROM st_valuecount('test', 'rast', ARRAY[-1, 3.1]::double precision[], 0.1);
+
+SELECT round(value::numeric, 3), count FROM st_valuecount('test', 'rast', 1, TRUE, NULL::double precision, 0);
+SELECT round(value::numeric, 3), count FROM st_valuecount('test', 'rast', 1, 3.14, 1);
+SELECT round(value::numeric, 3), count FROM st_valuecount('test', 'rast', 1, -1);
+SELECT round(value::numeric, 3), count FROM st_valuecount('test', 'rast', 3.1, 0.1);
+SELECT round(value::numeric, 3), count FROM st_valuecount('test', 'rast', -9);
+ROLLBACK;
--- /dev/null
+-10.000|1
+3.142|1
+-10.000|1
+3.142|1
+-10.000|1
+0.000|98
+3.142|1
+-10.000|1
+3.100|1
+3.100|1
+-10.000|1
+-10.000|1
+-10.000|1
+3.000|0
+3.140|1
+3.140|1
+-10.000|1
+-10.000|1
+-10.000|1
+3.142|1
+BEGIN
+-10.000|10
+0.000|980
+3.142|10
+-10.000|10
+1.000|0
+-10.000|10
+3.100|10
+-1.000|0
+3.100|10
+-10.000|10
+3.000|10
+-1.000|0
+3.100|10
+-9.000|0
+COMMIT
projects / postgis / commitdiff