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

[Stephen Crawford]

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