return rtn;
}
+struct rt_reclassexpr_t {
+ struct rt_reclassrange {
+ double min;
+ double max;
+ int inc_min; /* include min */
+ int inc_max; /* include max */
+ int exc_min; /* exceed min */
+ int exc_max; /* exceed max */
+ } src, dst;
+};
+
+/**
+ * Returns new band with values reclassified
+ *
+ * @param srcband : the band who's values will be reclassified
+ * @param pixtype : pixel type of the new band
+ * @param hasnodata : indicates if the band has a nodata value
+ * @param nodataval : nodata value for the new band
+ * @param exprset : array of rt_reclassexpr structs
+ * @param exprcount : number of elements in expr
+ *
+ * @return a new rt_band or 0 on error
+ */
+rt_band
+rt_band_reclass(rt_band srcband, rt_pixtype pixtype,
+ uint32_t hasnodata, double nodataval, rt_reclassexpr *exprset,
+ int exprcount) {
+ rt_band band = NULL;
+ int width = 0;
+ int height = 0;
+ int numval = 0;
+ int memsize = 0;
+ void *mem = NULL;
+ uint32_t src_hasnodata = 0;
+ double src_nodataval = 0.0;
+
+ int rtn;
+ int x;
+ int y;
+ int i;
+ double or = 0;
+ double ov = 0;
+ double nr = 0;
+ double nv = 0;
+ int do_nv = 0;
+ rt_reclassexpr expr = NULL;
+
+ assert(NULL != srcband);
+ assert(NULL != exprset);
+
+ /* source nodata */
+ src_hasnodata = rt_band_get_hasnodata_flag(srcband);
+ src_nodataval = rt_band_get_nodata(srcband);
+
+ /* size of memory block to allocate */
+ width = rt_band_get_width(srcband);
+ height = rt_band_get_height(srcband);
+ numval = width * height;
+ memsize = rt_pixtype_size(pixtype) * numval;
+ mem = (int *) rtalloc(memsize);
+ if (!mem) {
+ rterror("rt_band_reclass: Could not allocate memory for band");
+ return 0;
+ }
+
+ /* initialize to zero */
+ if (!hasnodata) {
+ memset(mem, 0, memsize);
+ }
+ /* initialize to nodataval */
+ else {
+ int32_t checkvalint = 0;
+ uint32_t checkvaluint = 0;
+ double checkvaldouble = 0;
+ float checkvalfloat = 0;
+
+ switch (pixtype) {
+ case PT_1BB:
+ {
+ uint8_t *ptr = mem;
+ uint8_t clamped_initval = rt_util_clamp_to_1BB(nodataval);
+ for (i = 0; i < numval; i++)
+ ptr[i] = clamped_initval;
+ checkvalint = ptr[0];
+ break;
+ }
+ case PT_2BUI:
+ {
+ uint8_t *ptr = mem;
+ uint8_t clamped_initval = rt_util_clamp_to_2BUI(nodataval);
+ for (i = 0; i < numval; i++)
+ ptr[i] = clamped_initval;
+ checkvalint = ptr[0];
+ break;
+ }
+ case PT_4BUI:
+ {
+ uint8_t *ptr = mem;
+ uint8_t clamped_initval = rt_util_clamp_to_4BUI(nodataval);
+ for (i = 0; i < numval; i++)
+ ptr[i] = clamped_initval;
+ checkvalint = ptr[0];
+ break;
+ }
+ case PT_8BSI:
+ {
+ int8_t *ptr = mem;
+ int8_t clamped_initval = rt_util_clamp_to_8BSI(nodataval);
+ for (i = 0; i < numval; i++)
+ ptr[i] = clamped_initval;
+ checkvalint = ptr[0];
+ break;
+ }
+ case PT_8BUI:
+ {
+ uint8_t *ptr = mem;
+ uint8_t clamped_initval = rt_util_clamp_to_8BUI(nodataval);
+ for (i = 0; i < numval; i++)
+ ptr[i] = clamped_initval;
+ checkvalint = ptr[0];
+ break;
+ }
+ case PT_16BSI:
+ {
+ int16_t *ptr = mem;
+ int16_t clamped_initval = rt_util_clamp_to_16BSI(nodataval);
+ for (i = 0; i < numval; i++)
+ ptr[i] = clamped_initval;
+ checkvalint = ptr[0];
+ break;
+ }
+ case PT_16BUI:
+ {
+ uint16_t *ptr = mem;
+ uint16_t clamped_initval = rt_util_clamp_to_16BUI(nodataval);
+ for (i = 0; i < numval; i++)
+ ptr[i] = clamped_initval;
+ checkvalint = ptr[0];
+ break;
+ }
+ case PT_32BSI:
+ {
+ int32_t *ptr = mem;
+ int32_t clamped_initval = rt_util_clamp_to_32BSI(nodataval);
+ for (i = 0; i < numval; i++)
+ ptr[i] = clamped_initval;
+ checkvalint = ptr[0];
+ break;
+ }
+ case PT_32BUI:
+ {
+ uint32_t *ptr = mem;
+ uint32_t clamped_initval = rt_util_clamp_to_32BUI(nodataval);
+ for (i = 0; i < numval; i++)
+ ptr[i] = clamped_initval;
+ checkvaluint = ptr[0];
+ break;
+ }
+ case PT_32BF:
+ {
+ float *ptr = mem;
+ float clamped_initval = rt_util_clamp_to_32F(nodataval);
+ for (i = 0; i < numval; i++)
+ ptr[i] = clamped_initval;
+ checkvalfloat = ptr[0];
+ break;
+ }
+ case PT_64BF:
+ {
+ double *ptr = mem;
+ for (i = 0; i < numval; i++)
+ ptr[i] = nodataval;
+ checkvaldouble = ptr[0];
+ break;
+ }
+ default:
+ {
+ rterror("rt_band_reclass: Unknown pixeltype %d", pixtype);
+ rtdealloc(mem);
+ return 0;
+ }
+ }
+
+#ifdef POSTGIS_RASTER_WARN_ON_TRUNCATION
+ /* Overflow checking */
+ rt_util_display_dbl_trunc_warning(nodataval, checkvalint, checkvaluint, checkvalfloat,
+ checkvaldouble, pixtype);
+#endif /* POSTGIS_RASTER_WARN_ON_TRUNCATION */
+ }
+ RASTER_DEBUGF(3, "rt_band_reclass: width = %d height = %d", width, height);
+
+ band = rt_band_new_inline(width, height, pixtype, hasnodata, nodataval, mem);
+ if (!band) {
+ rterror("rt_band_reclass: Could not create new band");
+ rtdealloc(mem);
+ return 0;
+ }
+ RASTER_DEBUGF(3, "rt_band_reclass: new band @ %p", band);
+
+ for (x = 0; x < width; x++) {
+ for (y = 0; y < height; y++) {
+ rtn = rt_band_get_pixel(srcband, x, y, &ov);
+
+ /* error getting value, skip */
+ if (rtn == -1) continue;
+
+ do {
+ do_nv = 0;
+
+ /* no data*/
+ if (src_hasnodata && hasnodata && ov == src_nodataval) {
+ do_nv = 1;
+ break;
+ }
+
+ for (i = 0; i < exprcount; i++) {
+ expr = exprset[i];
+
+ /* ov matches min and max*/
+ if (expr->src.min == ov && ov == expr->src.max) {
+ do_nv = 1;
+ break;
+ }
+
+ /* process min */
+ if (
+ (expr->src.exc_min && expr->src.min >= ov) ||
+ (expr->src.inc_min && expr->src.min <= ov) ||
+ (expr->src.min < ov)
+ ) {
+ /* process max */
+ if (
+ (expr->src.exc_max && ov >= expr->src.max) ||
+ (expr->src.inc_max && ov <= expr->src.max) ||
+ (ov < expr->src.max)
+ ) {
+ do_nv = 1;
+ break;
+ }
+ }
+ }
+ }
+ while (0);
+
+ /* no expression matched, do not continue */
+ if (!do_nv) continue;
+
+ /* converting a value from one range to another range
+ OldRange = (OldMax - OldMin)
+ NewRange = (NewMax - NewMin)
+ NewValue = (((OldValue - OldMin) * NewRange) / OldRange) + NewMin
+ */
+
+ /* nodata */
+ if (src_hasnodata && hasnodata && ov == src_nodataval) {
+ nv = nodataval;
+ }
+ /*
+ "src" min and max is the same, prevent division by zero
+ set nv to "dst" min, which should be the same as "dst" max
+ */
+ else if (expr->src.min == expr->src.max) {
+ nv = expr->dst.min;
+ }
+ else {
+ or = expr->src.max - expr->src.min;
+ nr = expr->dst.max - expr->dst.min;
+ nv = (((ov - expr->src.min) * nr) / or) + expr->dst.min;
+
+ if (nv < expr->dst.min)
+ nv = expr->dst.min;
+ else if (nv > expr->dst.max)
+ nv = expr->dst.max;
+ }
+
+ /* round the value for integers */
+ switch (pixtype) {
+ 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:
+ nv = round(nv);
+ break;
+ default:
+ break;
+ }
+
+ RASTER_DEBUGF(5, "(%d, %d) ov: %f or: %f - %f nr: %f - %f nv: %f"
+ , x
+ , y
+ , ov
+ , (NULL != expr) ? expr->src.min : 0
+ , (NULL != expr) ? expr->src.max : 0
+ , (NULL != expr) ? expr->dst.min : 0
+ , (NULL != expr) ? expr->dst.max : 0
+ , nv
+ );
+ rtn = rt_band_set_pixel(band, x, y, nv);
+ if (rtn == -1) {
+ rterror("rt_band_reclass: Could not assign value to new band");
+ rt_band_destroy(band);
+ rtdealloc(mem);
+ return 0;
+ }
+
+ expr = NULL;
+ }
+ }
+
+ return band;
+}
+
/*- rt_raster --------------------------------------------------------*/
struct rt_raster_serialized_t {
return rast;
}
+/**
+ * Replace band at provided index with new band
+ *
+ * @param raster: raster of band to be replaced
+ * @param band : new band to add to raster
+ * @param index : index of band to replace (1-based)
+ *
+ * @return 0 on error or replaced band
+ */
+int
+rt_raster_replace_band(rt_raster raster, rt_band band, int index) {
+ rt_band oldband = NULL;
+ assert(NULL != raster);
+
+ if (band->width != raster->width || band->height != raster->height) {
+ rterror("rt_raster_replace_band: Band does not match raster's dimensions: %dx%d band to %dx%d raster",
+ band->width, band->height, raster->width, raster->height);
+ return 0;
+ }
+
+ if (index > raster->numBands || index < 0) {
+ rterror("rt_raster_replace_band: Band index is not valid");
+ return 0;
+ }
+
+ oldband = rt_raster_get_band(raster, index);
+ RASTER_DEBUGF(3, "rt_raster_replace_band: old band at %p", oldband);
+ RASTER_DEBUGF(3, "rt_raster_replace_band: new band at %p", band);
+
+ raster->bands[index] = band;
+ RASTER_DEBUGF(3, "rt_raster_replace_band: new band at %p", raster->bands[index]);
+
+ rt_band_destroy(oldband);
+ return 1;
+}
+
rt_quantile rt_band_get_quantiles(rt_bandstats stats,
double *quantiles, int quantiles_count, int *rtn_count);
+/**
+ * Returns new band with values reclassified
+ *
+ * @param srcband : the band who's values will be reclassified
+ * @param pixtype : pixel type of the new band
+ * @param hasnodata : indicates if the band has a nodata value
+ * @param nodataval : nodata value for the new band
+ * @param exprset : array of rt_reclassexpr structs
+ * @param exprcount : number of elements in expr
+ *
+ * @return a new rt_band or 0 on error
+ */
+typedef struct rt_reclassexpr_t* rt_reclassexpr;
+rt_band rt_band_reclass(rt_band srcband, rt_pixtype pixtype,
+ uint32_t hasnodata, double nodataval,
+ rt_reclassexpr *exprset, int exprcount);
/*- rt_raster --------------------------------------------------------*/
rt_raster rt_raster_from_band(rt_raster raster, uint32_t *bandNums,
int count);
+/**
+ * Replace band at provided index with new band
+ *
+ * @param raster: raster of band to be replaced
+ * @param band : new band to add to raster
+ * @param index : index of band to replace (1-based)
+ *
+ * @return 0 on error or 1 on success
+ */
+int rt_raster_replace_band(rt_raster raster, rt_band band,
+ int index);
+
/*- utilities -------------------------------------------------------*/
#include <float.h>
#include <string.h>
#include <stdio.h>
+#include <stdlib.h> /* for strtod in RASTER_reclass */
#include <errno.h>
#include <assert.h>
+#include <ctype.h> /* for isspace */
#include <postgres.h> /* for palloc */
#include <access/gist.h>
#include <utils/elog.h>
#include <utils/builtins.h>
#include <executor/spi.h>
+#include <executor/executor.h> /* for GetAttributeByName in RASTER_reclass */
#include <funcapi.h>
/*#include "lwgeom_pg.h"*/
static char * replace(const char *str, const char *oldstr, const char *newstr,
int *count);
static char *strtoupper(char *str);
+static char *chartrim(char* input, char *remove); /* for RASTER_reclass */
+static char **strsplit(const char *str, const char *delimiter, int *n); /* for RASTER_reclass */
+static char *removespaces(char *str); /* for RASTER_reclass */
/***************************************************************
* Some rules for returning NOTICE or ERROR...
/* get quantiles */
Datum RASTER_quantile(PG_FUNCTION_ARGS);
+/* reclassify specified bands of a raster */
+Datum RASTER_reclass(PG_FUNCTION_ARGS);
+
/* Replace function taken from
* http://ubuntuforums.org/showthread.php?s=aa6f015109fd7e4c7e30d2fd8b717497&t=141670&page=3
}
+static char*
+chartrim(char *input, char *remove) {
+ char *start;
+ char *ptr;
+
+ if (!input) {
+ return NULL;
+ }
+ else if (!*input) {
+ return input;
+ }
+
+ start = (char *) palloc(sizeof(char) * strlen(input) + 1);
+ if (NULL == start) {
+ fprintf(stderr, "Not enough memory\n");
+ return NULL;
+ }
+ strcpy(start, input);
+
+ /* trim left */
+ while (strchr(remove, *start) != NULL) {
+ start++;
+ }
+
+ /* trim right */
+ ptr = start + strlen(start);
+ while (strchr(remove, *--ptr) != NULL);
+ *(++ptr) = '\0';
+
+ return start;
+}
+
+/* split a string based on a delimiter */
+static char**
+strsplit(const char *str, const char *delimiter, int *n) {
+ char *tmp = NULL;
+ char **rtn = NULL;
+ char *token = NULL;
+
+ if (!str)
+ return NULL;
+
+ /* copy str to tmp as strtok will mangle the string */
+ tmp = palloc(sizeof(char) * (strlen(str) + 1));
+ if (NULL == tmp) {
+ fprintf(stderr, "Not enough memory\n");
+ return NULL;
+ }
+ strcpy(tmp, str);
+
+ if (!strlen(tmp) || !delimiter || !strlen(delimiter)) {
+ *n = 1;
+ rtn = (char **) palloc(*n * sizeof(char *));
+ if (NULL == rtn) {
+ fprintf(stderr, "Not enough memory\n");
+ return NULL;
+ }
+ rtn[0] = (char *) palloc(sizeof(char) * (strlen(tmp) + 1));
+ if (NULL == rtn[0]) {
+ fprintf(stderr, "Not enough memory\n");
+ return NULL;
+ }
+ strcpy(rtn[0], tmp);
+ pfree(tmp);
+ return rtn;
+ }
+
+ *n = 0;
+ token = strtok(tmp, delimiter);
+ while (token != NULL) {
+ if (*n < 1) {
+ rtn = (char **) palloc(sizeof(char *));
+ }
+ else {
+ rtn = (char **) repalloc(rtn, (*n + 1) * sizeof(char *));
+ }
+ if (NULL == rtn) {
+ fprintf(stderr, "Not enough memory\n");
+ return NULL;
+ }
+
+ rtn[*n] = NULL;
+ rtn[*n] = (char *) palloc(sizeof(char) * (strlen(token) + 1));
+ if (NULL == rtn[*n]) {
+ fprintf(stderr, "Not enough memory\n");
+ return NULL;
+ }
+
+ strcpy(rtn[*n], token);
+ *n = *n + 1;
+
+ token = strtok(NULL, delimiter);
+ }
+
+ pfree(tmp);
+ return rtn;
+}
+
+static char *
+removespaces(char *str) {
+ char *rtn;
+
+ rtn = replace(str, " ", "", NULL);
+ rtn = replace(rtn, "\n", "", NULL);
+ rtn = replace(rtn, "\t", "", NULL);
+ rtn = replace(rtn, "\f", "", NULL);
+ rtn = replace(rtn, "\r", "", NULL);
+
+ return rtn;
+}
+
PG_FUNCTION_INFO_V1(RASTER_lib_version);
Datum RASTER_lib_version(PG_FUNCTION_ARGS)
{
}
}
+struct rt_reclassexpr_t {
+ struct rt_reclassrange {
+ double min;
+ double max;
+ int inc_min;
+ int inc_max;
+ int exc_min;
+ int exc_max;
+ } src, dst;
+};
+
+/**
+ * Reclassify the specified bands of the raster
+ */
+PG_FUNCTION_INFO_V1(RASTER_reclass);
+Datum RASTER_reclass(PG_FUNCTION_ARGS) {
+ rt_pgraster *pgrast = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ rt_raster raster = NULL;
+ rt_band band = NULL;
+ rt_band newband = NULL;
+ uint32_t numBands = 0;
+
+ ArrayType *array;
+ Oid etype;
+ Datum *e;
+ bool *nulls;
+ int16 typlen;
+ bool typbyval;
+ char typalign;
+ int ndims = 1;
+ int *dims;
+ int *lbs;
+ int n = 0;
+
+ int i = 0;
+ int j = 0;
+ int k = 0;
+
+ int a = 0;
+ int b = 0;
+ int c = 0;
+
+ rt_reclassexpr *exprset = NULL;
+ HeapTupleHeader tup;
+ bool isnull;
+ Datum tupv;
+ uint32_t nband = 0;
+ char *expr = NULL;
+ text *exprtext = NULL;
+ double val = 0;
+ char *junk = NULL;
+ int inc_val = 0;
+ int exc_val = 0;
+ char *pixeltype = NULL;
+ text *pixeltypetext = NULL;
+ rt_pixtype pixtype = PT_END;
+ double nodataval = 0;
+ bool hasnodata = FALSE;
+
+ char **comma_set = NULL;
+ int comma_n = 0;
+ char **colon_set = NULL;
+ int colon_n = 0;
+ char **dash_set = NULL;
+ int dash_n = 0;
+
+ POSTGIS_RT_DEBUG(3, "RASTER_reclass: Starting");
+
+ /* raster */
+ raster = rt_raster_deserialize(pgrast);
+ if (!raster) {
+ elog(ERROR, "RASTER_reclass: Could not deserialize raster");
+ PG_RETURN_NULL();
+ }
+ numBands = rt_raster_get_num_bands(raster);
+ POSTGIS_RT_DEBUGF(3, "RASTER_reclass: %d possible bands to be reclassified", n);
+
+ /* process set of reclassarg */
+ POSTGIS_RT_DEBUG(3, "RASTER_reclass: Processing Arg 1 (reclassargset)");
+ array = PG_GETARG_ARRAYTYPE_P(1);
+ etype = ARR_ELEMTYPE(array);
+ get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
+
+ ndims = ARR_NDIM(array);
+ dims = ARR_DIMS(array);
+ lbs = ARR_LBOUND(array);
+
+ deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
+ &nulls, &n);
+
+ if (!n) {
+ elog(NOTICE, "Invalid argument for reclassargset. Returning original raster");
+ rt_raster_destroy(raster);
+
+ SET_VARSIZE(pgrast, pgrast->size);
+ PG_RETURN_POINTER(pgrast);
+ }
+
+ /*
+ process each element of reclassarg
+ each element is the index of the band to process, the set
+ of reclass ranges and the output band's pixeltype
+ */
+ for (i = 0; i < n; i++) {
+ if (nulls[i]) continue;
+
+ /* each element is a tuple */
+ tup = (HeapTupleHeader) DatumGetPointer(e[i]);
+ if (NULL == tup) {
+ elog(NOTICE, "Invalid argument for reclassargset. Returning original raster");
+ rt_raster_destroy(raster);
+
+ SET_VARSIZE(pgrast, pgrast->size);
+ PG_RETURN_POINTER(pgrast);
+ }
+
+ /* band index (1-based) */
+ tupv = GetAttributeByName(tup, "nband", &isnull);
+ if (isnull) {
+ elog(NOTICE, "Invalid argument for reclassargset. Missing value of nband for reclassarg of index %d . Returning original raster", i);
+ rt_raster_destroy(raster);
+
+ SET_VARSIZE(pgrast, pgrast->size);
+ PG_RETURN_POINTER(pgrast);
+ }
+ nband = DatumGetInt32(tupv);
+ POSTGIS_RT_DEBUGF(3, "RASTER_reclass: expression for band %d", nband);
+
+ /* valid band index? */
+ if (nband < 1 || nband > numBands) {
+ elog(NOTICE, "Invalid argument for reclassargset. Invalid band index (must use 1-based) for reclassarg of index %d . Returning original raster", i);
+ rt_raster_destroy(raster);
+
+ SET_VARSIZE(pgrast, pgrast->size);
+ PG_RETURN_POINTER(pgrast);
+ }
+
+ /* reclass expr */
+ tupv = GetAttributeByName(tup, "reclassexpr", &isnull);
+ if (isnull) {
+ elog(NOTICE, "Invalid argument for reclassargset. Missing value of reclassexpr for reclassarg of index %d . Returning original raster", i);
+ rt_raster_destroy(raster);
+
+ SET_VARSIZE(pgrast, pgrast->size);
+ PG_RETURN_POINTER(pgrast);
+ }
+ exprtext = (text *) DatumGetPointer(tupv);
+ if (NULL == exprtext) {
+ elog(NOTICE, "Invalid argument for reclassargset. Missing value of reclassexpr for reclassarg of index %d . Returning original raster", i);
+ rt_raster_destroy(raster);
+
+ SET_VARSIZE(pgrast, pgrast->size);
+ PG_RETURN_POINTER(pgrast);
+ }
+ expr = text_to_cstring(exprtext);
+ POSTGIS_RT_DEBUGF(5, "RASTER_reclass: expr (raw) %s", expr);
+ expr = removespaces(expr);
+ POSTGIS_RT_DEBUGF(5, "RASTER_reclass: expr (clean) %s", expr);
+
+ /* split string to its components */
+ /* comma (,) separating rangesets */
+ comma_set = strsplit(expr, ",", &comma_n);
+ if (comma_n < 1) {
+ elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
+ rt_raster_destroy(raster);
+
+ SET_VARSIZE(pgrast, pgrast->size);
+ PG_RETURN_POINTER(pgrast);
+ }
+
+ /* set of reclass expressions */
+ POSTGIS_RT_DEBUGF(5, "RASTER_reclass: %d possible expressions", comma_n);
+ exprset = palloc(comma_n * sizeof(rt_reclassexpr));
+
+ for (a = 0, j = 0; a < comma_n; a++) {
+ POSTGIS_RT_DEBUGF(5, "RASTER_reclass: map %s", comma_set[a]);
+
+ /* colon (:) separating range "src" and "dst" */
+ colon_set = strsplit(comma_set[a], ":", &colon_n);
+ if (colon_n != 2) {
+ elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
+ for (k = 0; k < j; k++) pfree(exprset[k]);
+ pfree(exprset);
+ rt_raster_destroy(raster);
+
+ SET_VARSIZE(pgrast, pgrast->size);
+ PG_RETURN_POINTER(pgrast);
+ }
+
+ /* allocate mem for reclass expression */
+ exprset[j] = palloc(sizeof(struct rt_reclassexpr_t));
+
+ for (b = 0; b < colon_n; b++) {
+ POSTGIS_RT_DEBUGF(5, "RASTER_reclass: range %s", colon_set[b]);
+
+ /* dash (-) separating "min" and "max" */
+ dash_set = strsplit(colon_set[b], "-", &dash_n);
+ if (dash_n < 1 || dash_n > 3) {
+ elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
+ for (k = 0; k < j; k++) pfree(exprset[k]);
+ pfree(exprset);
+ rt_raster_destroy(raster);
+
+ SET_VARSIZE(pgrast, pgrast->size);
+ PG_RETURN_POINTER(pgrast);
+ }
+
+ for (c = 0; c < dash_n; c++) {
+ /* need to handle: (-9999-100 -> "(", "9999", "100" */
+ if (
+ c < 1 &&
+ strlen(dash_set[c]) == 1 && (
+ strchr(dash_set[c], '(') != NULL ||
+ strchr(dash_set[c], '[') != NULL ||
+ strchr(dash_set[c], ')') != NULL ||
+ strchr(dash_set[c], ']') != NULL
+ )
+ ) {
+ junk = palloc(sizeof(char) * (strlen(dash_set[c + 1]) + 2));
+ if (NULL == junk) {
+ elog(ERROR, "RASTER_reclass: Insufficient memory. Returning NULL");
+ for (k = 0; k <= j; k++) pfree(exprset[k]);
+ pfree(exprset);
+ rt_raster_destroy(raster);
+
+ PG_RETURN_NULL();
+ }
+
+ sprintf(junk, "%s%s", dash_set[c], dash_set[c + 1]);
+ c++;
+ dash_set[c] = repalloc(dash_set[c], sizeof(char) * (strlen(junk) + 1));
+ strcpy(dash_set[c], junk);
+ pfree(junk);
+
+ /* rebuild dash_set */
+ for (k = 1; k < dash_n; k++) {
+ dash_set[k - 1] = repalloc(dash_set[k - 1], (strlen(dash_set[k]) + 1) * sizeof(char));
+ strcpy(dash_set[k - 1], dash_set[k]);
+ }
+ dash_n--;
+ c--;
+ pfree(dash_set[dash_n]);
+ dash_set = repalloc(dash_set, sizeof(char *) * dash_n);
+ }
+
+ /* there shouldn't be more than two in dash_n */
+ if (c < 1 && dash_n > 2) {
+ elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
+ for (k = 0; k < j; k++) pfree(exprset[k]);
+ pfree(exprset);
+ rt_raster_destroy(raster);
+
+ SET_VARSIZE(pgrast, pgrast->size);
+ PG_RETURN_POINTER(pgrast);
+ }
+
+ /* check interval flags */
+ exc_val = 0;
+ inc_val = 1;
+ /* range */
+ if (dash_n != 1) {
+ /* min */
+ if (c < 1) {
+ if (
+ strchr(dash_set[c], ')') != NULL ||
+ strchr(dash_set[c], ']') != NULL
+ ) {
+ exc_val = 1;
+ inc_val = 1;
+ }
+ else if (strchr(dash_set[c], '(') != NULL){
+ inc_val = 0;
+ }
+ else {
+ inc_val = 1;
+ }
+ }
+ /* max */
+ else {
+ if (
+ strrchr(dash_set[c], '(') != NULL ||
+ strrchr(dash_set[c], '[') != NULL
+ ) {
+ exc_val = 1;
+ inc_val = 0;
+ }
+ else if (strrchr(dash_set[c], ']') != NULL) {
+ inc_val = 1;
+ }
+ else {
+ inc_val = 0;
+ }
+ }
+ }
+ POSTGIS_RT_DEBUGF(5, "RASTER_reclass: exc_val %d inc_val %d", exc_val, inc_val);
+
+ /* remove interval flags */
+ dash_set[c] = chartrim(dash_set[c], "()[]");
+ POSTGIS_RT_DEBUGF(5, "RASTER_reclass: min/max (char) %s", dash_set[c]);
+
+ /* value from string to double */
+ errno = 0;
+ val = strtod(dash_set[c], &junk);
+ if (errno != 0 || dash_set[c] == junk) {
+ elog(NOTICE, "Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
+ for (k = 0; k < j; k++) pfree(exprset[k]);
+ pfree(exprset);
+ rt_raster_destroy(raster);
+
+ SET_VARSIZE(pgrast, pgrast->size);
+ PG_RETURN_POINTER(pgrast);
+ }
+ POSTGIS_RT_DEBUGF(5, "RASTER_reclass: min/max (double) %f", val);
+
+ /* strsplit removes dash (a.k.a. negative sign), compare now to restore */
+ junk = strstr(colon_set[b], dash_set[c]);
+ /* not beginning of string */
+ if (junk != colon_set[b]) {
+ /* prior is a dash */
+ if (*(junk - 1) == '-') {
+ /* prior is beginning of string or prior - 1 char is dash, negative number */
+ if (
+ ((junk - 1) == colon_set[b]) ||
+ (*(junk - 2) == '-') ||
+ (*(junk - 2) == '[') ||
+ (*(junk - 2) == '(')
+ ) {
+ val *= -1.;
+ }
+ }
+ }
+ POSTGIS_RT_DEBUGF(5, "RASTER_reclass: min/max (double) %f", val);
+
+ /* src */
+ if (b < 1) {
+ /* singular value */
+ if (dash_n == 1) {
+ exprset[j]->src.exc_min = exprset[j]->src.exc_max = exc_val;
+ exprset[j]->src.inc_min = exprset[j]->src.inc_max = inc_val;
+ exprset[j]->src.min = exprset[j]->src.max = val;
+ }
+ /* min */
+ else if (c < 1) {
+ exprset[j]->src.exc_min = exc_val;
+ exprset[j]->src.inc_min = inc_val;
+ exprset[j]->src.min = val;
+ }
+ /* max */
+ else {
+ exprset[j]->src.exc_max = exc_val;
+ exprset[j]->src.inc_max = inc_val;
+ exprset[j]->src.max = val;
+ }
+ }
+ /* dst */
+ else {
+ /* singular value */
+ if (dash_n == 1)
+ exprset[j]->dst.min = exprset[j]->dst.max = val;
+ /* min */
+ else if (c < 1)
+ exprset[j]->dst.min = val;
+ /* max */
+ else
+ exprset[j]->dst.max = val;
+ }
+ }
+ }
+
+ POSTGIS_RT_DEBUGF(3, "RASTER_reclass: or: %f - %f nr: %f - %f"
+ , exprset[j]->src.min
+ , exprset[j]->src.max
+ , exprset[j]->dst.min
+ , exprset[j]->dst.max
+ );
+ j++;
+ }
+
+ /* pixel type */
+ tupv = GetAttributeByName(tup, "pixeltype", &isnull);
+ if (isnull) {
+ elog(NOTICE, "Invalid argument for reclassargset. Missing value of pixeltype for reclassarg of index %d . Returning original raster", i);
+ rt_raster_destroy(raster);
+
+ SET_VARSIZE(pgrast, pgrast->size);
+ PG_RETURN_POINTER(pgrast);
+ }
+ pixeltypetext = (text *) DatumGetPointer(tupv);
+ if (NULL == pixeltypetext) {
+ elog(NOTICE, "Invalid argument for reclassargset. Missing value of pixeltype for reclassarg of index %d . Returning original raster", i);
+ rt_raster_destroy(raster);
+
+ SET_VARSIZE(pgrast, pgrast->size);
+ PG_RETURN_POINTER(pgrast);
+ }
+ pixeltype = text_to_cstring(pixeltypetext);
+ POSTGIS_RT_DEBUGF(3, "RASTER_reclass: pixeltype %s", pixeltype);
+ pixtype = rt_pixtype_index_from_name(pixeltype);
+
+ /* nodata */
+ tupv = GetAttributeByName(tup, "nodataval", &isnull);
+ if (isnull) {
+ nodataval = 0;
+ hasnodata = FALSE;
+ }
+ else {
+ nodataval = DatumGetFloat8(tupv);
+ hasnodata = TRUE;
+ }
+ POSTGIS_RT_DEBUGF(3, "RASTER_reclass: nodataval %f", nodataval);
+ POSTGIS_RT_DEBUGF(3, "RASTER_reclass: hasnodata %d", hasnodata);
+
+ /* do reclass */
+ band = rt_raster_get_band(raster, nband - 1);
+ if (!band) {
+ elog(NOTICE, "Could not find raster band of index %d. Returning original raster", nband);
+ for (k = 0; k < j; k++) pfree(exprset[k]);
+ pfree(exprset);
+ rt_raster_destroy(raster);
+
+ SET_VARSIZE(pgrast, pgrast->size);
+ PG_RETURN_POINTER(pgrast);
+ }
+ newband = rt_band_reclass(band, pixtype, hasnodata, nodataval, exprset, j);
+ if (!newband) {
+ elog(ERROR, "RASTER_reclass: Could not reclassify raster band of index %d. Returning NULL", nband);
+ rt_raster_destroy(raster);
+ for (k = 0; k < j; k++) pfree(exprset[k]);
+ pfree(exprset);
+
+ PG_RETURN_NULL();
+ }
+
+ /* replace old band with new band */
+ if (rt_raster_replace_band(raster, newband, nband - 1) == 0) {
+ elog(ERROR, "RASTER_reclass: Could not replace raster band of index %d with reclassified band. Returning NULL", nband);
+ rt_band_destroy(newband);
+ rt_raster_destroy(raster);
+ for (k = 0; k < j; k++) pfree(exprset[k]);
+ pfree(exprset);
+
+ PG_RETURN_NULL();
+ }
+
+ /* free exprset */
+ for (k = 0; k < j; k++) pfree(exprset[k]);
+ pfree(exprset);
+ }
+
+ pgrast = rt_raster_serialize(raster);
+ rt_raster_destroy(raster);
+
+ if (!pgrast)
+ PG_RETURN_NULL();
+
+ POSTGIS_RT_DEBUG(3, "RASTER_reclass: Finished");
+ SET_VARSIZE(pgrast, pgrast->size);
+ PG_RETURN_POINTER(pgrast);
+}
+
/* ---------------------------------------------------------------- */
/* Memory allocation / error reporting hooks */
/* ---------------------------------------------------------------- */
AS $$ SELECT quantile, value FROM _st_quantile($1, 1, $2, 1, ARRAY[$3]::double precision[]) $$
LANGUAGE 'sql' IMMUTABLE STRICT;
+-----------------------------------------------------------------------
+-- ST_Reclass
+-----------------------------------------------------------------------
+CREATE TYPE reclassarg AS (
+ nband int,
+ reclassexpr text,
+ pixeltype text,
+ nodataval double precision
+);
+
+CREATE OR REPLACE FUNCTION _st_reclass(rast raster, VARIADIC reclassargset reclassarg[])
+ RETURNS raster
+ AS 'MODULE_PATHNAME', 'RASTER_reclass'
+ LANGUAGE 'C' IMMUTABLE STRICT;
+
+CREATE OR REPLACE FUNCTION st_reclass(rast raster, VARIADIC reclassargset reclassarg[])
+ RETURNS raster
+ AS $$
+ DECLARE
+ i int;
+ expr text;
+ BEGIN
+ -- for each reclassarg, validate elements as all except nodataval cannot be NULL
+ FOR i IN SELECT * FROM generate_subscripts($2, 1) LOOP
+ IF $2[i].nband IS NULL OR $2[i].reclassexpr IS NULL OR $2[i].pixeltype IS NULL THEN
+ RAISE WARNING 'Values are required for the nband, reclassexpr and pixeltype attributes.';
+ RETURN rast;
+ END IF;
+ END LOOP;
+
+ RETURN _st_reclass($1, VARIADIC $2);
+ END;
+ $$ LANGUAGE 'plpgsql' IMMUTABLE STRICT;
+
+CREATE OR REPLACE FUNCTION st_reclass(rast raster, nband int, reclassexpr text, pixeltype text, nodataval double precision)
+ RETURNS raster
+ AS $$ SELECT st_reclass($1, ROW($2, $3, $4, $5)) $$
+ LANGUAGE 'SQL' IMMUTABLE STRICT;
+
+CREATE OR REPLACE FUNCTION st_reclass(rast raster, nband int, reclassexpr text, pixeltype text)
+ RETURNS raster
+ AS $$ SELECT st_reclass($1, ROW($2, $3, $4, NULL)) $$
+ LANGUAGE 'SQL' IMMUTABLE STRICT;
+
+CREATE OR REPLACE FUNCTION st_reclass(rast raster, reclassexpr text, pixeltype text)
+ RETURNS raster
+ AS $$ SELECT st_reclass($1, ROW(1, $2, $3, NULL)) $$
+ LANGUAGE 'SQL' IMMUTABLE STRICT;
+
-----------------------------------------------------------------------
-- MapAlgebra
-----------------------------------------------------------------------
rt_raster_destroy(raster);
}
+static void testRasterReplaceBand() {
+ rt_raster raster;
+ rt_band band;
+ void* mem;
+ size_t datasize;
+ uint16_t width;
+ uint16_t height;
+ double nodata;
+ int rtn;
+
+ raster = rt_raster_new(10, 10);
+ assert(raster); /* or we're out of virtual memory */
+ band = addBand(raster, PT_8BUI, 0, 0);
+ CHECK(band);
+ band = addBand(raster, PT_8BUI, 1, 255);
+ CHECK(band);
+
+ width = rt_raster_get_width(raster);
+ height = rt_raster_get_height(raster);
+
+ datasize = rt_pixtype_size(PT_8BUI) * width * height;
+ mem = malloc(datasize);
+ band = rt_band_new_inline(width, height, PT_8BUI, 1, 1, mem);
+ assert(band);
+
+ rtn = rt_raster_replace_band(raster, band, 0);
+ CHECK(rtn);
+ nodata = rt_band_get_nodata(rt_raster_get_band(raster, 0));
+ CHECK((nodata == 1));
+
+ rt_band_destroy(band);
+ rt_raster_destroy(raster);
+ free(mem);
+}
+
+struct rt_reclassexpr_t {
+ struct rt_reclassrange {
+ double min;
+ double max;
+ int inc_min; /* include min */
+ int inc_max; /* include max */
+ int exc_min; /* exceed min */
+ int exc_max; /* exceed max */
+ } src, dst;
+};
+static void testBandReclass() {
+ rt_reclassexpr *exprset;
+
+ rt_raster raster;
+ rt_band band;
+ uint16_t x;
+ uint16_t y;
+ double nodata;
+ int cnt = 2;
+ int i = 0;
+ int rtn;
+ rt_band newband;
+ double val;
+
+ raster = rt_raster_new(100, 10);
+ assert(raster); /* or we're out of virtual memory */
+ band = addBand(raster, PT_16BUI, 0, 0);
+ CHECK(band);
+ rt_band_set_nodata(band, 0);
+
+ for (x = 0; x < 100; x++) {
+ for (y = 0; y < 10; y++) {
+ rtn = rt_band_set_pixel(band, x, y, x * y + (x + y));
+ CHECK((rtn != -1));
+ }
+ }
+
+ nodata = rt_band_get_nodata(band);
+ CHECK_EQUALS(nodata, 0);
+
+ exprset = malloc(cnt * sizeof(rt_reclassexpr));
+ assert(exprset);
+
+ for (i = 0; i < cnt; i++) {
+ exprset[i] = malloc(sizeof(struct rt_reclassexpr_t));
+ assert(exprset[i]);
+
+ if (i == 0) {
+ /* nodata */
+ exprset[i]->src.min = 0;
+ exprset[i]->src.inc_min = 0;
+ exprset[i]->src.exc_min = 0;
+
+ exprset[i]->src.max = 0;
+ exprset[i]->src.inc_max = 0;
+ exprset[i]->src.exc_max = 0;
+
+ exprset[i]->dst.min = 0;
+ exprset[i]->dst.max = 0;
+ }
+ else {
+ /* range */
+ exprset[i]->src.min = 0;
+ exprset[i]->src.inc_min = 0;
+ exprset[i]->src.exc_min = 0;
+
+ exprset[i]->src.max = 1000;
+ exprset[i]->src.inc_max = 1;
+ exprset[i]->src.exc_max = 0;
+
+ exprset[i]->dst.min = 1;
+ exprset[i]->dst.max = 255;
+ }
+ }
+
+ newband = rt_band_reclass(band, PT_8BUI, 0, 0, exprset, cnt);
+ CHECK(newband);
+
+ rtn = rt_band_get_pixel(newband, 0, 0, &val);
+ CHECK((rtn != -1));
+ CHECK_EQUALS(val, 0);
+
+ rtn = rt_band_get_pixel(newband, 49, 5, &val);
+ CHECK((rtn != -1));
+ CHECK_EQUALS(val, 77);
+
+ rtn = rt_band_get_pixel(newband, 99, 9, &val);
+ CHECK((rtn != -1));
+ CHECK_EQUALS(val, 255);
+
+
+ for (i = cnt - 1; i >= 0; i--)
+ free(exprset[i]);
+ free(exprset);
+ rt_band_destroy(band);
+ rt_band_destroy(newband);
+ rt_raster_destroy(raster);
+}
+
int
main()
{
testBandStats();
printf("Successfully tested band stats\n");
+ printf("Testing rt_raster_replace_band\n");
+ testRasterReplaceBand();
+ printf("Successfully tested rt_raster_replace_band\n");
+
+ printf("Testing rt_band_reclass\n");
+ testBandReclass();
+ printf("Successfully tested rt_band_reclass\n");
+
deepRelease(raster);
rt_utility.sql \
create_rt_mapalgebra_test.sql \
rt_mapalgebra.sql \
+ rt_reclass.sql \
$(NULL)
TEST_GIST = \
--- /dev/null
+SELECT
+ ST_Value(t.rast, 1, 1, 1),
+ ST_Value(t.rast, 1, 10, 10),
+ ST_Value(t.rast, 1, 2, 10),
+ ST_BandNoDataValue(rast, 1)
+FROM (
+ SELECT ST_Reclass(
+ ST_SetValue(
+ ST_SetValue(
+ ST_AddBand(
+ ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0, -1),
+ 1, '32BUI', 0, 0
+ ),
+ 1, 1, 499
+ ),
+ 10, 10, 12
+ ),
+ ROW(1, '0-100:1-10, 101-500:11-150,501 - 10000: 151-254', '8BUI', 255)
+ ) AS rast
+) AS t;
+SELECT
+ ST_Value(t.rast, 1, 1, 1),
+ ST_Value(t.rast, 1, 10, 10),
+ ST_Value(t.rast, 1, 2, 10),
+ ST_BandNoDataValue(rast, 1),
+ ST_Value(t.rast, 2, 1, 1),
+ ST_Value(t.rast, 2, 10, 10),
+ ST_Value(t.rast, 2, 2, 10),
+ ST_BandNoDataValue(rast, 2)
+FROM (
+ SELECT ST_Reclass(
+ ST_SetValue(
+ ST_SetValue(
+ ST_AddBand(
+ ST_AddBand(
+ ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0, -1),
+ 1, '8BUI', 0, 1
+ ),
+ 2, '32BUI', 0, 0
+ ),
+ 2, 1, 1, 499
+ ),
+ 2, 10, 10, 12
+ ),
+ ROW(2, '0-100:1-10, 101-500:11-150,501 - 10000: 151-254', '8BUI', 255)
+ ) AS rast
+) AS t;
+SELECT
+ ST_Value(t.rast, 1, 1, 1),
+ ST_Value(t.rast, 1, 10, 10),
+ ST_Value(t.rast, 1, 2, 10),
+ ST_BandNoDataValue(rast, 1),
+ ST_Value(t.rast, 2, 1, 1),
+ ST_Value(t.rast, 2, 10, 10),
+ ST_Value(t.rast, 2, 2, 10),
+ ST_BandNoDataValue(rast, 2)
+FROM (
+ SELECT ST_Reclass(
+ ST_SetValue(
+ ST_SetValue(
+ ST_AddBand(
+ ST_AddBand(
+ ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0, -1),
+ 1, '8BUI', 0, 1
+ ),
+ 2, '32BUI', 0, 0
+ ),
+ 2, 1, 1, 499
+ ),
+ 2, 10, 10, 12
+ ),
+ ROW(1, '0:1', '8BUI', 255),
+ ROW(2, '0-100:1-10, 101-500:11-150,501 - 10000: 151-254', '8BUI', 255)
+ ) AS rast
+) AS t;
+SELECT
+ ST_Value(t.rast, 1, 1, 1),
+ ST_Value(t.rast, 1, 10, 10),
+ ST_BandNoDataValue(rast, 1)
+FROM (
+ SELECT ST_Reclass(
+ ST_SetValue(
+ ST_SetValue(
+ ST_AddBand(
+ ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0, -1),
+ 1, '8BUI', 0, 0
+ ),
+ 1, 1, 1, 255
+ ),
+ 1, 10, 10, 0
+ ),
+ ROW(1, '0-100]:200-255,(100-200]:100-200,(200-255]:0-100', '8BUI', NULL)
+ ) AS rast
+) AS t;
+SELECT
+ ST_Value(t.rast, 1, 1, 1),
+ ST_Value(t.rast, 1, 10, 10),
+ ST_BandNoDataValue(rast, 1)
+FROM (
+ SELECT ST_Reclass(
+ ST_SetValue(
+ ST_SetValue(
+ ST_AddBand(
+ ST_MakeEmptyRaster(100, 100, 10, 10, 2, 2, 0, 0, -1),
+ 1, '32BF', 1, 0
+ ),
+ 1, 1, 1, 3.14159
+ ),
+ 1, 10, 10, 2.71828
+ ),
+ ROW(1, '-10000--100]:1-50,(-100-1000]:50-150,(1000-10000]:150-254', '8BUI', 0)
+ ) AS rast
+) AS t;
--- /dev/null
+150|2||255
+0|0|0|1|150|2||255
+1|1|1|255|150|2||255
+100|200|
+59|59|0
projects / postgis / commitdiff