“Raster” data in PostGIS?

PostGIS does not support raster data as yet. Hopefully there should be raster support soon. Raster layers can be visualised from a local or remote server using geoserver. However having raster layers integrated into a PostGIS data base would be extremely useful for tasks such as overlaying herbarium records with continuous extrapolated climate data or elevation.

I have previously experimented with importing the points of a raster layer. This works quite well. It can be used for simple spatial overlays using proximity operations. However the visual impression produced by the points is not very satisfactory and a query using the contains operator would be more efficient.

QGIS has a plugin for graticule generation. This forms a grid of individual polygons. So if a grid is produced using the graticule plugin of QGIS and then imported into PostGIS using SPIT it is possible to partially mimic a raster layer, although with some computational overhead.

The document mnp-workshop3.pdf showed how worldclim data (http://www.worldclim.org/) can be imported into GRASS and R.

If R is run from GRASS the following lines can be used to import the appropriate raster layers at a chosen resolution into R.

Most of the lines of the script are concerned with calculations on the basic worldclim data. The dataset consists of 36 grids of monthly mean minimum temperature, maximum temperature and precipitation. Some mean annual values may be derived for simplicity.

a<-c(paste(“tmin”,1:12,sep=””),paste(“tmax”,1:12,sep=””),paste(“prec”,1:12,sep=””) ,”alt”)
d<-readRAST6(a)
d@data[,1:24]<-lapply(d@data[,1:24],function(x)x/10)
fullgrid(d)<-FALSE

f<-function(x)mean(x[1:12])
d@data[[“meanmim”]]<-apply(d@data,1,f)

f<-function(x)mean(x[13:24])
d@data[[“meanmax”]]<-apply(d@data,1,f)

f<-function(x)mean(x[1:24])
d@data[[“meantemp”]]<-apply(d@data,1,f)

f<-function(x)sum(x[25:36])
d@data[[“totprec”]]<-apply(d@data,1,f)

The trick for importing the layers to PostGIS from R is shown below. If the data is in R as a spatialPixelDataFrame the coordinates can be added to the data to form a table.

dd<-data.frame(coordinates(d),d@data)

This table can then be saved to a PostGIS database using ODBC.

library(RODBC)
con<-odbcConnect(“mydb”)
odbcQuery(con,”SET search_path =chiapas, pg_catalog;”)

sqlSave(con,dd,”worldclim”)

###########################################

Now either in pgadmin or running the queries directly through R by placing them inside an odbcQuery function.

ALTER TABLE chiapas.worldclim add column gid serial NOT NULL;
ALTER TABLE chiapas.worldclim ADD CONSTRAINT worldclim_pkey PRIMARY KEY(gid);
ALTER TABLE chiapas.worldclim ADD COLUMN the_geom geometry;

UPDATE chiapas.worldclim
SET the_geom = PointFromText(‘POINT(‘ || s1 || ‘ ‘ || s2 || ‘)’,4326);

CREATE INDEX
grat30arc_idx_geo on grat30arc USING GIST(the_geom GIST_GEOMETRY_OPS);

CREATE INDEX
worldclim_idx_geo on chiapas.worldclim USING GIST(the_geom GIST_GEOMETRY_OPS);

If a polygon graticule called grat30arc covering the same region with a 30 arch second grid has been produced using the graticule plugin in Qgis and imported into the public schema then the average values of the points within each grid square can then be calculated for the graticule and a new table produced using this as a model query.

CREATE TABLE chiapas.worldclim2 with oids as
SELECT a.*, g.totprec, g.meanmin, g.meanmax, g.elev from
grat30arc a,

(SELECT AVG(c.totprec) as totprec,
AVG(c.meanmim) as meanmin,
AVG(c.meanmax) as meanmax,
AVG(c.alt) as elev,d.gid
FROM chiapas.worldclim c, grat30arc d
WHERE c.the_geom && d.the_geom
AND contains(d.the_geom, c.the_geom) group by d.gid) g

WHERE a.gid=g.gid;

The result is very similar to the visualisation of an arcgrid in arcview when it is viewed in QGis over a shaded relief map. In this case the pixels are approximately 1 km x 1 km (30 arc seconds or 0.008333333 degrees) and the elevation lines up well with the contours derived from a 1:250,000 map.

Spatial indices on the herbarium table and the worldclim2 layer with the following

CREATE INDEX
herb_idx_geo on herb.herbario USING GIST(the_geom GIST_GEOMETRY_OPS);

CREATE INDEX
worldclim2_idx_geo on chiapas.worldclim2 USING GIST(the_geom GIST_GEOMETRY_OPS);
Then the following query runs in less than four seconds to provide 6,000 herbarium records with extrapolated climate data to the nearest km.

SELECT h.*, w.totprec, w.elev,w.meanmin,w.meanmax
FROM herb.herbario h, chiapas.worldclim2 w
WHERE h.the_geom && w.the_geom
AND contains(w.the_geom, h.the_geom);

Or

select h.*, w.totprec, w.elev,w.meanmin,w.meanmax
from herb.herbario h, chiapas.worldclim2 w
where ST_contains(w.the_geom, h.the_geom);