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

[Bborie Park]

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
>

-- 
Bborie Park
Programmer
Center for Vectorborne Diseases
UC Davis
530-752-8380
bkpark at ucdavis.edu