Pracujesz z danymi PostGIS w R?

27

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?
nanounanue
źródło

Odpowiedzi:

34

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ć.

Spacedman
źródło
Świetna odpowiedź, ale jak włączyć tę funkcję (sterownik Postgis) w rgadl? Jestem w Ubuntu 13.04 ...
nanounanue
Masz to? Funkcja ogrDrivers () powinna gdzieś ci powiedzieć. Jeśli nie, to jest to zupełnie inne pytanie (prawdopodobnie najlepiej google na najpierw, a następnie zadane na R-sig-geo)
Spacedman
W Ubuntu sterownik jest instalowany domyślnie. W MacOS X tak nie jest. Dzięki!
nanounanue
Czy w powyższym kodzie jest możliwe readOGRużycie metody sql zamiast pełnej tabeli?
nanounanue
Obecnie myślę, że nie. Około 2,5 lat temu na r-sig-geo było trochę gadania na ten temat, ale wydaje się, że nic nie zostało zrobione. Dodanie whereklauzuli i przekazanie jej do OGR wydaje się proste, setAttributeFilterale wszystko to musi zostać wykonane w kodzie C i C ++ ...
Spacedman
8

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 .

Alexandre Neto
źródło
Dobra świetna rada, ale ponieważ chcę zautomatyzować wszystko w R (w tym wykresy), przechodząc do QGis, przerywa przepływ, prawda?
nanounanue
W takim przypadku, jeśli czujesz się swobodnie, po prostu użyj R, aby wykreślić swoje mapy. Mimo to, mając przygotowane szablony projektów qgis na podstawie danych postgis (zaktualizowanych), również by się zaktualizowały. Sądzę, że ostatecznie będzie to osobisty wybór, czy użyć R czy QGIS.
Alexandre Neto,
Dziękuję za szybką odpowiedź, ale jak mogę stworzyć fabułę za pomocą R ze stołu w Postgis?
nanounanue
Nie mam dużego doświadczenia z R i nie wiem, w jaki sposób wykreślisz dane wektorowe przy użyciu tych danych (tak jak powiedziałem, że używam do tego QGIS), w jaki sposób wykreślasz pliki shap w R? Do łączenia się z PostgresSQL z RI wcześniej korzystałem z RPostgreSQL . Myślę, że rgdal ]. Powodzenia!
Alexandre Neto,
5

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.

andschar
źródło
1
To najprostsze rozwiązanie, jest winieta cran.r-project.org/web/packages/sf/vignettes/sf2.html wyjaśniająca, jak używać sf
Cedric
2

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

nickves
źródło
Czy mógłbyś podać przykład ostatniego punktu, to znaczy, jak mogę zrobić, jeśli chcę narysować mapę z R ze stołu w Postgis?
nanounanue
@nanounanue sure: biblioteka ("rgdal") mydata = readOGR (dsn = "PG: dbname = <mydb>", layer = "schema.table") działka (mydata, axes = TRUE) tytuł („My Plot”).
nickves
spójrz również na tę stronę: wiki.intamap.org/index.php/PostGIS
nickves
2

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
Fred Moser
źródło
2

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

Guillaume
źródło
1

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")
Guillaume
źródło
0

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.

ideamotor
źródło