Pracujesz z danymi PostGIS w R?

Pracuję z R prawie cały czas, a teraz używam go do eksploracji danych przestrzennych.

Mam bazę danych PostGIS z (oczywiście) danymi GIS.

Jeśli chcę przeprowadzić statystyczną analizę przestrzenną i mapy działek, to lepszy sposób:

• eksportuj tabele jako pliki kształtów lub;
• pracować bezpośrednio w bazie danych?

Jeśli masz możliwość sterownika PostGIS w pakiecie rgdal, to jest tylko kwestia utworzenia ciągu połączenia i korzystania z niego. Tutaj łączę się z moją lokalną bazą danych gisprzy użyciu domyślnych poświadczeń, więc moje DSN jest raczej proste. Konieczne może być dodanie hosta, nazwy użytkownika lub hasła. Zobacz dokumentację gdal, aby uzyskać informacje.

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

Jakie tabele są w tej bazie danych?

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

Uzyskać jeden:

> 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

Co mam

> 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  

W przeciwnym razie możesz użyć funkcji bazy danych R i bezpośrednio przeszukać tabele.

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

Zwraca geometrię elementu, w df$geomktórej będziesz musiał przekonwertować na spobiekty klasy (SpatialPolygons, SpatialPoints, SpatialLines), aby cokolwiek zrobić. Pomaga w tym funkcja readWKT w rgeos.

Rzeczy, na które należy uważać, to zwykle takie rzeczy, jak kolumny bazy danych, których nie można zmapować na typy danych R. Do zapytania możesz dołączyć SQL, aby wykonywać konwersje, filtrowanie lub ograniczanie. To powinno wystartować.

Jeśli masz dane w Postgis, nie eksportuj ich do pliku kształtu. Z mojego punktu widzenia jest to pewien krok wstecz.

Możesz wysłać zapytanie do bazy danych Postgis z poziomu R, używając instrukcji SQL, importując je jako ramki danych i, ponieważ znasz R, wykonaj wszystkie potrzebne dane geostatystyczne. Wierzę, że możesz również wyeksportować swój wynik geostatystyczny z powrotem do postgis.

Korzystając z SQL z funkcjami Postgis, możesz także wykonywać wszelkiego rodzaju analizy przestrzenne, takie jak operacje nakładania, odległości itp.

Do tworzenia map użyłbym QGIS , oprogramowania OpenSource GIS, które może bezpośrednio czytać postgis (o ile wiem, że był to początkowy cel projektu), a nadchodząca wersja 2.0 zawiera wiele funkcji do tworzenia świetnie wyglądających map .

Nowo wprowadzony pakiet sf (następca sp) zapewnia funkcje st_read()i st_read_db(). Po tym samouczku i z mojego doświadczenia jest szybszy niż wspomniane już sposoby. Ponieważ pewnego dnia sf prawdopodobnie zastąpi sp, warto teraz zajrzeć;)

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

możesz również uzyskać dostęp do bazy danych za pomocą 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)

Dzięki st_write()niemu możesz przesyłać dane.

Możesz korzystać ze wszystkich narzędzi jednocześnie dla każdego kroku swojego rozwiązania.

• Jeśli chcesz przeprowadzić analizę geostatyczną, użyj pakietów R. R. jest bardziej solidny i pozwala uzyskać bardziej analityczny wynik. Możesz importować dane na podstawie zapytań SQL.
• Jeśli chcesz agregować swoje dane w oparciu o logikę, możesz skorzystać z PostGIS. Czy potrafisz odpowiedzieć na złożone pytania, np. Które punkty są w moich przepisanych granicach? Ale na wielką skalę.
• Do mapowania możesz użyć R lub QGIS. QGIS jest bardziej prosty, z R możesz walczyć o osiągnięcie pożądanego rezultatu.

Możemy udzielić bardziej szczegółowej odpowiedzi, jeśli podasz nam więcej szczegółów na temat swojego problemu

Wybrałbym także kombinację rgdal i RPostgreSQL. Zatem taki sam kod jak @Guillaume, z wyjątkiem tryCatch, który obsługuje więcej wierszy, pseudolosowej nazwy tabeli i użycia niezalogowanej tabeli dla lepszej wydajności. (Uwaga dla siebie: nie możemy użyć tabeli TEMP, ponieważ nie jest widoczna z 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)
  })
}

Stosowanie:

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

Ale wciąż jest to boleśnie powolne:

W przypadku małego zestawu wielokątów (6 obiektów, 22 pola):

część postgis:

user  system elapsed
0.001   0.000   0.008

część readOGR:

user  system elapsed
0.313   0.021   1.436

Dostępny jest teraz pakiet RPostGIS , który może importować geom PostGIS do R za pomocą zapytań SQL.

Możesz także łączyć rgdal i RPostreSQL. Ta przykładowa funkcja tworzy tymczasową tabelę z RPostgreSQL i wysyła ją do readOGR w celu uzyskania wyniku obiektu przestrzennego. Jest to naprawdę nieefektywne i brzydkie, ale działa całkiem dobrze. Pamiętaj, że zapytanie musi być zapytaniem SELECT, a użytkownik musi mieć dostęp do zapisu do bazy danych.

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

Możesz to nazwać:

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

Jeśli zwrócisz zapytanie z „ST_AsText (geom) jako geomwkt” i pobierzesz wynik do danych, możesz użyć:

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

Wciąż boleśnie powolne .... 1 sekunda na 100 geomów na teście.

Geotuple - https://github.com/rhansson/geotuple to aplikacja internetowa, która łączy R-Server i PostGIS (za pomocą RPostgreSQL)