Arbeiten mit PostGIS-Daten in R?

27

Ich arbeite fast die ganze Zeit mit R und nutze es jetzt für das räumliche Data Mining.

Ich habe eine PostGIS-Datenbank mit (offensichtlich) GIS-Daten.

Wenn ich statistische räumliche Analysen und Karten machen möchte, ist der bessere Weg:

  • Exportieren Sie die Tabellen als Shapefiles oder;
  • direkt an der Datenbank arbeiten?
nanounanue
quelle

Antworten:

34

Wenn Sie PostGIS-Treiber im rgdal-Paket haben, müssen Sie nur eine Verbindungszeichenfolge erstellen und diese verwenden. Hier verbinde ich mich mit meinen gisStandardanmeldeinformationen mit meiner lokalen Datenbank , daher ist mein DSN ziemlich einfach. Möglicherweise müssen Sie einen Host, einen Benutzernamen oder ein Kennwort hinzufügen. Weitere Informationen finden Sie in den GDAL-Dokumenten.

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

Welche Tabellen befinden sich in dieser Datenbank?

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

Holen Sie sich eines:

> 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

Was habe ich

> 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  

Andernfalls können Sie die Datenbankfunktionalität von R verwenden und die Tabellen direkt abfragen.

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

Das gibt die Feature-Geometrie zurück df$geom, in die Sie konvertieren müssen, um mit spKlassenobjekten (SpatialPolygons, SpatialPoints, SpatialLines) etwas zu tun. Die readWKT-Funktion in rgeos kann dabei helfen.

Die Dinge, auf die Sie achten müssen, sind normalerweise Dinge wie Datenbankspalten, die nicht auf R-Datentypen abgebildet werden können. Sie können SQL in die Abfrage einbeziehen, um Konvertierungen, Filterungen oder Einschränkungen vorzunehmen. Damit sollten Sie jedoch beginnen.

Raumfahrer
quelle
Tolle Antwort, aber wie aktiviere ich die Fähigkeit (den Postgis-Treiber) in rgadl? Ich bin in Ubuntu 13.04 ...
Nanounanue
Hast Du es? Die Funktion ogrDrivers () sollte Ihnen irgendwo Bescheid geben. Wenn nicht, dann ist das eine ganz andere Frage (wahrscheinlich am besten zuerst gegoogelt und dann auf R-Sig-Geo gestellt)
Spacedman
In Ubuntu ist der Treiber standardmäßig installiert. Dies ist unter MacOS X nicht der Fall. Danke!
Nanounanue
Ist es in Ihrem obigen Code möglich, in der readOGRMethode eine SQL anstelle einer vollständigen Tabelle zu verwenden?
Nanounanue
Derzeit denke ich nicht. Es gab vor ungefähr 2,5 Jahren ein Geschwätz über r-sig-geo, aber nichts scheint getan worden zu sein. Es sieht einfach aus, eine whereKlausel hinzuzufügen und diese über an OGR zu übergeben, setAttributeFilteraber alles muss in C- und C ++ - Code erfolgen ...
Spacedman,
8

Wenn Sie Daten in Postgis haben, exportieren Sie diese nicht in Shapefile. Aus meiner Sicht ist es eine Art Rückschritt.

Sie können Ihre Postgis-Datenbank mit SQL-Anweisungen aus R abfragen, sie als Datenrahmen importieren und, da Sie mit R vertraut sind, alle Geostatistiken ausführen, die Sie von dort aus benötigen. Ich glaube, Sie können Ihr geostatistisches Ergebnis auch wieder in Postgis exportieren.

Mithilfe von SQL mit Postgis-Funktionen können Sie auch alle Arten von räumlichen Analysen durchführen, z. B. Überlagerungsoperationen, Entfernungen usw.

Zum Zeichnen von Karten würde ich QGIS verwenden , eine OpenSource-GIS-Software, die Postgis direkt lesen kann (soweit ich weiß, dass dies das ursprüngliche Ziel des Projekts war), und die kommende Version 2.0 wird mit vielen Funktionen geliefert, um großartig aussehende Karten zu erstellen .

Alexandre Neto
quelle
Ok, toller Rat, aber da ich alles in R automatisieren möchte (einschließlich Plots), gehe ich zu QGis, unterbricht den Fluss, nicht wahr?
Nanounanue
In diesem Fall können Sie mit R Ihre Karten zeichnen, wenn Sie damit vertraut sind. Wenn die Layouts der qgis-Projekte auf der Grundlage der (aktualisierten) Postgis-Daten erstellt werden, werden sie auch aktualisiert. Ich denke, dass es am Ende eine persönliche Entscheidung sein wird, ob Sie R oder QGIS verwenden.
Alexandre Neto
Vielen Dank für Ihre schnelle Antwort, aber wie kann ich mit R aus einer Tabelle in Postgis einen Plot erstellen?
Nanounanue
Ich bin nicht sehr erfahren mit R und ich weiß nicht, wie Sie Vektordaten damit zeichnen (wie ich bereits sagte, ich verwende QGIS dafür), wie Sie Shapfiles in R zeichnen. Für die Verbindung zu PostgresSQL von RI aus haben Sie zuvor RPostgreSQL verwendet . Ich denke rgdal ]. Viel Glück!
Alexandre Neto
5

Das neu eingeführte sf-Paket (Nachfolger von sp) bietet die Funktionen st_read()und st_read_db(). Nach diesem Tutorial und meiner Erfahrung nach ist es schneller als die bereits erwähnten Möglichkeiten. Da sf wahrscheinlich eines Tages sp ersetzen wird, ist es auch ein guter Anruf, jetzt nachzuschauen;)

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

Sie können auch mit RPostgreSQL auf die Datenbank zugreifen:

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)

Mit können st_write()Sie Daten hochladen.

andschar
quelle
1
Dies ist die einfachste Lösung. Es gibt eine Vignette unter cran.r-project.org/web/packages/sf/vignettes/sf2.html, in der erklärt wird, wie sf
Cedric
2

Sie können alle Tools für jeden Schritt Ihrer Lösung gleichzeitig verwenden.

  • Wenn Sie eine geostatische Analyse durchführen möchten, verwenden Sie die Pakete von R. R, die robuster sind und Ihnen ein besseres Analyseergebnis ermöglichen. Sie können Daten basierend auf SQL-Abfragen importieren.
  • Wenn Sie Ihre Daten logisch aggregieren möchten, können Sie PostGIS verwenden. Sie können komplexe Fragen beantworten, bei denen sich viele Punkte innerhalb meiner vorgegebenen Grenzen befinden? Aber im großen Stil.
  • Für das Mapping können Sie entweder R oder QGIS verwenden. QGIS ist einfacher, mit R könnten Sie Schwierigkeiten haben, das gewünschte Ergebnis zu erzielen.

Wir können Ihnen eine genauere Antwort geben, wenn Sie uns weitere Einzelheiten zu Ihrem Problem mitteilen

nickves
quelle
Könnten Sie ein Beispiel für den letzten Punkt geben, ich meine, wie kann ich vorgehen, wenn ich eine Karte mit R aus einer Tabelle in Postgis zeichnen möchte?
Nanounanue
@nanounanue sure: library ("rgdal") mydata = readOGR (dsn = "PG: dbname = <mydb>", layer = "schema.table") plot (mydata, axes = TRUE) title ("My Plot").
nickves
Schauen Sie sich auch diese Seite an: wiki.intamap.org/index.php/PostGIS
nickves
2

Ich würde auch eine Kombination aus rgdal und RPostgreSQL verwenden. Also derselbe Code wie bei @Guillaume, mit Ausnahme eines tryCatch, der mehr Zeilen verarbeitet, eines pseudozufälligen Tabellennamens und der Verwendung einer nicht protokollierten Tabelle für eine bessere Leistung. (Anmerkung für mich: Wir können die TEMP-Tabelle nicht verwenden, da sie in readOGR nicht sichtbar ist.)

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

Verwendung:

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

Aber das ist immer noch schmerzhaft langsam:

Für eine kleine Gruppe von Polygonen (6 Features, 22 Felder):

Postgis-Teil:

user  system elapsed
0.001   0.000   0.008

readOGR-Teil:

user  system elapsed
0.313   0.021   1.436
Fred Moser
quelle
2

Es gibt jetzt ein RPostGIS- Paket, mit dem PostGIS-Geoms mit SQL-Abfragen in R importiert werden können.

Guillaume
quelle
1

Sie können auch rgdal und RPostreSQL kombinieren. Diese Beispielfunktion erstellt eine temporäre Tabelle mit RPostgreSQL und sendet sie zur Ausgabe eines Geo-Objekts an readOGR. Das ist wirklich ineffizient und hässlich, aber es funktioniert ganz gut. Beachten Sie, dass die Abfrage eine SELECT-Abfrage sein muss und der Benutzer Schreibzugriff auf die Datenbank haben muss.

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

Sie können es mit so etwas wie aufrufen:

> 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
quelle
0

Wenn Sie eine Abfrage mit 'ST_AsText (geom) as geomwkt' zurückgeben und das Ergebnis in Daten abrufen, können Sie Folgendes verwenden:

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

Immer noch schmerzhaft langsam .... 1 Sekunde für 100 Geoms bei einem Test.

ideamotor
quelle