From: Bborie Park Date: Mon, 16 May 2011 19:43:10 +0000 (+0000) Subject: Added ST_Reclass function X-Git-Tag: 2.0.0alpha1~1657 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=6eb30cd26241cabb7c933e28dd0f62de41da7445;p=postgis Added ST_Reclass function - added rt_band_reclass and rt_raster_replace_band to rt_core/rt_api.c and rt_api.h - added test case to test/core/testapi.c - added RASTER_reclass to rt_pg/rt_pg.c - added SQL functions for ST_Reclass - added regression tests Associated ticket is #903 git-svn-id: http://svn.osgeo.org/postgis/trunk@7154 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c index fdbbe8f11..1c7776a28 100644 --- a/raster/rt_core/rt_api.c +++ b/raster/rt_core/rt_api.c @@ -2119,6 +2119,323 @@ rt_band_get_quantiles(rt_bandstats stats, 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 { @@ -4595,3 +4912,39 @@ rt_raster_from_band(rt_raster raster, uint32_t *bandNums, int count) { 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; +} + diff --git a/raster/rt_core/rt_api.h b/raster/rt_core/rt_api.h index c97ad173d..bd87d6b3b 100644 --- a/raster/rt_core/rt_api.h +++ b/raster/rt_core/rt_api.h @@ -489,6 +489,22 @@ typedef struct rt_quantile_t* rt_quantile; 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 --------------------------------------------------------*/ @@ -817,6 +833,18 @@ int32_t rt_raster_copy_band(rt_raster torast, 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 -------------------------------------------------------*/ diff --git a/raster/rt_pg/rt_pg.c b/raster/rt_pg/rt_pg.c index f3baf0a7b..22008948b 100644 --- a/raster/rt_pg/rt_pg.c +++ b/raster/rt_pg/rt_pg.c @@ -30,8 +30,10 @@ #include #include #include +#include /* for strtod in RASTER_reclass */ #include #include +#include /* for isspace */ #include /* for palloc */ #include @@ -40,6 +42,7 @@ #include #include #include +#include /* for GetAttributeByName in RASTER_reclass */ #include /*#include "lwgeom_pg.h"*/ @@ -66,6 +69,9 @@ PG_MODULE_MAGIC; 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... @@ -196,6 +202,9 @@ Datum RASTER_histogram(PG_FUNCTION_ARGS); /* 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 @@ -267,6 +276,117 @@ strtoupper(char * str) } +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) { @@ -3537,6 +3657,466 @@ Datum RASTER_quantile(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 */ /* ---------------------------------------------------------------- */ diff --git a/raster/rt_pg/rtpostgis.sql.in.c b/raster/rt_pg/rtpostgis.sql.in.c index 503e2ed0d..e71233085 100644 --- a/raster/rt_pg/rtpostgis.sql.in.c +++ b/raster/rt_pg/rtpostgis.sql.in.c @@ -973,6 +973,55 @@ CREATE OR REPLACE FUNCTION st_approxquantile(rast raster, hasnodata boolean, qua 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 ----------------------------------------------------------------------- diff --git a/raster/test/core/testapi.c b/raster/test/core/testapi.c index c878ec8b5..a1cc396a7 100644 --- a/raster/test/core/testapi.c +++ b/raster/test/core/testapi.c @@ -1127,6 +1127,140 @@ static void testBandStats() { 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() { @@ -1441,6 +1575,14 @@ 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); diff --git a/raster/test/regress/Makefile.in b/raster/test/regress/Makefile.in index c9ed69ea5..9f594651f 100644 --- a/raster/test/regress/Makefile.in +++ b/raster/test/regress/Makefile.in @@ -78,6 +78,7 @@ TEST_UTILITY = \ rt_utility.sql \ create_rt_mapalgebra_test.sql \ rt_mapalgebra.sql \ + rt_reclass.sql \ $(NULL) TEST_GIST = \ diff --git a/raster/test/regress/rt_reclass.sql b/raster/test/regress/rt_reclass.sql new file mode 100644 index 000000000..ac9400c5d --- /dev/null +++ b/raster/test/regress/rt_reclass.sql @@ -0,0 +1,113 @@ +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; diff --git a/raster/test/regress/rt_reclass_expected b/raster/test/regress/rt_reclass_expected new file mode 100644 index 000000000..84d63610d --- /dev/null +++ b/raster/test/regress/rt_reclass_expected @@ -0,0 +1,5 @@ +150|2||255 +0|0|0|1|150|2||255 +1|1|1|255|150|2||255 +100|200| +59|59|0