Travailler avec des données PostGIS dans R?

Je travaille avec R presque tout le temps, et maintenant je l'utilise pour faire de l'exploration de données spatiales.

J'ai une base de données PostGIS avec (évidemment) des données SIG.

Si je veux faire une analyse spatiale statistique et tracer des cartes est la meilleure façon de:

• exporter les tables en tant que fichiers de formes ou;
• travailler directement dans la base de données?

Si vous avez une capacité de pilote PostGIS dans le package rgdal, il suffit de créer une chaîne de connexion et de l'utiliser. Ici, je me connecte à ma base de données locale en gisutilisant des informations d'identification par défaut, donc mon DSN est plutôt simple. Vous devrez peut-être ajouter un hôte, un nom d'utilisateur ou un mot de passe. Voir la documentation de gdal pour plus d'informations.

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

Quelles tables se trouvent dans cette base de données?

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

Prends en un:

> 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

Qu'est-ce que j'ai?

> 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  

Sinon, vous pouvez utiliser la fonctionnalité de base de données de R et interroger directement les tables.

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

Cela renvoie la géométrie de l' entité dans df$geomlaquelle vous devrez convertir en spobjets de classe (SpatialPolygons, SpatialPoints, SpatialLines) pour faire quoi que ce soit avec. La fonction readWKT dans rgeos peut vous aider.

Les éléments à prendre en compte sont généralement des éléments tels que des colonnes de base de données qui ne peuvent pas être mappés à des types de données R. Vous pouvez inclure SQL dans la requête pour effectuer des conversions, filtrer ou limiter. Cela devrait cependant vous aider à démarrer.

Si vous avez des données dans Postgis, ne les exportez pas vers le fichier de formes. De mon point de vue, c'est une sorte de recul.

Vous pouvez interroger votre base de données Postgis à partir de R à l'aide d'instructions SQL, les importer en tant que cadres de données et, comme vous connaissez R, faire toutes les géostatistiques dont vous avez besoin à partir de là. Je pense que vous pouvez également réexporter votre résultat géostatistique vers postgis.

En utilisant SQL avec les fonctions Postgis, vous pouvez également effectuer toutes sortes d'analyses spatiales, telles que des opérations de superposition, des distances, etc.

Pour le traçage de cartes, j'utiliserais QGIS , un logiciel SIG OpenSource, qui peut lire directement les postgis (pour autant que je sache que c'était l'objectif initial du projet), et la prochaine version 2.0 est livrée avec de nombreuses fonctionnalités pour produire de superbes cartes .

Le nouveau sf-package (successeur de sp) fournit les fonctions st_read()et st_read_db(). Après ce tutoriel et d'après mon expérience, c'est plus rapide que les moyens déjà mentionnés. Comme sf remplacera probablement sp un jour c'est aussi un bon appel à regarder maintenant;)

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

vous pouvez également accéder à la base de données en utilisant 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)

Avec st_write()vous pouvez télécharger des données.

Vous pouvez utiliser tous les outils en même temps en fonction de chaque étape de votre solution.

• Si vous voulez faire une analyse géostatique, utilisez les packages de R. R qui sont plus robustes et vous permettent un résultat plus analytique. Vous pouvez importer des données basées sur des requêtes SQL.
• Si vous souhaitez agréger vos données sur une base logique, vous pouvez utiliser PostGIS. Vous pouvez répondre à des questions complexes comme combien de points sont dans mes limites prescrites? Mais à grande échelle.
• Pour le mappage, vous pouvez utiliser R ou QGIS. QGIS est plus simple, avec R, vous pourriez avoir du mal à atteindre le résultat souhaité.

Nous pourrions vous fournir une réponse plus spécifique si vous nous donniez plus de détails sur votre problème

J'irais également sur une combinaison de rgdal et RPostgreSQL. Donc, même code que @Guillaume, sauf avec un tryCatch qui gère plus de lignes, un nom de table pseudo-aléatoire et l'utilisation d'une table non enregistrée pour de meilleures performances. (Remarque pour moi-même: nous ne pouvons pas utiliser la table TEMP, car elle n'est pas visible depuis 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")

Mais, c'est encore douloureusement lent:

Pour un petit ensemble de polygones (6 entités, 22 champs):

partie postgis:

user  system elapsed
0.001   0.000   0.008

partie readOGR:

user  system elapsed
0.313   0.021   1.436

Il existe maintenant un package RPostGIS qui peut importer des géomètres PostGIS dans R avec des requêtes SQL.

Vous pouvez également combiner rgdal et RPostreSQL. Cet exemple de fonction crée une table temporaire avec RPostgreSQL et l'envoie à readOGR pour la sortie d'un objet spatial. C'est vraiment inefficace et moche, mais cela fonctionne assez bien. Notez que la requête doit être une requête SELECT et que l'utilisateur doit avoir un accès en écriture à la base de données.

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

Vous pouvez l'appeler avec quelque chose comme:

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

Si vous renvoyez une requête avec 'ST_AsText (geom) as geomwkt' et récupérez le résultat dans les données, vous pouvez utiliser:

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

Toujours douloureusement lent .... 1 seconde pour 100 geoms sur un test.

Geotuple - https://github.com/rhansson/geotuple est une application web qui connecte R-Server et PostGIS (en utilisant RPostgreSQL)