[postgis-users] Aggregating rasters by adding and other confusions
David M. Kaplan
david.kaplan at ird.fr
Mon Mar 12 04:48:13 PDT 2018
Hi,
I did some tests to try to figure this out and it seems like most of
the confusion was due to integer mathematics (though I am not sure that
the current behavior is the "correct" behavior). First, I created some
test rasters:
CREATE TEMP TABLE test_rasters ASSELECT 1 AS rid, ST_SetValues(
ST_AddBand( ST_MakeEmptyRaster(2, 2, 0, 0, 1, -1, 0, 0, 0),
1, '8BUI', 0 ), 1, 1, 1, ARRAY[[1,2],
[0,0]]::double precision[][] ) AS rastUNIONSELECT 2 AS rid,
ST_SetValues( ST_AddBand( ST_MakeEmptyRaster(2, 2, 0, 0,
1, -1, 0, 0, 0), 1, '8BUI', 0 ), 1, 1, 1,
ARRAY[[1,0], [3,0]]::double precision[][] ) AS rastUNIONSELECT 3
AS rid, ST_SetValues( ST_AddBand( ST_MakeEmptyRaster(2,
2, 0, 0, 1, -1, 0, 0, 0), 1, '8BUI', 0 ),
1, 1, 1, ARRAY[[0,0], [3,0]]::double precision[][] ) AS rast;
This leads to the following rasters:
# SELECT rid, ST_DumpValues(rast) FROM test_rasters; rid
| st_dumpvalues -----+--------------------- 3 |
(1,"{{0,0},{3,0}}") 1 | (1,"{{1,2},{0,0}}") 2 |
(1,"{{1,0},{3,0}}")(3 rows)
Then I did the following calculations on those rasters:
# SELECT ST_DumpValues(ST_Union(rast,1,'COUNT')) FROM
test_rasters; st_dumpvalues ---------------------
(1,"{{3,3},{3,3}}")
# SELECT ST_DumpValues(ST_Union(rast,1,'SUM')) FROM
test_rasters; st_dumpvalues ----------------------
-- (1,"{{2,2},{6,NULL}}")(1 row)
# SELECT ST_DumpValues(ST_Union(rast,1,'MEAN')) FROM
test_rasters; st_dumpvalues ----------------------
-- (1,"{{1,1},{2,NULL}}")(1 row)
# SELECT
ST_DumpValues(ST_MapAlgebra(ST_Union(rast,1,'SUM'),'8BUI','[rast] /
3')) FROM test_rasters; st_dumpvalues ----------------------
-- (1,"{{1,1},{2,NULL}}")(1 row)
# SELECT
ST_DumpValues(ST_MapAlgebra(ST_Union(rast,1,'SUM'),'64BF','[rast] /
3')) FROM
test_rasters; st_dumpvalues
-----------------------------------------------------
--- (1,"{{0.666666666666667,0.666666666666667},{2,NULL}}")(1 row)
So ST_Union is counting the zero values, but when it calculates the
MEAN, it does integer mathematics that leads to unexpected results.
Whereas ST_MapAlgebra allows you to select the output raster pixel
type, ST_Union does not, so I don't think there is a good work around
for this. Is there a way to "cast" a raster from one pixel type to
another? I imagine that with appropriate calls to ST_DumpValues and
ST_SetValues one could make it happen.
One could argue whether or not the current behavior is correct. Based
on my experience with the standard avg() aggregate function, I expected
ST_Union to select an output pixel type that is appropriate for the
calculation being carried out, i.e., that it would output a '32BF' or
'64BF' raster as integer means are a bit unusual. Perhaps ST_Union
should have an output pixel type argument or by default return '64BF'?
On a related note, I noticed that ST_SetBandNoDataValue has somewhat
unexpected behavior when you set the no data value to something that is
outside the range of the pixel type, basically setting the no data
value to the max or min of the data range:
# SELECT (ST_BandMetadata(ST_SetBandNoDataValue(rast,1,-
99))).nodatavalue,# ST_DumpValues(ST_SetBandNoDataValue(rast,1,-
99)) AS values# FROM test_rasters WHERE rid=1; nodatavalue
| values -------------+------------------------
--- 0 | (1,"{{1,2},{NULL,NULL}}")(1 row)
# SELECT
(ST_BandMetadata(ST_SetBandNoDataValue(rast,1,99))).nodatavalue,#
ST_DumpValues(ST_SetBandNoDataValue(rast,1,99)) AS values# FROM
test_rasters WHERE rid=1; nodatavalue | values -----------
--+--------------------- 99 | (1,"{{1,2},{0,0}}")(1 row)
# SELECT
(ST_BandMetadata(ST_SetBandNoDataValue(rast,1,999))).nodatavalue,#
ST_DumpValues(ST_SetBandNoDataValue(rast,1,999)) AS values# FROM
test_rasters WHERE rid=1; nodatavalue | values -----------
--+--------------------- 255 | (1,"{{1,2},{0,0}}")(1 row)
I would have expected it to return an error when the demanded no data
value was outside the range. Should this behavior be changed?
Thanks,David
-----------------------------------------------------------------------
On Sat, 2018-03-10 at 12:00 -0800, postgis-users-request at lists.osgeo.or
g wrote:
> Message: 1
> Date: Fri, 9 Mar 2018 23:22:56 -0500
> From: "Regina Obe" <lr at pcorp.us>
> To: "'PostGIS Users Discussion'" <postgis-users at lists.osgeo.org>
> Subject: Re: [postgis-users] Aggregating rasters by adding and other
> confusions
> Message-ID: <000301d3b827$77f497b0$67ddc710$@pcorp.us>
> Content-Type: text/plain; charset="utf-8"
>
> David,
>
>
>
> I took a cursory glance at the union code we have in place.
>
> What it seems to do is two passes
>
>
>
> 1 pass does an ST_Union using 'COUNT' (so in your case you'd get
> numbers between 0 and 3, 0 being no rasters considered having any
> data)
>
> 2nd pass does an ST_Union using 'SUM'
>
> And returns a new raster where each pixel is SUM/COUNT
>
>
>
> Basic code is here:
>
>
>
> http://postgis.net/docs/doxygen/2.4/da/dde/rtpg__mapalgebra_8c_a1d940
> 65e6cef5d5d61417b82b2cf4fb6.html#a1d94065e6cef5d5d61417b82b2cf4fb6
>
> (note it only computes a value if the first band (which I presume to
> be the count band) value > 0 and has no nodata value. )
>
> I don't know why a count band would have a nodata value aside from
> when it's 0
>
>
>
> What does COUNT return for you?
>
>
>
> Thanks,
>
> Regina
>
>
>
> From: postgis-users [mailto:postgis-users-bounces at lists.osgeo.org] On
> Behalf Of David M. Kaplan
> Sent: Tuesday, March 06, 2018 10:29 AM
> To: postgis-users at lists.osgeo.org
> Subject: Re: [postgis-users] Aggregating rasters by adding and other
> confusions
>
>
>
> Hi,
>
>
>
> Yes, this does seem to be the solution. For completeness, the
> following did the same thing as my aggregate function:
>
>
>
> SELECT ST_Union(rast,1,'SUM') FROM blah;
>
> This includes replacing 0's with NULL's - ST_Union had the same
> behavior at my aggregate function.
>
>
>
> Can you explain this? For my sum, this isn't a big deal, but
> ST_Union(... 'MEAN') seems to be treating 0's at NULL values, thereby
> removing them from the mean calculation:
>
>
>
> db=# SELECT (ST_SummaryStats(ST_Union(rast,1,'SUM'))).sum FROM blah;
> sum
> ------
> 2792
> (1 row)
>
> db=# SELECT (ST_SummaryStats(ST_Union(rast,1,'MEAN'))).sum FROM blah;
> sum
> ------
> 1346
> (1 row)
>
>
> The table "blah" has 3 rasters with no null values, so I would expect
> the sum of the mean to be 1/3 the sum of the sum, but clearly this is
> not the case.
>
>
>
> I tried playing around with the no data value for my rasters to see
> if this would change things. It does change things, but not in any
> way I can explain:
>
>
>
> fads=# SELECT
> (ST_SummaryStats(ST_Union(ST_SetBandNoDataValue(rast,1,0),1,'MEAN')))
> .sum FROM blah;
> sum
> ------
> 2025
> (1 row)
>
> fads=# SELECT
> (ST_SummaryStats(ST_Union(ST_SetBandNoDataValue(rast,1,-
> 999),1,'MEAN'))).sum FROM blah;
> sum
> ------
> 2025
> (1 row)
>
>
> Setting the no data value has no effect on the sum of the sum.
>
>
>
> Does any of this make sense?
>
>
>
> Thanks,
>
> David
>
>
>
>
>
> On Mon, 2018-03-05 at 12:00 -0800, postgis-users-request at lists.osgeo.
> org <mailto:postgis-users-request at lists.osgeo.org> wrote:
>
> Date: Mon, 5 Mar 2018 12:05:40 -0500
> From: "Regina Obe" <lr at pcorp.us <mailto:lr at pcorp.us> >
> To: "'PostGIS Users Discussion'" <postgis-users at lists.osgeo.org <mail
> to:postgis-users at lists.osgeo.org> >
> Subject: Re: [postgis-users] Aggregating rasters by adding and other
> confusions
> Message-ID: <001f01d3b4a4$30d1fa20$9275ee60$@pcorp.us <mailto:001f01d
> 3b4a4$30d1fa20$9275ee60$@pcorp.us> >
> Content-Type: text/plain; charset="utf-8"
>
> I think what you are looking for is ST_Union (.. SUM) note this has
> union types – FIRST, MIN, MAX, COUNT, SUM, MEAN, RANGE
>
>
>
> http://postgis.net/docs/manual-2.4/RT_ST_Union.html
>
>
>
>
>
>
>
> From: postgis-users [mailto:postgis-users-bounces at lists.osgeo.org] On
> Behalf Of David M. Kaplan
> Sent: Monday, March 05, 2018 10:03 AM
> To: postgis-users at lists.osgeo.org <mailto:postgis-users at lists.osgeo.o
> rg>
> Subject: [postgis-users] Aggregating rasters by adding and other
> confusions
>
>
>
> Hi,
>
>
>
> I have recently started working with the postgis raster
> functionality. In general, I have found this really useful and have
> been able to do some neat things fairly simply with this raster
> functionality. Nevertheless, there are a few basic things that I am
> confused about and I was hoping someone could give me a hand.
>
>
>
> (1) First of all, I have a table with a bunch of rasters that have
> the same extent, alignment, scale, etc. and I want to aggregate them
> together into a single raster using pixel-by-pixel addition. It seems
> like there should be a function to do this, but I can't find one. Is
> there an aggregate "ST_MapAlgebra" function?
>
>
>
> Given that I couldn't find one, I defined an aggregate function as
> follows:
>
>
>
> CREATE OR REPLACE FUNCTION AddRasters(r1 raster, r2 raster)
> RETURNS raster AS
> $BODY$
> SELECT ST_MapAlgebra($1,$2,'[rast1]+[rast2]');
> $BODY$ LANGUAGE 'sql' IMMUTABLE STRICT;
>
> CREATE AGGREGATE sum (raster)
> (
> sfunc = AddRasters,
> stype = raster
> );
>
>
> (2) This seems to work, but it has the unexpected behavior that it
> replaces 0 values with NULL. In my case, this is fine, but I am
> wondering why it does this? I can't find anything that indicates that
> it should be replacing zeros with NULL. Here is the metadata
> associated with one of my rasters (the others are similar):
>
>
>
> # SELECT ST_BandMetadata(rast), ST_Metadata(rast),
> ST_SummaryStats(rast) FROM blah;
> -[ RECORD 1 ]---+--------------------------------------------------
> -----
> st_bandmetadata | (16BUI,,f,)
> st_metadata | (-180,90,360,180,1,-1,0,0,4326,1)
> st_summarystats |
> (64800,417,0.00643518518518519,0.223617719977485,0,46)
>
>
> I have not defined a nodataval for these layers and the original
> layers have no NULL values.
>
>
>
> (3) Is there a postgis command to turn all the NULL values back into
> zeros?
>
>
>
> (4) I was also considering just defining the '+' operator for raster
> + raster to be pixel-by-pixel addition. Is there any reason that I
> wouldn't want to do this?
>
>
>
> (5) Finally, I have been visualizing results with QGIS using the DB
> Manager. However, I don't see how to select a row from a raster table
> and incorporate just that row into the canvas. Is there a way to do
> this?
>
>
>
> Thanks for the assistance,
>
> David
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20180312/515e2cbb/attachment.html>
More information about the postgis-users
mailing list