Load PostGIS geometries in R without rgdal

Uncategorized,

As I said in my last post, rgdal lacks some of the features of GDAL, including the ability to subset columns and rows the source layer, and I demonstrated a workaround. The workaround relied upon the RPostgreSQL package, and this raises a question: Is it possible to transfer geographic data from PostGIS to R just using RPostgreSQL?

It is, and in fact, this may be necessary if you are working on Windows.

The rgdal package comes with its own statically linked GDAL, which unfortunately has a very limited set of format drivers—the ubiquitous ESRI Shapefile, but not the PostgreSQL/PostGIS driver. If you try to use dbGetQuery to return a recordset with a geometry column, R will choke. But you can use the ST_AsText function in SQL to convert the PostGIS geometry to Well-Known Text, read it into R, and then use the readWKT function in the rgeos package to convert it into a SpatialPolygons, SpatialLines, or SpatialPoints object.

The readWKT function is not vectorized, so I build the set of SpatialPolygons and rbind them inside a for loop. I’m sure there’s a better way to do this using the *apply family of functions, but I haven’t been able to figure that one out.

library(RPostgreSQL)
library(rgeos)
library(sp)
 
# Load data from the PostGIS server
conn = dbConnect(
  dbDriver("PostgreSQL"), dbname=dbname, host=host, port=5432, 
  user=user, password=password
  )
 
strSQL = "
  SELECT gid, ST_AsText(geom) AS wkt_geometry, attr1, attr2[, ...]
  FROM geo_layer"
dfTemp = dbGetQuery(conn, strSQL)
row.names(dfTemp) = dfTemp$gid
 
# Create spatial polygons
# To set the PROJ4 string, enter the EPSG SRID and uncomment the 
# following two lines:
# EPSG = make_EPSG()
# p4s = EPSG[which(EPSG$code == SRID), "prj4"]
for (i in seq(nrow(dfTemp))) {
  if (i == 1) {
    spTemp = readWKT(dfTemp$wkt_geometry[i], dfTemp$gid[i])
    # If the PROJ4 string has been set, use the following instead
    # spTemp = readWKT(dfTemp$wkt_geometry[i], dfTemp$gid[i], p4s)
  }
  else {
    spTemp = rbind(
      spTemp, readWKT(dfTemp$wkt_geometry[i], dfTemp$gid[i])
      # If the PROJ4 string has been set, use the following instead
      # spTemp, readWKT(dfTemp$wkt_geometry[i], dfTemp$gid[i], p4s)
    )
  }
}
 
# Create SpatialPolygonsDataFrame, drop WKT field from attributes
spdfFinal = SpatialPolygonsDataFrame(spTemp, dfTemp[-2])

Created by Pretty R at inside-R.org