R/PostGIS functions

by Mathieu Basille.

Here is a set of R functions that allow some interaction between R and PostGIS. All functions need the RODBC package. For every function, the first argument db is the connection identifier from RODBC:

library(RODBC)
db <- odbcConnect("connection_id")

pgAddColumn creates a new column with the name colname and the type coltype (integer by default) to the table table.

pgAddColumn <- function(db, table, colname, coltype = "integer")
### Mathieu Basille, basille@ase-research.org
### Last modified: 2011-10-30
{
    str <- paste("ALTER TABLE", table, "ADD COLUMN", colname,
        coltype, ";")
    sqlQuery(db, str)
}

pgComment replaces the comment of the table table with comment.

pgComment <- function(db, table, comment, type = c("TABLE", "VIEW"))
### Mathieu Basille, basille@ase-research.org
### Last modified: 2011-10-30
{
    str <- paste("COMMENT ON ", type, " ", table, " IS '", comment,
        "';", sep = "")
    sqlQuery(db, str)
}

pgAsDate converts the column date of a table table into a timestamp with timezone. The column date must be in the format 'YYYY-MM-DD HH:MM:SS TZ' (note the TZ, e.g. 'EST5EDT').

pgAsDate <- function(db, table, date)
### Mathieu Basille, basille@ase-research.org
### Last modified: 2011-10-30
{
    str <- paste("ALTER TABLE ", table, " ALTER COLUMN ", date,
        " TYPE timestamptz USING ", date, "::timestamptz;", sep = "")
    sqlQuery(db, str)
}

pgIntersectPts intersects endpoints (or startpoints, use type = "Start") of a multisteps steps with a raster rast. If create, it stores the result in the table table, otherwise send the table to R.

pgIntersectPts <- function(db, geom, raster, type = c("Point",
    "StartPoint", "EndPoint"), colname = "val", create = FALSE,
    table)
### Mathieu Basille, basille@ase-research.org
### Last modified: 2011-10-30
{
    type <- match.arg(type)
    geom <- switch(type, Point = geom, StartPoint = paste("(SELECT gid, ST_StartPoint(the_geom) AS the_geom FROM",
        geom, ")"), EndPoint = paste("(SELECT gid, ST_EndPoint(the_geom) AS the_geom FROM",
        geom, ")"))
    str <- paste("SELECT gid, (pix).val AS ", colname, " ", "FROM (SELECT pts.gid, ",
        "  ST_Intersection(rast, the_geom) AS pix ", "  FROM ",
        raster, ", ", geom, " AS pts ", "  WHERE ST_Intersects(rast, the_geom)",
        ") AS inter ORDER BY gid;", sep = "")
    if (create)
        str <- paste("CREATE TABLE ", table, " AS ", str)
    sqlQuery(db, str)
}

pgIntersectSts intersects steps of a multisteps steps with a raster rast. If type = "cont", it computes the mean value of the raster for each step and returns it a column colname; if type = "fac", it computes the total length of each class with names colnameN (with N from min to max values) for each step. If create, it stores the result in the table table, otherwise send the table to R.

pgIntersectSts <- function(db, steps, raster, type = c("cont",
    "fac"), colname = "val", min = 1, max, create = FALSE, table)
### Mathieu Basille, basille@ase-research.org
### Last modified: 2011-10-30
{
    if (type == "cont")
        str <- paste("SELECT gid, SUM(val * ST_Length(the_geom))/SUM(ST_Length(the_geom)) AS ", colname, " ",
        "FROM (",
        "  SELECT gid,",
        "    (rv).geom AS the_geom,",
        "    (rv).val",
        "  FROM (",
        "    SELECT sts.gid, ST_Intersection(rast, the_geom) AS rv",
        "    FROM ", raster, ",", steps, " AS sts",
        "    WHERE ST_Intersects(rast, the_geom)",
        "  ) AS interseg) as intermean ",
        "GROUP BY gid ",
        "ORDER BY gid;", sep = "")
    else if (type == "fac")
        str <- paste("SELECT * FROM crosstab(",
        "  'SELECT gid, val, SUM(ST_Length(the_geom))",
        "  FROM (",
        "    SELECT gid,",
        "      (rv).geom AS the_geom,",
	"      (rv).val",
        "    FROM (",
        "      SELECT sts.gid, ST_Intersection(rast, the_geom) AS rv",
        "      FROM ", raster, ", ", steps, " AS sts ",
	"      WHERE ST_Intersects(rast, the_geom)",
	"    ) AS interseg) AS interfac",
        "  GROUP BY gid, val",
        "  ORDER BY gid',",
        "  'SELECT generate_series(", min, ",", max, ")') ",
        "AS ct(gid int", paste(", ", colname, min:max, " numeric",
                               sep = "", collapse = ""), ") ",
        "ORDER BY gid;", sep = "")
    if (create)
        str <- paste("CREATE TABLE ", table, " AS ", str)
    sqlQuery(db, str)
}