I work with R almost all the time, and now I am using it for doing spatial data mining.

I have a PostGIS database with (obviously) GIS data.

If I want to make statistical spatial analysis and plot maps is the better way to:

  • export the tables as shapefiles or;
  • work directly to the database?
up vote 32 down vote accepted

If you have PostGIS driver capability in the rgdal package then its just a question of creating a connection string and using that. Here I'm connecting to my local database gis using default credentials, so my DSN is rather simple. You might need to add a host, username, or password. See gdal docs for info.

> require(rgdal)
> dsn="PG:dbname='gis'"

What tables are in that database?

> ogrListLayers(dsn)
 [1] "ccsm_polygons"         "nongp"                 "WrldTZA"              
 [4] "nongpritalin"          "ritalinmerge"          "metforminmergev"      

Get one:

> polys = readOGR(dsn="PG:dbname='gis'","ccsm_polygons")
OGR data source with driver: PostgreSQL 
Source: "PG:dbname='gis'", layer: "ccsm_polygons"
with 32768 features and 4 fields
Feature type: wkbMultiPolygon with 2 dimensions

What have I got?

> summary(polys)
Object of class SpatialPolygonsDataFrame
Coordinates:
        min      max
x -179.2969 180.7031
y  -90.0000  90.0000
Is projected: NA 
proj4string : [NA]
Data attributes:
      area         perimeter       ccsm_polys      ccsm_pol_1   
 Min.   :1.000   Min.   :5.000   Min.   :    2   Min.   :    1  
 1st Qu.:1.000   1st Qu.:5.000   1st Qu.: 8194   1st Qu.: 8193  
 Median :1.000   Median :5.000   Median :16386   Median :16384  
 Mean   :1.016   Mean   :5.016   Mean   :16386   Mean   :16384  
 3rd Qu.:1.000   3rd Qu.:5.000   3rd Qu.:24577   3rd Qu.:24576  
 Max.   :2.000   Max.   :6.000   Max.   :32769   Max.   :32768  

Otherwise you can use R's database functionality and query the tables directly.

> require(RPostgreSQL)
Loading required package: RPostgreSQL
Loading required package: DBI
> m <- dbDriver("PostgreSQL")
> con <- dbConnect(m, dbname="gis")
> q="SELECT ST_AsText(the_geom) AS geom from ccsm_polygons LIMIT 10;"
> rs = dbSendQuery(con,q)
> df = fetch(rs,n=-1)

That returns the feature geometry in df$geom, which you'll need to convert to sp class objects (SpatialPolygons, SpatialPoints, SpatialLines) to do anything with. The readWKT function in rgeos can help with that.

The things to beware of are usually stuff like database columns that can't be mapped to R data types. You can include SQL in the query to do conversions, filtering, or limiting. This should get you started though.

  • Great answer, but how I enable the capability (the Postgis driver) in rgadl? I am in Ubuntu 13.04... – nanounanue Jul 2 '13 at 20:15
  • Do you have it? The ogrDrivers() function should tell you somewhere. If not then that's a whole other question (probably best googled for first and then asked on R-sig-geo) – Spacedman Jul 2 '13 at 22:47
  • In Ubuntu, the driver is installed by default. That is not the case in MacOS X. Thanks! – nanounanue Jul 2 '13 at 23:28
  • In your code above, Is it possible in the readOGR method use a sql instead of a full table? – nanounanue Jul 2 '13 at 23:28
  • Currently I think not. There was some chatter on r-sig-geo about 2.5y ago about this, but nothing seems to have been done. It looks simple to add a where clause and pass that to OGR via setAttributeFilter but that all has to be done in C and C++ code... – Spacedman Jul 3 '13 at 8:14

If you have data in Postgis, don't export it to shapefile. From my point of view, it's kind of a step back.

You can query your postgis database from R using SQL statements, importing them as dataframes and, since you are familiar with R, do all the geostatistics you need from there. I believe you can also export your geostatistic result back to postgis.

Using SQL with Postgis Functions you can also do all kind of spatial analysis, like overlay operations, distances,and so on.

For map plotting I would use QGIS, a OpenSource GIS software, that can read postgis directly (as far as I know that was the initial goal of the project), and the upcoming version 2.0 comes with lots of features to produce great looking maps.

  • Ok great advice, but since I want to automate everything in R (including plots) going to QGis, breaks the flow, isn´t it? – nanounanue Jul 1 '13 at 23:26
  • In that case, if you are comfortable with it, just use R to plot your maps. Even so, having the qgis projects layouts prepared based on the postgis (updated) data, they would update as well. I guess that in the end it will be a personal choice whether to use R or QGIS. – Alexandre Neto Jul 1 '13 at 23:48
  • Thank you for your quick response, but how can I make a plot using R, from a table in Postgis? – nanounanue Jul 1 '13 at 23:50
  • I'm not very experienced with R and I don't know how would you plot vector data using it (like I said I use QGIS for that), how do you plot shapfiles in R? For connecting to PostgresSQL from R I have used RPostgreSQL before. I think rgdal]. Good luck! – Alexandre Neto Jul 2 '13 at 9:11

You can use all the tools at the same time based for each step for your solution.

  • If you want to do geostaticstical analysis, use R. R's packages are more robust and allows you for a more analytical result. You can import data based on SQL querries.
  • If you want to aggregate your data based on a logical basis you can use PostGIS. You can answer complex queries like which many points are within my prescribed limits? But on grand scale.
  • For mapping, you can use either R or QGIS. QGIS is more straight forward, with R you might strugle for achieving the desired result.

We could provide you for a more specific answer if you'd give us more details from your problem

  • Could you provide an example of the last point, I mean, how can I do if I want to plot a map with R from a table in Postgis? – nanounanue Jul 2 '13 at 3:55
  • @nanounanue sure: library("rgdal") mydata = readOGR(dsn="PG:dbname=<mydb>",layer="schema.table") plot(mydata,axes=TRUE) title("My Plot"). – nickves Jul 2 '13 at 9:17
  • also take a look at this page: wiki.intamap.org/index.php/PostGIS – nickves Jul 2 '13 at 9:20

I would also go on a combination of rgdal and RPostgreSQL. So, same code as @Guillaume, except with a tryCatch that handle more lines, a pseudo-random table name and the use of an unlogged table for better performance. (Note to myself: we can't use TEMP table, because it's not visible from readOGR )

dbGetSp <- function(dbInfo,query) {
 if(!require('rgdal')|!require(RPostgreSQL))stop('missing rgdal or RPostgreSQL')
  d <- dbInfo
  tmpTbl <- sprintf('tmp_table_%s',round(runif(1)*1e5))
  dsn <- sprintf("PG:dbname='%s' host='%s' port='%s' user='%s' password='%s'",
    d$dbname,d$host,d$port,d$user,d$password
    )
  drv <- dbDriver("PostgreSQL")
  con <- dbConnect(drv, dbname=d$dbname, host=d$host, port=d$port,user=d$user, password=d$password)
  tryCatch({
    sql <- sprintf("CREATE UNLOGGED TABLE %s AS %s",tmpTbl,query)
    res <- dbSendQuery(con,sql)
    nr <- dbGetInfo(res)$rowsAffected
    if(nr<1){
      warning('There is no feature returned.');
      return()
    }
    sql <- sprintf("SELECT f_geometry_column from geometry_columns WHERE f_table_name='%s'",tmpTbl)
    geo <- dbGetQuery(con,sql)
    if(length(geo)>1){
      tname <- sprintf("%s(%s)",tmpTbl,geo$f_geometry_column[1])
    }else{
      tname <- tmpTbl;
    }
    out <- readOGR(dsn,tname)
    return(out)
  },finally={
    sql <- sprintf("DROP TABLE %s",tmpTbl)
    dbSendQuery(con,sql)
    dbClearResult(dbListResults(con)[[1]])
    dbDisconnect(con)
  })
}

Usage:

d=list(host='localhost', dbname='spatial_db', port='5432', user='myusername', password='mypassword')
spatialObj<-dbGetSp(dbInfo=d,"SELECT * FROM spatial_table")

But, this is still painfully slow:

For a small set of polygons (6 features, 22 fields) :

postgis part:

user  system elapsed
0.001   0.000   0.008

readOGR part:

user  system elapsed
0.313   0.021   1.436

You can also combine rgdal and RPostreSQL. This example function creates a temporary table with RPostgreSQL and sends it to readOGR for output of a spatial object. This is really inefficient and ugly, but it works quite well. Note that the query has to be a SELECT query and the user needs to have write access to the database.

RPostGIS <- function(coninfo,query) {
  dsn=paste("PG:dbname='",coninfo$dbname,"' host='",coninfo$host,"' port='",coninfo$port,"' user='",coninfo$user,"' password='",coninfo$password,"'", sep='')
  drv <- dbDriver("PostgreSQL")
  con <- dbConnect(drv, user=coninfo$user, password=coninfo$password, dbname=coninfo$dbname)
  res <- dbSendQuery(con,paste('CREATE TABLE tmp1209341251dva1 AS ',query,sep=''))
  geo <- dbGetQuery(con,"SELECT f_geometry_column from geometry_columns WHERE f_table_name='tmp1209341251dva1'")
  if(length(geo)>1){
    tname=paste("tmp1209341251dva1(",geo$f_geometry_column[1],")")
  }else{
    tname="tmp1209341251dva1";
  }
  out <- tryCatch(readOGR(dsn,tname), finally=dbSendQuery(con,'DROP TABLE tmp1209341251dva1'))
  dbDisconnect(con)
  return(out)
}

You can call it with something like:

> require('rgdal')
> require('RPostgreSQL')
> coninfo=list(host='localhost',dbname='spatial_db',port='5432',user='myusername',password='mypassword')
> spatial_obj<-RPostGIS(coninfo,"SELECT * FROM spatial_table")

There is now a RPostGIS package which can import PostGIS geoms into R with SQL queries.

The newly introduced sf-package (succesor of sp) provides the st_read() and st_read_db() functions. After this tutorial and from my experience it's faster than the already mentioned ways. As sf will probably replace sp one day it's also a good call to have look now ;)

require(sf)
dsn = "PG:dbname='dbname' host='host' port='port' user='user' password='pw'"
st_read(dsn, "schema.table")

you can also access the DB using RPostgreSQL:

require(sf)
require(RPostgreSQL)
drv <- dbDriver("PostgreSQL")
con <- dbConnect(drv, dbname = dbname, user = user, host = host, port = port, password = pw)

st_read_db(con, table = c("schema", "table"))
# or:
st_read_db(con, query = "SELECT * FROM schema.table")

dbDisconnect(con)
dbUnloadDriver(drv)

With st_write() you can upload data.

If you return a query with 'ST_AsText(geom) as geomwkt' and fetch the result into data, you can use:

library(rgeos);library(sp)
wkt_to_sp <- function(data) {
  #data is data.frame from postgis with geomwkt as only geom
  SpP <- SpatialPolygons(lapply(1:length(data$geomwkt), 
           function(x) Polygons(list(Polygon(readWKT(data$geomwkt[x]))),x)))
  data <- data[,!(names(data) == "geomwkt")]
  return(SpatialPolygonsDataFrame(SpP, data))
}

Still painfully slow.... 1 second for 100 geoms on a test.

Geotuple - https://github.com/rhansson/geotuple is a web app that connects R-Server and PostGIS (using RPostgreSQL)

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.