[raster] Add a series of ST_SetValues functions

[pracine]

The first one takes the values from another raster extent.

The second sets the extent to a single value.

The third sets the shape of a geometry to a single value.

[pracine]

The first and the second one have plpgsql prototypes in script/plpgsql/st_setvalues.sql

[Bborie Park]

The third variant (geometry-based) could make use of the geomval data type...

ST_SetValues(
	rast raster,
	nband integer,
	geomvalset geomval[]
) -> raster

The geomval array allows the user to set many pixels at once in one call to ST_SetValues.

[Bborie Park]

Another variant to add allows setting many pixels at once based upon a one or two dimension array...

The following would set the values for a pixel-line using the x, y coordinate as the left-most pixel to set. If the number of elements in valset exceeds (width - x), the remaining values are discarded.

ST_SetValues(
	rast raster,
	nband integer,
	x integer,
	y integer,
	valset double precision[]
) -> raster

The following would set the values for an area using the x, y coordinate as the upper-left pixel to set. Each inner array would be a pixel-line. If the number of elements in a inner array exceeds (width - x), the remaining values of that array are discarded. The same is true for the outer array where the number of elements exceeds (height - y).

ST_SetValues(
	rast raster,
	nband integer,
	x integer,
	y integer,
	valset double precision[][]
) -> raster

[Bborie Park]

Looking at the variant taking values from a second raster in the plpgsql prototype, I had a question. Is the second raster a simple data source since there is no validation of pixel types, SRID, width, height?

[Bborie Park]

Added variant #1 from raster/scripts/plpgsql/st_setvalues.sql in r10147. Added variant #4 using an array in r10146.

[pracine]

Replying to dustymugs:

Looking at the variant taking values from a second raster in the plpgsql prototype, I had a question. Is the second raster a simple data source since there is no validation of pixel types, SRID, width, height?

This variant works in world coordinate so if the SRID, the width, the height are different it should not matter. In the worst case nothing is copied. If the pixel type are different, the source value should be casted to the destination pixel type.

This variant will be, with ST_AddBand(from another raster), our only way to copy values from another raster.

[Bborie Park]

Are you sure about Variant 1 (from scripts/plpgsql) working in world coordinates? For example, with a skewed raster...

WITH foo AS (
	SELECT ST_SetValues(
		ST_AddBand(ST_MakeEmptyRaster(3, 3, 0, 0, 1, -1, 1, 1, 0), 1, '8BUI', 1, 0),
		1,
		1, 1,
		ARRAY[[1, 2, 3], [4, 5, 6], [7, 8, 9]]
	) AS rast
),
bar AS (
	SELECT ST_AddBand(ST_MakeEmptyRaster(3, 3, 1, -1, 1, -1, 0, 0, 0), 1, '8BUI', 0, 0) AS rast
)
SELECT
	(pap).x,
	(pap).y,
	(pap).val
FROM (
	SELECT
		ST_PixelAsPolygons(ST_SetValues(
			bar.rast, 1,
			NULL, NULL,
			NULL, NULL,
			foo.rast, 1,
			FALSE, FALSE
		)) AS pap
	FROM bar
	CROSS JOIN foo
) t

results

psql:test.sql:27: NOTICE:  aaa oldy=1, newheight=<NULL>
psql:test.sql:27: NOTICE:  bbb newx=1, newy=1
psql:test.sql:27: NOTICE:  ccc newwidth=3, newheight=3
psql:test.sql:27: NOTICE:  111 x2=0, y2=0
psql:test.sql:27: NOTICE:  222 newx=1, newy=1
psql:test.sql:27: NOTICE:  333 newwidth=2, newheight=2
psql:test.sql:27: NOTICE:  555 newx + newwidth - 1=2, newy + newheight - 1=2
psql:test.sql:27: NOTICE:  666 x=1, y=1
psql:test.sql:27: NOTICE:  666 x=1, y=2
psql:test.sql:27: NOTICE:  666 x=2, y=1
psql:test.sql:27: NOTICE:  666 x=2, y=2
 x | y | val 
---+---+-----
 1 | 1 |   5
 2 | 1 |   6
 3 | 1 |    
 1 | 2 |   8
 2 | 2 |   9
 3 | 2 |    
 1 | 3 |    
 2 | 3 |    
 3 | 3 |

That is completely incorrect especially when visualized in the attached image (flawed.png).

Another simpler example, different scales...

WITH foo AS (
	SELECT ST_SetValues(
		ST_AddBand(ST_MakeEmptyRaster(3, 3, 0, 0, 0.1, -0.1, 0, 0, 0), 1, '8BUI', 1, 0),
		1,
		1, 1,
		ARRAY[[1, 2, 3], [4, 5, 6], [7, 8, 9]]
	) AS rast
),
bar AS (
	SELECT ST_AddBand(ST_MakeEmptyRaster(3, 3, 0, 0, 1, -1, 0, 0, 0), 1, '8BUI', 0, 0) AS rast
)
SELECT
	(pap).x,
	(pap).y,
	(pap).val
FROM (
	SELECT
		ST_PixelAsPolygons(ST_SetValues(
			bar.rast, 1,
			NULL, NULL,
			NULL, NULL,
			foo.rast, 1,
			FALSE, FALSE
		)) AS pap
	FROM bar
	CROSS JOIN foo
) t

results

psql:test.sql:27: NOTICE:  aaa oldy=1, newheight=<NULL>
psql:test.sql:27: NOTICE:  bbb newx=1, newy=1
psql:test.sql:27: NOTICE:  ccc newwidth=3, newheight=3
psql:test.sql:27: NOTICE:  111 x2=1, y2=1
psql:test.sql:27: NOTICE:  222 newx=1, newy=1
psql:test.sql:27: NOTICE:  333 newwidth=3, newheight=3
psql:test.sql:27: NOTICE:  555 newx + newwidth - 1=3, newy + newheight - 1=3
psql:test.sql:27: NOTICE:  666 x=1, y=1
psql:test.sql:27: NOTICE:  666 x=1, y=2
psql:test.sql:27: NOTICE:  666 x=1, y=3
psql:test.sql:27: NOTICE:  666 x=2, y=1
psql:test.sql:27: NOTICE:  666 x=2, y=2
psql:test.sql:27: NOTICE:  666 x=2, y=3
psql:test.sql:27: NOTICE:  666 x=3, y=1
psql:test.sql:27: NOTICE:  666 x=3, y=2
psql:test.sql:27: NOTICE:  666 x=3, y=3
 x | y | val 
---+---+-----
 1 | 1 |   1
 2 | 1 |   2
 3 | 1 |   3
 1 | 2 |   4
 2 | 2 |   5
 3 | 2 |   6
 1 | 3 |   7
 2 | 3 |   8
 3 | 3 |   9

By the looks of it, the function is working in raster coordinates instead.

[pracine]

Right. Looking at the code it does not fetch values from rast2 using world coordinates. It should... That would make more sense no? Otherwise the function would support only well aligned rasters.

Maybe I wrote "Both ST_SetValues functions found in this file are ready to be being implemented in C" too fast :)

[Bborie Park]

If the rasters aren't aligned, there isn't much that can be done. A function that copies values from one raster to another in raster coordinates is still useful though as you could be stitching together a set of tiles into one large raster (union basically but a specific case).

[pracine]

Replying to dustymugs:

If the rasters aren't aligned, there isn't much that can be done.

Why? Just convert the first raster pixel x and y into world coordinates (the centroid acutally) and then those world coordinates into the second raster x and y coordinates to look for a value. If the world coordinates fall outside the second raster extent and keepsourcetnodata is false, just leave the destination value unchanged. If keepsourcetnodata is true, the destination becomes nodata. In other word, threat any out of bound coordinate as nodata. That's what we should always do.

[Bborie Park]

So treat pixels as points? It doesn't work if the pixels are treated as areas.

[pracine]

Why not? I agree that when the pixel size of the second raster is significantly smaller than the pixel size of the first, the copied value is not really representative and that the main goal of this function is to provide a fast way to copy raster to raster (avoiding the complexity of determining aggregated value for one pixel) and to be used by ST_MapAlgebra which does not work on unaligned rasters...

Taking the centroid does not change the result when the rasters are aligned and it make the function more flexible for other use cases. I would just document the fact that using ST_SetValues with raster of different pixel size can results in not so representative pixel values.

If you want to threat the pixel as a polygon, use ST_MapAlgebraFct with a special function converting the shape of the pixel to a polygon, intersecting with the second raster and offering a series of methods to select the pixel value: COUNT, DISTINCT_COUNT, COUNT_MOST_FREQUENT, COUNT_LEAST_FREQUENT, MOST_FREQUENT_VALUE, LEAST_FREQUENT_VALUE, MAXIMUM, MINIMUM, RANGE, SUM, MEAN, STDDEV, etc... See Objective FV.27 in the specs.

[bnordgren]

I smell scope creep. You're starting to describe a hybrid of ST_Resample and ST_Transform (but backwards). I think it would be better to allow ST_Transform and ST_Resample to accept a target raster, rather than allow this "pixel accessor" to accept a source raster with potentially a completely different geotransform and SRID. Keep a distinction between pixel accessor and geoprocessing in the API.

If the above methods of selecting the pixel value are desirable, I would vote for adding them to ST_Resample/ST_Transform.

[Bborie Park]

Variant #3 added in r10465.

[Bborie Park]

Closing ticket as the variant getting values from another raster is confusing, especially after receiving comments.

[pracine]

Can we create another ticket so we do not forget about a variant copying pixels from another band?

[Bborie Park]

Another ticket for that one variant would be great. This ticket is just too cluttered and too many elements.