Re: [postgis-users] Aggregating rasters by adding and other confusions

2018-03-12 Thread David M. Kaplan
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.667,0.667},{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-requ...@lists.osgeo.or
g wrote:
> Message: 1
> Date: Fri, 9 Mar 2018 23:22:56 -0500
> From: "Regina Obe" <l...@pcorp.us>
> To: "'PostGIS Users Discussion'" <postgis-users@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 i

Re: [postgis-users] Aggregating rasters by adding and other confusions

2018-03-09 Thread Regina Obe
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_a1d94065e6cef5d5d61417b82b2cf4fb6.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-boun...@lists.osgeo.org] On Behalf Of 
David M. Kaplan
Sent: Tuesday, March 06, 2018 10:29 AM
To: postgis-users@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-requ...@lists.osgeo.org 
<mailto:postgis-users-requ...@lists.osgeo.org>  wrote:

Date: Mon, 5 Mar 2018 12:05:40 -0500
From: "Regina Obe" <l...@pcorp.us <mailto:l...@pcorp.us> >
To: "'PostGIS Users Discussion'" <postgis-users@lists.osgeo.org 
<mailto:postgis-users@lists.osgeo.org> >
Subject: Re: [postgis-users] Aggregating rasters by adding and other
confusions
Message-ID: <001f01d3b4a4$30d1fa20$9275ee60$@pcorp.us 
<mailto:001f01d3b4a4$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-boun...@lists.osgeo.org] On Behalf Of 
David M. Kaplan
Sent: Monday, March 05, 2018 10:03 AM
To: postgis-users@lists.osgeo.org <mailto:postgis-users@lists.osgeo.org> 
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;

Re: [postgis-users] Aggregating rasters by adding and other confusions

2018-03-06 Thread David M. Kaplan
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-
requ...@lists.osgeo.org wrote:
> Date: Mon, 5 Mar 2018 12:05:40 -0500
> From: "Regina Obe" <l...@pcorp.us>
> To: "'PostGIS Users Discussion'" <postgis-users@lists.osgeo.org>
> Subject: Re: [postgis-users] Aggregating rasters by adding and other
> confusions
> Message-ID: <001f01d3b4a4$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-boun...@lists.osgeo.org] On Behalf 
> Of David M. Kaplan
> Sent: Monday, March 05, 2018 10:03 AM
> To: postgis-users@lists.osgeo.org
> 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___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Aggregating rasters by adding and other confusions

2018-03-05 Thread Regina Obe
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-boun...@lists.osgeo.org] On Behalf Of 
David M. Kaplan
Sent: Monday, March 05, 2018 10:03 AM
To: postgis-users@lists.osgeo.org
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

 

 

 

-- 

**

David M. Kaplan

Charge de Recherche 1

 

Institut de Recherche pour le Developpement (IRD)

UMR MARBEC (IRD/Ifremer/CNRS/UMII)

av. Jean Monnet

CS 30171

34203 Sete cedex

France

 

Email: david.kap...@ird.fr  

Phone: +33 (0)4 99 57 32 25

Fax: +33 (0)4 99 57 32 95

 

http://www.umr-marbec.fr/kaplan-david.html

http://www.davidmkaplan.fr/

**

___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users