[SCM] PostGIS branch master updated. 3.5.0-427-gc8de602d6
git at osgeo.org
git at osgeo.org
Fri Jul 4 13:38:01 PDT 2025
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "PostGIS".
The branch, master has been updated
via c8de602d6025adc6c16644d348e90e039dc5629f (commit)
from 784f21a1fd90246aba1c250a2729ccf2de29b7f1 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit c8de602d6025adc6c16644d348e90e039dc5629f
Author: Paul Ramsey <paul.ramsey at snowflake.com>
Date: Fri Jul 4 13:35:44 2025 -0700
Add ST_ReclassExact(raster, float8[], float8[]) returns raster
function to do 1:1 value mapping of rasters. Ordinary reclass
supports mapping ranges to ranges, which is complex and slow.
For the case where the mapping is 1:1, this function is 100x
faster.
diff --git a/NEWS b/NEWS
index f331baa3b..7952cfacb 100644
--- a/NEWS
+++ b/NEWS
@@ -37,3 +37,5 @@ Daniel Nylander (Swedish Team)
- #5841, Change interrupt handling to remove use of pqsignal to support PG 18 (Paul Ramsey)
- Add ST_CoverageClean to edge match and gap remove polygonal
coverages (Paul Ramsey) from GEOS 3.14 (Martin Davis)
+ - Add ST_ReclassExact to quickly remap values in raster (Paul Ramsey)
+
diff --git a/doc/reference_raster.xml b/doc/reference_raster.xml
index e283dd77c..46efd0632 100644
--- a/doc/reference_raster.xml
+++ b/doc/reference_raster.xml
@@ -5,7 +5,7 @@
<para><varname>raster</varname> is a PostGIS type for storing and analyzing raster data. </para>
<para>For loading rasters from raster files please refer to <xref linkend="RT_Loading_Rasters"/></para>
- <para>Some examples in this reference use a raster table of dummy data,
+ <para>Some examples in this reference use a raster table of dummy data,
created with the following code:</para>
<programlisting language="sql">CREATE TABLE dummy_rast(rid integer, rast raster);
INSERT INTO dummy_rast(rid, rast)
@@ -40,7 +40,7 @@ VALUES (1,
'41000000007719564100000000000000000000000000000000FFFFFFFF050005000400FDFEFDFEFEFDFEFEFDF9FAFEF' ||
'EFCF9FBFDFEFEFDFCFAFEFEFE04004E627AADD16076B4F9FE6370A9F5FE59637AB0E54F58617087040046566487A1506CA2E3FA5A6CAFFBFE4D566DA4CB3E454C5665')::raster);</programlisting>
- <para>The functions below are ones which a user of PostGIS Raster is likely to need.
+ <para>The functions below are ones which a user of PostGIS Raster is likely to need.
There are other raster support functions which are not of interest to a general
user.</para>
@@ -12392,13 +12392,13 @@ SELECT ST_MapAlgebraFctNgb(rast, 1, '8BUI', 4,4,
<refsection>
<title>Description</title>
- <para>Creates a new raster formed by applying a reclassification operation defined by the <varname>reclassexpr</varname> on the input raster (<varname>rast</varname>).
+ <para>Creates a new raster formed by applying a reclassification operation defined by the <varname>reclassexpr</varname> on the input raster (<varname>rast</varname>).
Refer to <xref linkend="reclassarg"/> for the description of reclassification expressions.
- If no <varname>band</varname> is specified band 1 is assumed.
+ If no <varname>band</varname> is specified band 1 is assumed.
</para>
- <para>The new raster will have the same georeference, width, and height as the original raster.
- The bands of the new raster have pixel type of <varname>pixeltype</varname>.
+ <para>The new raster will have the same georeference, width, and height as the original raster.
+ The bands of the new raster have pixel type of <varname>pixeltype</varname>.
If <varname>reclassargset</varname> is specified then each <varname>reclassarg</varname> defines the type of the target band.
Bands not designated are returned unchanged.
</para>
@@ -12494,6 +12494,117 @@ UPDATE wind
</refsection>
</refentry>
+ <refentry xml:id="RT_ST_ReclassExact">
+ <refnamediv>
+ <refname>ST_ReclassExact</refname>
+ <refpurpose>Creates a new raster composed of bands reclassified from original, using a 1:1 mapping from values in the original band to new values in the destination band.</refpurpose>
+ </refnamediv>
+
+ <refsynopsisdiv>
+ <funcsynopsis>
+ <funcprototype>
+ <funcdef>raster <function>ST_ReclassExact</function></funcdef>
+ <paramdef><type>raster </type> <parameter>rast</parameter></paramdef>
+ <paramdef><type>double precision[] </type> <parameter>inputvalues</parameter></paramdef>
+ <paramdef><type>double precision[] </type> <parameter>outputvalues</parameter></paramdef>
+ <paramdef choice="opt"><type>integer </type> <parameter>bandnumber=1</parameter></paramdef>
+ <paramdef choice="opt"><type>text </type> <parameter>pixeltype=32BF</parameter></paramdef>
+ <paramdef choice="opt"><type>double precision </type> <parameter>nodatavalue=NULL</parameter></paramdef>
+ </funcprototype>
+ </funcsynopsis>
+ </refsynopsisdiv>
+
+ <refsection>
+ <title>Description</title>
+
+ <para>Creates a new raster formed by applying a reclassification operation
+ defined by the <varname>inputvalues</varname> and
+ <varname>outputvalues</varname> arrays. Pixel values found in the
+ input array are mapped to the corresponding value in the output array.
+ All other pixel values are mapped to the <varname>nodatavalue</varname>.
+ </para>
+ <para>The output pixel type defaults to float, but can be specified
+ using the <varname>pixeltype</varname> parameter.
+ If no <varname>bandnumber</varname> is specified band 1 is assumed.
+ </para>
+
+ <para>The new raster will have the same georeference, width,
+ and height as the original raster.
+ Bands not designated are returned unchanged.
+ </para>
+
+ <para role="availability" conformance="3.6.0">Availability: 3.6.0 </para>
+ </refsection>
+
+ <refsection>
+ <title>Example: Basic</title>
+ <para>Create a small raster and map its pixels to new values.</para>
+ <programlisting>
+CREATE TABLE reclassexact (
+ id integer,
+ rast raster
+);
+
+--
+-- Create a raster with just four pixels
+-- [1 2]
+-- [3 4]
+--
+INSERT INTO reclassexact (id, rast)
+SELECT 1, ST_SetValues(
+ ST_AddBand(
+ ST_MakeEmptyRaster(
+ 2, -- width in pixels
+ 2, -- height in pixels
+ 0, -- upper-left x-coordinate
+ 0, -- upper-left y-coordinate
+ 1, -- pixel size in x-direction
+ -1, -- pixel size in y-direction (negative for north-up)
+ 0, -- skew in x-direction
+ 0, -- skew in y-direction
+ 4326 -- SRID (e.g., WGS 84)
+ ),
+ '32BUI'::text, -- pixel type (e.g., '32BF' for float, '8BUI' for unsigned 8-bit int)
+ 0.0, -- initial value for the band (e.g., 0.0 or a no-data value)
+ -99 -- nodatavalue
+ ),
+ 1, -- band number (usually 1 for single-band rasters)
+ 1, -- x origin for setting values (usually 1)
+ 1, -- y origin for setting values (usually 1)
+ ARRAY[
+ ARRAY[1, 2],
+ ARRAY[3, 4]
+ ]::double precision[][] -- 2D array of values
+ );
+
+-- Reclass the values to new values
+-- and dump the values of the new raster for display
+WITH rc AS (
+ SELECT ST_ReclassExact(
+ rast, -- input raster
+ ARRAY[4,3,2,1], -- input map
+ ARRAY[14,13,12,11], -- output map
+ 1, -- band number to remap
+ '32BUI' -- output raster pixtype
+ ) AS rast
+ FROM reclassexact
+ WHERE id = 1
+ )
+SELECT 'rce-1', (ST_DumpValues(rc.rast)).*
+FROM rc;
+ </programlisting>
+ </refsection>
+ <refsection>
+ <title>See Also</title>
+ <para>
+ <xref linkend="RT_ST_Reclass"/>,
+ <xref linkend="RT_ST_AddBand"/>,
+ <xref linkend="RT_ST_Band"/>,
+ <xref linkend="RT_ST_MakeEmptyRaster"/>
+ </para>
+ </refsection>
+ </refentry>
+
<refentry xml:id="RT_ST_Union">
<refnamediv>
<refname>ST_Union</refname>
diff --git a/raster/rt_core/librtcore.h b/raster/rt_core/librtcore.h
index beced79c4..399f9d8e2 100644
--- a/raster/rt_core/librtcore.h
+++ b/raster/rt_core/librtcore.h
@@ -154,6 +154,7 @@ typedef struct rt_quantile_t* rt_quantile;
typedef struct rt_valuecount_t* rt_valuecount;
typedef struct rt_gdaldriver_t* rt_gdaldriver;
typedef struct rt_reclassexpr_t* rt_reclassexpr;
+typedef struct rt_reclassmap_t* rt_reclassmap;
typedef struct rt_iterator_t* rt_iterator;
typedef struct rt_iterator_arg_t* rt_iterator_arg;
@@ -435,6 +436,18 @@ rt_band rt_band_new_inline(
uint8_t* data
);
+/**
+ * Fill in the cells of a band with a starting value
+ * frequently used to init with nodata value
+ *
+ * @param band : band to initialize
+ * @param initval : value to initialize with
+ */
+void rt_band_init_value(
+ rt_band band,
+ double initval
+);
+
/**
* Create an out-db rt_band
*
@@ -997,6 +1010,20 @@ rt_band rt_band_reclass(
rt_reclassexpr *exprset, int exprcount
);
+/**
+ * Returns new band with values reclassified
+ * @param srcband : the band who's values will be reclassified
+ * @param map : the src and dst values to remapping
+ * @param hasnodata : indicates if the band has a nodata value
+ * @param nodataval : nodata value for the new band
+ *
+ * @return a new rt_band or NULL on error
+ */
+rt_band rt_band_reclass_exact(
+ rt_band srcband, rt_reclassmap map,
+ uint32_t hasnodata, double nodataval
+);
+
/*- rt_raster --------------------------------------------------------*/
/**
@@ -2582,6 +2609,20 @@ struct rt_reclassexpr_t {
} src, dst;
};
+/* src/dst mapping for rt_classmap */
+struct rt_classpair_t {
+ double src;
+ double dst;
+};
+
+/* exact reclassification expression */
+struct rt_reclassmap_t {
+ uint32_t count;
+ rt_pixtype srctype;
+ rt_pixtype dsttype;
+ struct rt_classpair_t *pairs;
+};
+
/* raster iterator */
struct rt_iterator_t {
rt_raster raster;
diff --git a/raster/rt_core/rt_band.c b/raster/rt_core/rt_band.c
index f633bc822..4578609a2 100644
--- a/raster/rt_core/rt_band.c
+++ b/raster/rt_core/rt_band.c
@@ -101,6 +101,157 @@ rt_band_new_inline(
return band;
}
+/**
+ * Fill in the cells of a band with a starting value
+ * frequently used to init with nodata value
+ *
+ * @param band : band to initialize
+ * @param initval : value to initialize with
+ */
+void
+rt_band_init_value(
+ rt_band band,
+ double initval
+) {
+ assert(band != NULL);
+ assert(band->data.mem != NULL);
+ rt_pixtype pixtype = band->pixtype;
+ uint32_t width = band->width;
+ uint32_t height = band->height;
+ uint32_t numval = width * height;
+ void *mem = band->data.mem;
+ size_t memsize = numval * rt_pixtype_size(pixtype);
+
+ /* initialize to nodataval */
+ int32_t checkvalint = 0;
+ uint32_t checkvaluint = 0;
+ double checkvaldouble = 0;
+ float checkvalfloat = 0;
+
+ /* initialize to zero */
+ if (FLT_EQ(initval, 0.0)) {
+ memset(mem, 0, memsize);
+ return;
+ }
+
+ switch (pixtype) {
+ case PT_1BB:
+ {
+ uint8_t *ptr = mem;
+ uint8_t clamped_initval = rt_util_clamp_to_1BB(initval);
+ for (uint32_t 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(initval);
+ for (uint32_t 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(initval);
+ for (uint32_t 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(initval);
+ for (uint32_t 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(initval);
+ for (uint32_t 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(initval);
+ for (uint32_t 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(initval);
+ for (uint32_t 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(initval);
+ for (uint32_t 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(initval);
+ for (uint32_t 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(initval);
+ for (uint32_t i = 0; i < numval; i++)
+ ptr[i] = clamped_initval;
+ checkvalfloat = ptr[0];
+ break;
+ }
+ case PT_64BF:
+ {
+ double *ptr = mem;
+ for (uint32_t i = 0; i < numval; i++)
+ ptr[i] = initval;
+ checkvaldouble = ptr[0];
+ break;
+ }
+ default:
+ {
+ rterror("%s: Unknown pixeltype %d", __func__, pixtype);
+ return;
+ }
+ }
+
+ /* Overflow checking */
+ rt_util_dbl_trunc_warning(
+ initval,
+ checkvalint, checkvaluint,
+ checkvalfloat, checkvaldouble,
+ pixtype
+ );
+
+ return;
+}
+
+
/**
* Create an out-db rt_band
*
@@ -2080,4 +2231,3 @@ rt_band_corrected_clamped_value(
return ES_NONE;
}
-
diff --git a/raster/rt_core/rt_mapalgebra.c b/raster/rt_core/rt_mapalgebra.c
index 4d34520a0..ec56a5ada 100644
--- a/raster/rt_core/rt_mapalgebra.c
+++ b/raster/rt_core/rt_mapalgebra.c
@@ -34,6 +34,25 @@
* rt_band_reclass()
******************************************************************************/
+static inline double
+rt_band_reclass_round_integer(rt_pixtype pixtype, double nv)
+{
+ 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:
+ return round(nv);
+ default:
+ return nv;
+ }
+}
+
/**
* Returns new band with values reclassified
*
@@ -94,134 +113,6 @@ rt_band_reclass(
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;
- }
- }
-
- /* Overflow checking */
- rt_util_dbl_trunc_warning(
- nodataval,
- checkvalint, checkvaluint,
- checkvalfloat, checkvaldouble,
- pixtype
- );
- }
- 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");
@@ -229,6 +120,8 @@ rt_band_reclass(
return 0;
}
rt_band_set_ownsdata_flag(band, 1); /* we DO own this data!!! */
+ rt_band_init_value(band, hasnodata ? nodataval : 0.0);
+
RASTER_DEBUGF(3, "rt_band_reclass: new band @ %p", band);
for (x = 0; x < width; x++) {
@@ -338,21 +231,7 @@ rt_band_reclass(
}
/* 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;
- }
+ nv = rt_band_reclass_round_integer(pixtype, nv);
RASTER_DEBUGF(4, "(%d, %d) ov: %f or: %f - %f nr: %f - %f nv: %f"
, x
@@ -378,6 +257,138 @@ rt_band_reclass(
return band;
}
+
+/*
+ * Pre-sorting the classpairs allows us to use bsearch
+ * on the array of classpairs in the pixel reclassification
+ * step later on.
+ */
+static int
+rt_classpair_cmp(const void *aptr, const void *bptr)
+{
+ struct rt_classpair_t *a = (struct rt_classpair_t*)aptr;
+ struct rt_classpair_t *b = (struct rt_classpair_t*)bptr;
+ if (a->src < b->src) return -1;
+ else if (a->src > b->src) return 1;
+ else return 0;
+}
+
+
+/**
+ * Returns new band with values reclassified
+ *
+ * @param srcband : the band who's values will be reclassified
+ * @param map : rt_reclassmap with src->dst mappings
+ * @param hasnodata : indicates if the user is supplying a nodata
+ * @param nodataval : user supplied nodata value for the new band
+ *
+ * @return a new rt_band or NULL on error
+ */
+rt_band
+rt_band_reclass_exact(
+ rt_band srcband, rt_reclassmap map,
+ uint32_t hasnodata, double nodataval
+) {
+ rt_band band = NULL;
+ uint32_t width = 0;
+ uint32_t height = 0;
+ int numval = 0;
+ int memsize = 0;
+ void *mem = NULL;
+ uint32_t src_hasnodata = 0;
+ rt_pixtype pixtype;
+
+ assert(NULL != srcband);
+ assert(NULL != map);
+
+ /* Validate that the map has consistent array lengths */
+ if (map->count == 0)
+ rterror("%s: reclassification map must contain at least one class pair", __func__);
+
+ /* Need a nodataval from somewhere */
+ src_hasnodata = rt_band_get_hasnodata_flag(srcband);
+ if (!(src_hasnodata || hasnodata))
+ rterror("%s: source band missing nodata value and nodata value not supplied", __func__);
+
+ /*
+ * Use the nodata value of the input in the case where user
+ * does not supply a nodataval
+ */
+ if (!hasnodata && src_hasnodata) {
+ rt_band_get_nodata(srcband, &nodataval);
+ hasnodata = src_hasnodata;
+ }
+
+ /* size of memory block to allocate */
+ width = rt_band_get_width(srcband);
+ height = rt_band_get_height(srcband);
+ numval = width * height;
+ pixtype = map->dsttype;
+ memsize = rt_pixtype_size(pixtype) * numval;
+ mem = (int *) rtalloc(memsize);
+ if (!mem)
+ rterror("%s: Could not allocate memory for band", __func__);
+
+ band = rt_band_new_inline(width, height, pixtype, hasnodata, nodataval, mem);
+ if (! band) {
+ rterror("%s: Could not add band to raster. Aborting", __func__);
+ rtdealloc(mem);
+ return NULL;
+ }
+ rt_band_set_ownsdata_flag(band, 1); /* we own this data */
+
+ /*
+ * apply nodata as baseline value and then only
+ * map in values that match our map
+ */
+ rt_band_init_value(band, nodataval);
+
+ /*
+ * Put the rt_classpairs in order of source value
+ * so we can bsearch them
+ */
+ qsort(map->pairs, map->count, sizeof(struct rt_classpair_t), rt_classpair_cmp);
+
+ for (uint32_t x = 0; x < width; x++) {
+ for (uint32_t y = 0; y < height; y++) {
+ int isnodata;
+ double nv = nodataval;
+ struct rt_classpair_t query = {0.0, 0.0};
+ struct rt_classpair_t *rslt;
+
+ /* get pixel, skip on error */
+ if (rt_band_get_pixel(srcband, x, y, &query.src, &isnodata) != ES_NONE)
+ continue;
+
+ /* output was already initialized to nodataval */
+ if (isnodata)
+ continue;
+ /*
+ * Look for pixel value in map.
+ * If not in map, skip, leaving the nodata value in place.
+ * Otherwise continue and replace with new value
+ */
+ rslt = bsearch(&query, map->pairs, map->count, sizeof(struct rt_classpair_t), rt_classpair_cmp);
+ if (!rslt)
+ continue;
+ else
+ nv = rslt->dst;
+
+ /* round the new value for integer pixel types */
+ nv = rt_band_reclass_round_integer(pixtype, nv);
+
+ if (rt_band_set_pixel(band, x, y, nv, NULL) != ES_NONE) {
+ rterror("%s: Could not assign value to new band", __func__);
+ rt_band_destroy(band);
+ rtdealloc(mem);
+ return NULL;
+ }
+ }
+ }
+
+ return band;
+}
+
/******************************************************************************
* rt_raster_iterator()
******************************************************************************/
diff --git a/raster/rt_core/rt_raster.c b/raster/rt_core/rt_raster.c
index 5ef18357b..d98591a86 100644
--- a/raster/rt_core/rt_raster.c
+++ b/raster/rt_core/rt_raster.c
@@ -499,12 +499,6 @@ rt_raster_generate_new_band(
int oldnumbands = 0;
int numbands = 0;
void * mem = NULL;
- int32_t checkvalint = 0;
- uint32_t checkvaluint = 0;
- double checkvaldouble = 0;
- float checkvalfloat = 0;
- int i;
-
assert(NULL != raster);
@@ -527,143 +521,27 @@ rt_raster_generate_new_band(
return -1;
}
- if (FLT_EQ(initialvalue, 0.0))
- memset(mem, 0, datasize);
- else {
- switch (pixtype)
- {
- case PT_1BB:
- {
- uint8_t *ptr = mem;
- uint8_t clamped_initval = rt_util_clamp_to_1BB(initialvalue);
- 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(initialvalue);
- 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(initialvalue);
- 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(initialvalue);
- 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(initialvalue);
- 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(initialvalue);
- 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(initialvalue);
- 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(initialvalue);
- 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(initialvalue);
- 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(initialvalue);
- 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] = initialvalue;
- checkvaldouble = ptr[0];
- break;
- }
- default:
- {
- rterror("rt_raster_generate_new_band: Unknown pixeltype %d", pixtype);
- rtdealloc(mem);
- return -1;
- }
- }
- }
-
- /* Overflow checking */
- rt_util_dbl_trunc_warning(
- initialvalue,
- checkvalint, checkvaluint,
- checkvalfloat, checkvaldouble,
- pixtype
- );
-
band = rt_band_new_inline(width, height, pixtype, hasnodata, nodatavalue, mem);
if (! band) {
rterror("rt_raster_generate_new_band: Could not add band to raster. Aborting");
rtdealloc(mem);
return -1;
}
- rt_band_set_ownsdata_flag(band, 1); /* we DO own this data!!! */
- index = rt_raster_add_band(raster, band, index);
- numbands = rt_raster_get_num_bands(raster);
- if (numbands == oldnumbands || index == -1) {
- rterror("rt_raster_generate_new_band: Could not add band to raster. Aborting");
- rt_band_destroy(band);
- }
- /* set isnodata if hasnodata = TRUE and initial value = nodatavalue */
- if (hasnodata && FLT_EQ(initialvalue, nodatavalue))
- rt_band_set_isnodata_flag(band, 1);
+ rt_band_init_value(band, initialvalue);
+
+ rt_band_set_ownsdata_flag(band, 1); /* we DO own this data!!! */
+ index = rt_raster_add_band(raster, band, index);
+ numbands = rt_raster_get_num_bands(raster);
+ if (numbands == oldnumbands || index == -1) {
+ rterror("rt_raster_generate_new_band: Could not add band to raster. Aborting");
+ rt_band_destroy(band);
+ return -1;
+ }
+
+ /* set isnodata if hasnodata = TRUE and initial value = nodatavalue */
+ if (hasnodata && FLT_EQ(initialvalue, nodatavalue))
+ rt_band_set_isnodata_flag(band, 1);
return index;
}
diff --git a/raster/rt_pg/rtpg_mapalgebra.c b/raster/rt_pg/rtpg_mapalgebra.c
index d898c63af..c1cac17c3 100644
--- a/raster/rt_pg/rtpg_mapalgebra.c
+++ b/raster/rt_pg/rtpg_mapalgebra.c
@@ -59,6 +59,7 @@ Datum RASTER_clip(PG_FUNCTION_ARGS);
/* reclassify specified bands of a raster */
Datum RASTER_reclass(PG_FUNCTION_ARGS);
+Datum RASTER_reclass_exact(PG_FUNCTION_ARGS);
/* apply colormap to specified band of a raster */
Datum RASTER_colorMap(PG_FUNCTION_ARGS);
@@ -3998,6 +3999,143 @@ Datum RASTER_reclass(PG_FUNCTION_ARGS) {
PG_RETURN_POINTER(pgrtn);
}
+
+/**
+ * ST_ReclassExact(rast raster,
+ * inputvalues float8[],
+ * outputvalues float8[],
+ * bandnumber integer DEFAULT 1,
+ * outputpixeltype text DEFAULT '32BF',
+ * nodatavalue float8 DEFAULT NULL
+ * ) RETURNS raster
+ */
+PG_FUNCTION_INFO_V1(RASTER_reclass_exact);
+Datum RASTER_reclass_exact(PG_FUNCTION_ARGS) {
+ rt_pgraster *pgraster = NULL;
+ rt_pgraster *pgrtn = NULL;
+ rt_raster raster = NULL;
+ rt_band band = NULL;
+ rt_band newband = NULL;
+ uint32_t bandnum = 0, numbands = 0;
+ ArrayType *arraysrc, *arraydst;
+ bool hasnodata = false;
+ double nodataval = 0.0;
+ text *pixeltype = NULL;
+ uint32_t szsrc, szdst;
+ ArrayIterator itersrc, iterdst;
+ Datum dsrc, ddst;
+ bool nullsrc, nulldst;
+ rt_reclassmap reclassmap;
+ rt_pixtype pixtypsrc, pixtypdst;
+
+ /* Check null on all arguments since function is not strict */
+ if (PG_ARGISNULL(0) || PG_ARGISNULL(1) || PG_ARGISNULL(2) || PG_ARGISNULL(3) || PG_ARGISNULL(4))
+ PG_RETURN_NULL();
+
+ /* Read SQL arguments */
+ pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ arraysrc = PG_GETARG_ARRAYTYPE_P(1);
+ arraydst = PG_GETARG_ARRAYTYPE_P(2);
+ bandnum = PG_GETARG_INT32(3); /* One-based band number */
+ pixeltype = PG_GETARG_TEXT_P(4);
+ if (!PG_ARGISNULL(5)) {
+ hasnodata = true;
+ nodataval = PG_GETARG_FLOAT8(5);
+ }
+
+ /* Check arrays have same size */
+ szsrc = ArrayGetNItems(ARR_NDIM(arraysrc), ARR_DIMS(arraysrc));
+ szdst = ArrayGetNItems(ARR_NDIM(arraydst), ARR_DIMS(arraydst));
+ if (szsrc != szdst)
+ elog(ERROR, "array lengths must be the same");
+
+ /* Read the raster */
+ raster = rt_raster_deserialize(pgraster, FALSE);
+ if (!raster) {
+ PG_FREE_IF_COPY(pgraster, 0);
+ elog(ERROR, "%s: Could not deserialize raster", __func__);
+ }
+
+ /* Check that this raster has bands we can work with */
+ numbands = rt_raster_get_num_bands(raster);
+ if (numbands < 1)
+ elog(ERROR, "Raster has no bands");
+ if (bandnum < 1 || bandnum > numbands)
+ elog(ERROR, "Invalid band index %d, input raster has %d bands. Band indexes are one-based.",
+ bandnum, numbands);
+
+ /* Read the band */
+ band = rt_raster_get_band(raster, bandnum-1);
+ if (!band)
+ elog(ERROR, "Could not find raster band of index %d", bandnum);
+
+ /* Get array element types */
+ pixtypsrc = rt_band_get_pixtype(band);
+ pixtypdst = rt_pixtype_index_from_name(text_to_cstring(pixeltype));
+
+ if (pixtypdst == PT_END)
+ elog(ERROR, "Unknown output pixel type '%s'", text_to_cstring(pixeltype));
+
+ /* Error out on unreadable pixeltype */
+ if (pixtypsrc == PT_END)
+ elog(ERROR, "Unsupported pixtype");
+
+ /* Create the rt_reclassmap */
+ reclassmap = palloc(sizeof(struct rt_reclassmap_t));
+ reclassmap->count = 0;
+ reclassmap->srctype = pixtypsrc;
+ reclassmap->dsttype = pixtypdst;
+ reclassmap->pairs = palloc(sizeof(struct rt_classpair_t) * szdst);
+
+ if (!hasnodata)
+ nodataval = rt_pixtype_get_min_value(pixtypdst);
+
+ /* Build up rt_reclassmap.pairs from arrays */
+ itersrc = array_create_iterator(arraysrc, 0, NULL);
+ iterdst = array_create_iterator(arraydst, 0, NULL);
+ while(array_iterate(itersrc, &dsrc, &nullsrc) &&
+ array_iterate(iterdst, &ddst, &nulldst))
+ {
+ double valsrc = nullsrc ? nodataval : (double)DatumGetFloat8(dsrc);
+ double valdst = nulldst ? nodataval : (double)DatumGetFloat8(ddst);
+ Assert(szdst > reclassmap->count);
+ reclassmap->pairs[reclassmap->count].src = valsrc;
+ reclassmap->pairs[reclassmap->count].dst = valdst;
+ reclassmap->count++;
+ }
+ array_free_iterator(itersrc);
+ array_free_iterator(iterdst);
+
+ /* Carry out reclassification */
+ newband = rt_band_reclass_exact(band, reclassmap, hasnodata, nodataval);
+ /* Clean up finished map */
+ pfree(reclassmap->pairs);
+ pfree(reclassmap);
+ if (!newband)
+ elog(ERROR, "Band reclassification failed");
+
+ /* replace old band with new band */
+ if (rt_raster_replace_band(raster, newband, bandnum-1) == NULL) {
+ rt_band_destroy(newband);
+ rt_raster_destroy(raster);
+ PG_FREE_IF_COPY(pgraster, 0);
+ elog(ERROR, "Could not replace raster band of index %d with reclassified band", bandnum);
+ }
+
+ pgrtn = rt_raster_serialize(raster);
+ rt_raster_destroy(raster);
+ PG_FREE_IF_COPY(pgraster, 0);
+ if (!pgrtn)
+ PG_RETURN_NULL();
+
+ SET_VARSIZE(pgrtn, pgrtn->size);
+ PG_RETURN_POINTER(pgrtn);
+}
+
+
+
+
+
/* ---------------------------------------------------------------- */
/* apply colormap to specified band of a raster */
/* ---------------------------------------------------------------- */
diff --git a/raster/rt_pg/rtpostgis.sql.in b/raster/rt_pg/rtpostgis.sql.in
index c4f897484..9760723d2 100644
--- a/raster/rt_pg/rtpostgis.sql.in
+++ b/raster/rt_pg/rtpostgis.sql.in
@@ -1325,6 +1325,16 @@ CREATE OR REPLACE FUNCTION _st_reclass(rast raster, VARIADIC reclassargset recla
AS 'MODULE_PATHNAME', 'RASTER_reclass'
LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE;
+CREATE OR REPLACE FUNCTION ST_ReclassExact(rast raster,
+ inputvalues float8[], outputvalues float8[],
+ bandnumber integer DEFAULT 1,
+ outputpixeltype text DEFAULT '32BF',
+ nodatavalue float8 DEFAULT NULL
+ )
+ RETURNS raster
+ AS 'MODULE_PATHNAME', 'RASTER_reclass_exact'
+ LANGUAGE 'c' IMMUTABLE PARALLEL SAFE;
+
CREATE OR REPLACE FUNCTION st_reclass(rast raster, VARIADIC reclassargset reclassarg[])
RETURNS raster
AS $$
diff --git a/raster/test/regress/rt_reclassexact.sql b/raster/test/regress/rt_reclassexact.sql
new file mode 100644
index 000000000..8002b88bd
--- /dev/null
+++ b/raster/test/regress/rt_reclassexact.sql
@@ -0,0 +1,94 @@
+DROP TABLE IF EXISTS reclassexact;
+CREATE TABLE reclassexact (
+ id integer,
+ rast raster
+);
+
+INSERT INTO reclassexact (id, rast)
+SELECT 1, ST_SetValues(
+ ST_AddBand(
+ ST_MakeEmptyRaster(
+ 2, -- width in pixels
+ 2, -- height in pixels
+ 0, -- upper-left x-coordinate
+ 0, -- upper-left y-coordinate
+ 1, -- pixel size in x-direction
+ -1, -- pixel size in y-direction (negative for north-up)
+ 0, -- skew in x-direction
+ 0, -- skew in y-direction
+ 4326 -- SRID (e.g., WGS 84)
+ ),
+ '32BUI'::text, -- pixel type (e.g., '32BF' for float, '8BUI' for unsigned 8-bit int)
+ 0.0, -- initial value for the band (e.g., 0.0 or a no-data value)
+ -99 -- nodatavalue
+ ),
+ 1, -- band number (usually 1 for single-band rasters)
+ 1, -- x origin for setting values (usually 1)
+ 1, -- y origin for setting values (usually 1)
+ ARRAY[
+ ARRAY[1, 2],
+ ARRAY[3, 4]
+ ]::double precision[][] -- 2D array of values
+ );
+
+-- Reclass and resort input arrays
+WITH rc AS (
+ SELECT ST_ReclassExact(
+ rast, -- input raster
+ ARRAY[4,3,2,1], -- input map
+ ARRAY[14,13,12,11], -- output map
+ 1, -- band number to remap
+ '32BUI' -- output raster pixtype
+ ) AS rast
+ FROM reclassexact
+ WHERE id = 1
+ )
+SELECT 'rce-1', (ST_DumpValues(rc.rast)).*
+FROM rc;
+
+-- Reclass to a new pixel type
+WITH rc AS (
+ SELECT ST_ReclassExact(
+ rast, -- input raster
+ ARRAY[4,3,2,1], -- input map
+ ARRAY[14,13,12,11], -- output map
+ 1, -- band number to remap
+ '32BF' -- output raster pixtype
+ ) AS rast
+ FROM reclassexact
+ WHERE id = 1
+ )
+SELECT 'rce-2', (ST_DumpValues(rc.rast)).*
+FROM rc;
+
+-- Mismatched array sizes
+WITH rc AS (
+ SELECT ST_ReclassExact(
+ rast, -- input raster
+ ARRAY[4,3,2,1,0], -- input map
+ ARRAY[14,13,12,11], -- output map
+ 1, -- band number to remap
+ '32BF' -- output raster pixtype
+ ) AS rast
+ FROM reclassexact
+ WHERE id = 1
+ )
+SELECT 'err-1', (ST_DumpValues(rc.rast)).*
+FROM rc;
+
+-- Bogus pixel type
+WITH rc AS (
+ SELECT ST_ReclassExact(
+ rast, -- input raster
+ ARRAY[4,3,2,1], -- input map
+ ARRAY[14,13,12,11], -- output map
+ 1, -- band number to remap
+ '32BFF' -- output raster pixtype
+ ) AS rast
+ FROM reclassexact
+ WHERE id = 1
+ )
+SELECT 'err-2', (ST_DumpValues(rc.rast)).*
+FROM rc;
+
+DROP TABLE IF EXISTS reclassexact;
diff --git a/raster/test/regress/rt_reclassexact_expected b/raster/test/regress/rt_reclassexact_expected
new file mode 100644
index 000000000..691c3931a
--- /dev/null
+++ b/raster/test/regress/rt_reclassexact_expected
@@ -0,0 +1,5 @@
+NOTICE: table "reclassexact" does not exist, skipping
+rce-1|1|{{11,12},{13,14}}
+rce-2|1|{{11,12},{13,14}}
+ERROR: array lengths must be the same
+ERROR: Unknown output pixel type '32BFF'
diff --git a/raster/test/regress/tests.mk b/raster/test/regress/tests.mk
index ac17fda9e..9d75e50e1 100644
--- a/raster/test/regress/tests.mk
+++ b/raster/test/regress/tests.mk
@@ -84,6 +84,7 @@ RASTER_TEST_UTILITY = \
$(top_srcdir)/raster/test/regress/rt_asjpeg \
$(top_srcdir)/raster/test/regress/rt_aspng \
$(top_srcdir)/raster/test/regress/rt_reclass \
+ $(top_srcdir)/raster/test/regress/rt_reclassexact \
$(top_srcdir)/raster/test/regress/rt_gdalwarp \
$(top_srcdir)/raster/test/regress/rt_gdalcontour \
$(top_srcdir)/raster/test/regress/rt_asraster \
-----------------------------------------------------------------------
Summary of changes:
NEWS | 2 +
doc/reference_raster.xml | 123 ++++++++++-
raster/rt_core/librtcore.h | 41 ++++
raster/rt_core/rt_band.c | 152 +++++++++++++-
raster/rt_core/rt_mapalgebra.c | 297 ++++++++++++++-------------
raster/rt_core/rt_raster.c | 150 ++------------
raster/rt_pg/rtpg_mapalgebra.c | 138 +++++++++++++
raster/rt_pg/rtpostgis.sql.in | 10 +
raster/test/regress/rt_reclassexact.sql | 94 +++++++++
raster/test/regress/rt_reclassexact_expected | 5 +
raster/test/regress/tests.mk | 1 +
11 files changed, 727 insertions(+), 286 deletions(-)
create mode 100644 raster/test/regress/rt_reclassexact.sql
create mode 100644 raster/test/regress/rt_reclassexact_expected
hooks/post-receive
--
PostGIS
More information about the postgis-tickets
mailing list