Работаете с данными PostGIS в R?

Я работаю с R почти все время, и теперь я использую его для анализа пространственных данных.

У меня есть база данных PostGIS с (очевидно) данными ГИС.

Если я хочу сделать статистический пространственный анализ и построить карты, то это лучший способ:

• экспортировать таблицы как шейп-файлы или;
• работать напрямую с базой данных?

Если у вас есть возможность драйвера PostGIS в пакете rgdal, тогда просто вопрос создания строки подключения и ее использования. Здесь я подключаюсь к своей локальной базе данных, gisиспользуя учетные данные по умолчанию, поэтому мой DSN довольно прост. Вам может понадобиться добавить хост, имя пользователя или пароль. Смотрите gdal документы для информации.

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

Какие таблицы в этой базе данных?

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

Получить один:

> 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

Что у меня есть?

> 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  

В противном случае вы можете использовать функциональность базы данных R и запрашивать таблицы напрямую.

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

Это возвращает геометрию объекта df$geom, которую вам нужно преобразовать в spобъекты класса (SpatialPolygons, SpatialPoints, SpatialLines), чтобы иметь с чем угодно. Функция readWKT в rgeos может помочь с этим.

Остерегайтесь вещей, таких как столбцы базы данных, которые нельзя сопоставить с типами данных R. Вы можете включить SQL в запрос для выполнения преобразований, фильтрации или ограничения. Это должно начать вас, хотя.

Если у вас есть данные в Postgis, не экспортируйте их в шейп-файл. С моей точки зрения, это своего рода шаг назад.

Вы можете запросить вашу базу данных postgis из R, используя операторы SQL, импортируя их как кадры данных, и, поскольку вы знакомы с R, выполните всю необходимую вам геостатистику. Я верю, что вы также можете экспортировать свой геостатистический результат обратно в postgis.

Используя SQL с функциями Postgis, вы также можете выполнять все виды пространственного анализа, такие как операции наложения, расстояния и так далее.

Для построения карт я бы использовал QGIS , программное обеспечение OpenSource GIS, которое может считывать postgis напрямую (насколько я знаю, что было первоначальной целью проекта), а в следующей версии 2.0 будет множество функций для создания великолепно выглядящих карт .

Нововведенная SF-пакет (succesor из зр) обеспечивает st_read()и st_read_db()функции. После этого урока и из моего опыта это быстрее, чем уже упомянутые способы. Поскольку sf, вероятно, заменит sp однажды, это также хороший вызов, чтобы посмотреть сейчас;)

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

Вы также можете получить доступ к БД с помощью 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)

С помощью st_write()вы можете загружать данные.

Вы можете использовать все инструменты одновременно для каждого шага вашего решения.

• Если вы хотите провести геостатический анализ, используйте пакеты R. R, которые более надежны и позволяют получить более аналитический результат. Вы можете импортировать данные на основе запросов SQL.
• Если вы хотите агрегировать свои данные на логической основе, вы можете использовать PostGIS. Вы можете ответить на сложные вопросы, например, сколько точек находятся в моих установленных пределах? Но в большом масштабе.
• Для отображения вы можете использовать R или QGIS. QGIS более прост, с R вы можете бороться за достижение желаемого результата.

Мы могли бы предоставить вам более конкретный ответ, если вы дадите нам более подробную информацию о вашей проблеме

Я также пошел бы на комбинацию rgdal и RPostgreSQL. Таким образом, тот же код, что и у @Guillaume, за исключением tryCatch, который обрабатывает больше строк, псевдослучайное имя таблицы и использование незафиксированной таблицы для повышения производительности. (Примечание для себя: мы не можем использовать таблицу TEMP, потому что она не видна из 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)
  })
}

Использование:

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

Но это все еще мучительно медленно:

Для небольшого набора полигонов (6 объектов, 22 поля):

Postgis часть:

user  system elapsed
0.001   0.000   0.008

часть readOGR:

user  system elapsed
0.313   0.021   1.436

Теперь есть пакет RPostGIS , который может импортировать геомы PostGIS в R с помощью SQL-запросов.

Вы также можете комбинировать rgdal и RPostreSQL. В этом примере функция создает временную таблицу с RPostgreSQL и отправляет ее в readOGR для вывода пространственного объекта. Это действительно неэффективно и некрасиво, но работает довольно хорошо. Обратите внимание, что запрос должен быть запросом SELECT, а пользователь должен иметь доступ для записи в базу данных.

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

Вы можете назвать это чем-то вроде:

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

Если вы возвращаете запрос с ST_AsText (geom) как geomwkt и извлекаете результат в данные, вы можете использовать:

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

Все еще мучительно медленно ... 1 секунда для 100 гемов на тесте.

Geotuple - https://github.com/rhansson/geotuple - это веб-приложение, которое соединяет R-Server и PostGIS (используя RPostgreSQL)