Travailler avec des données PostGIS dans R?

27

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?
nanounanue
la source

Réponses:

34

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.

Spacedman
la source
Excellente réponse, mais comment activer la fonctionnalité (le pilote Postgis) dans rgadl? Je suis dans Ubuntu 13.04 ...
nanounanue
L'avez vous? La fonction ogrDrivers () devrait vous indiquer quelque part. Sinon, c'est une toute autre question (probablement mieux googlé pour la première fois, puis posé sur R-sig-geo)
Spacedman
Dans Ubuntu, le pilote est installé par défaut. Ce n'est pas le cas sous MacOS X. Merci!
nanounanue
Dans votre code ci-dessus, est-il possible dans la readOGRméthode d'utiliser un sql au lieu d'une table complète?
nanounanue
Actuellement, je ne pense pas. Il y a eu des discussions sur r-sig-geo il y a environ 2,5 ans à ce sujet, mais rien ne semble avoir été fait. Il semble simple d'ajouter une whereclause et de la transmettre à OGR via setAttributeFiltermais que tout doit être fait en code C et C ++ ...
Spacedman
8

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 .

Alexandre Neto
la source
Bon conseil, mais comme je veux automatiser tout dans R (y compris les tracés) en allant vers QGis, ça brise le flux, n'est-ce pas?
nanounanue
Dans ce cas, si cela vous convient, utilisez simplement R pour tracer vos cartes. Même ainsi, une fois les mises en page des projets qgis préparées sur la base des données postgis (mises à jour), elles seraient également mises à jour. Je suppose qu'en fin de compte, ce sera un choix personnel d'utiliser R ou QGIS.
Alexandre Neto
Merci pour votre réponse rapide, mais comment puis-je faire un tracé en utilisant R, à partir d'une table dans Postgis?
nanounanue
Je ne suis pas très expérimenté avec R et je ne sais pas comment traceriez-vous des données vectorielles en l'utilisant (comme je l'ai dit, j'utilise QGIS pour cela), comment tracer des fichiers de formes dans R? Pour vous connecter à PostgresSQL depuis RI, vous avez déjà utilisé RPostgreSQL . Je pense rgdal ]. Bonne chance!
Alexandre Neto
5

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.

andschar
la source
1
Ceci est la solution la plus simple, il y a une vignette cran.r-project.org/web/packages/sf/vignettes/sf2.html expliquant comment utiliser sf
Cedric
2

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

nickves
la source
Pourriez-vous fournir un exemple du dernier point, je veux dire, comment puis-je faire si je veux tracer une carte avec R à partir d'une table dans Postgis?
nanounanue le
@nanounanue sûr: bibliothèque ("rgdal") mydata = readOGR (dsn = "PG: dbname = <mydb>", layer = "schema.table") plot (mydata, axes = TRUE) title ("My Plot").
nickves
jetez également un œil à cette page: wiki.intamap.org/index.php/PostGIS
nickves
2

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
Fred Moser
la source
2

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

Guillaume
la source
1

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")
Guillaume
la source
0

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.

ideamotor
la source