[postgis-users] raster map algebra question(s)

Paragon Corporation lr at pcorp.us
Thu Dec 22 17:05:43 PST 2011


Steve,

Yah I was assuming you had a structure similar to what Bborie had, but with
the further simplistic assumption that they weren't chunked.  If they are
chunked you'd do something like this insteady

So you might have rast, observation_date and the query would be 

 SELECT ST_Union(rast, 'MEAN' )  OVER (PARTITION BY ST_ConvexHull(rast) ROWS
BETWEEN 1 PRECEDING AND CURRENT
 ROW ORDER BY temp_date)  As newrast, grid_num, temp_date
 FROM temperatures;

-- NOTE that you can easily turn this to a 7 day moving average by changing
the PRECEDING etc.

I would go with storing all your like data in a single table with
observation_date similar to what Bborie has described.
FWIW:  The approach that Bborie has below is pretty much what 

ST_Union(rast,'MEAN')  does except ST_Union can do it on more than two
rasters at a time.

Bborie,

I am curious if you get the same answers using ST_Union as what you are
doing. Seems like a good test to stress test ST_Union MEAN.  

Is your data freely downloadable from somewhere.  I was going to pull some
from weather.gov to experiment with as I suspect moving averages will be a
common use case and I suspect no one has tried raster in combination with
window agg constructions. In theory it should work because all aggregates
are supposed to work in a window agg construct without programming them to
do so specifically.  At least as I recall all the sql aggs I've custom built
have worked fine in window constructs.

Hope that helps,
Regina
http://www.postgis.us



> -----Original Message-----
> From: postgis-users-bounces at postgis.refractions.net 
> [mailto:postgis-users-bounces at postgis.refractions.net] On 
> Behalf Of Stephen Crawford
> Sent: Thursday, December 22, 2011 4:10 PM
> To: PostGIS Users Discussion
> Subject: Re: [postgis-users] raster map algebra question(s)
> 
> That's great.  Thank You very much.
> 
> On 12/22/2011 3:53 PM, Bborie Park wrote:
> > Steve,
> >
> > I actually use similar datasets for my work.  I've structured my 
> > tables where a year's worth of rasters sit in one table.  So, 
> > tmin_2009, tmin_2010, tmin_2011.  Same goes with different 
> variables, 
> > tmax_2009, tmax_2010, tmax_2011.
> >
> > Each table has a column for the observation_date with the table 
> > structure looking like:
> >
> > rid serial
> > observation_date date
> > rast raster
> >
> > Each raster when loaded in is broken into 100x100 pixel 
> tiles.  One of 
> > the calculations I've been testing is the average temperature for a 
> > half month.  The query I use is:
> >
> > /* using ST_MapAlgebraExpr */
> > SELECT sum(sum) / sum(count)
> > FROM (
> > SELECT
> >     (st_summarystats(rast)).*
> > FROM (
> > SELECT
> >     ST_MapAlgebraExpr(n.tile, x.tile, '(rast1 + rast2) / 
> 2.', '32BF',
> > 'INTERSECTION') AS rast
> > FROM tmin_2010 n
> > JOIN tmax_2010 x
> >     ON n.observation_date = x.observation_date
> >         AND ST_Contains(ST_ConvexHull(n.tile), 
> ST_ConvexHull(x.tile))
> >         AND ST_Contains(ST_ConvexHull(x.tile), 
> ST_ConvexHull(n.tile)) 
> > WHERE n.observation_date BETWEEN '2010-09-01' AND '2010-09-15'
> > ) foo
> > ) bar;
> >
> > /* using ST_MapAlgebraFct */
> > CREATE OR REPLACE FUNCTION mapalgebra_avg(rast1 double precision,
> > rast2 double precision, VARIADIC args text[])
> >     RETURNS double precision
> >     AS $$
> >     DECLARE
> >     BEGIN
> >         IF rast1 IS NOT NULL AND rast2 IS NOT NULL THEN
> >             RETURN (rast1 + rast2) / 2.;
> >         END IF;
> >         RETURN NULL;
> >     END;
> >     $$ LANGUAGE 'plpgsql' IMMUTABLE;
> >
> > SELECT sum(sum) / sum(count)
> > FROM (
> > SELECT
> >     (st_summarystats(rast)).*
> > FROM (
> > SELECT
> >     ST_MapAlgebraFct(n.tile, x.tile, 'mapalgebra_avg(double 
> precision, 
> > double precision, text[])'::regprocedure, '32BF', 
> 'INTERSECTION') AS 
> > rast FROM tmin_2010 n JOIN tmax_2010 x
> >     ON n.observation_date = x.observation_date
> >         AND ST_Contains(ST_ConvexHull(n.tile), 
> ST_ConvexHull(x.tile))
> >         AND ST_Contains(ST_ConvexHull(x.tile), 
> ST_ConvexHull(n.tile)) 
> > WHERE n.observation_date BETWEEN '2010-09-01' AND '2010-09-15'
> > ) foo
> > ) bar;
> >
> > Hopefully, that gives you an idea of how to use ST_MapAlgebra.
> >
> > -bborie
> >
> > On 12/22/2011 12:00 PM, Stephen Crawford wrote:
> >> Hi Regina,
> >>
> >> Actually it's more like I have 7 raster tables, one for 
> each day, and 
> >> i would like to create a new raster as the average of 
> those seven days.
> >>
> >> AT this point I'm just trying to determine the structure and logic 
> >> for my db. Since rasters are new to me, I also have 
> fundamental questions.
> >> The data are always over the same domain, a number of 
> daily weather 
> >> variables (temp, rainfall, humidity, etc). Is it best to 
> break these 
> >> out by date and variable as I am now, with tables such as 
> >> "temp_12-22-2011", "temp_12-21-2011", "precip_12-22-2011", 
> >> "precip_12-21-2011", etc......or should I put all the vars for one 
> >> date in one table? In my my test db I have them by date, 
> by variable, 
> >> so the query I want to try first is the one I mention above.
> >>
> >> Thanks,
> >> Steve
> >>
> >> On 12/22/2011 2:44 PM, Paragon Corporation wrote:
> >>> Steve,
> >>>
> >>> Depends how your data is setup -- I'm assuming you are 
> looking for a 
> >>> moving average?
> >>>
> >>> If for example you have a raster record for each date in the same 
> >>> table, with a field denoting the date, then the best bet 
> is probably 
> >>> use ST_Union with the optional MEAN expression.
> >>>
> >>> I haven't stress tested that yet.
> >>>
> >>> So say you have a table
> >>>
> >>> Rast, temp_date
> >>>
> >>> Then you can do a
> >>>
> >>> SELECT ST_Union(rast, 'MEAN' ) OVER (ROWS BETWEEN 1 PRECEDING AND 
> >>> CURRENT ROW ORDER BY temp_date) As newrast, temp_date FROM 
> >>> temperatures;
> >>>
> >>>
> >>> Haven't tried using with window aggregates , but the 
> above should work.
> >>> Also note the ROWS with number ranges was introduced in 
> PostgreSQL 
> >>> 9.0 so that particular syntax won't work with 8.4.
> >>>
> >>> Hope that helps,
> >>> Regina
> >>>
> >>> -----Original Message-----
> >>> From: postgis-users-bounces at postgis.refractions.net
> >>> [mailto:postgis-users-bounces at postgis.refractions.net] On 
> Behalf Of 
> >>> Stephen Crawford
> >>> Sent: Thursday, December 22, 2011 2:23 PM
> >>> To: PostGIS Users Discussion
> >>> Subject: [postgis-users] raster map algebra question(s)
> >>>
> >>> Thanks to all who helped me on my previous raster output 
> questions; 
> >>> I now have a better handle on what I can (and can't) do.
> >>>
> >>> How can I do map algebra across raster tables? I have a 
> db of raster 
> >>> tables (all the same extend) with weather variables by 
> date. I have 
> >>> a number of ways i want to use these data, but in the 
> initial case I 
> >>> would like to, as an example, create a new raster table where the 
> >>> data value at each cell is the average of the previous 
> day's values. 
> >>> Any thoughts?
> >>>
> >>> Thanks,
> >>> Steve
> >>>
> >>> --
> >>> Stephen Crawford
> >>> Center for Environmental Informatics The Pennsylvania State 
> >>> University
> >>>
> >>>
> >>>
> >>> _______________________________________________
> >>> postgis-users mailing list
> >>> postgis-users at postgis.refractions.net
> >>> http://postgis.refractions.net/mailman/listinfo/postgis-users
> >>>
> >>>
> >>> _______________________________________________
> >>> postgis-users mailing list
> >>> postgis-users at postgis.refractions.net
> >>> http://postgis.refractions.net/mailman/listinfo/postgis-users
> >>
> >
> 
> --
> Stephen Crawford
> Center for Environmental Informatics
> The Pennsylvania State University
> src176 at psu.edu
> 814.865.9905
> 
> 
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
> 





More information about the postgis-users mailing list