**Open Geospatial Technologies » R**, and kindly contributed to R-bloggers)

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])

**leave a comment**for the author, please follow the link and comment on their blog:

**Open Geospatial Technologies » R**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...