Trabalhando com dados PostGIS no R?

27

Trabalho com R quase o tempo todo e agora estou usando-o para fazer mineração de dados espaciais.

Eu tenho um banco de dados PostGIS com (obviamente) dados GIS.

Se eu quiser fazer análises espaciais estatísticas e mapas de plotagem, é a melhor maneira de:

  • exportar as tabelas como shapefiles ou;
  • trabalhar diretamente para o banco de dados?
nanounanue
fonte

Respostas:

34

Se você possui a capacidade do driver PostGIS no pacote rgdal, é apenas uma questão de criar uma string de conexão e usá-la. Aqui estou me conectando ao meu banco de dados local gisusando credenciais padrão, portanto meu DSN é bastante simples. Pode ser necessário adicionar um host, nome de usuário ou senha. Veja os documentos da gdal para obter informações.

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

Quais tabelas estão nesse banco de dados?

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

Pegue um:

> 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

O que eu tenho?

> 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  

Caso contrário, você pode usar a funcionalidade do banco de dados de R e consultar as tabelas diretamente.

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

Isso retorna a geometria do recurso df$geom, na qual você precisará converter em spobjetos de classe (SpatialPolygons, SpatialPoints, SpatialLines) para fazer qualquer coisa. A função readWKT nos rgeos pode ajudar com isso.

O que você deve observar são coisas como colunas de banco de dados que não podem ser mapeadas para tipos de dados R. Você pode incluir SQL na consulta para fazer conversões, filtrar ou limitar. Isso deve ajudá-lo a começar.

Spacedman
fonte
Ótima resposta, mas como habilito o recurso (o driver Postgis) rgadl? Estou no Ubuntu 13.04 ...
nanounanue 02/07/2013
Você tem? A função ogrDrivers () deve indicar algum lugar. Se não, então isso é uma outra questão (provavelmente melhor pesquisei para o primeiro e, em seguida, perguntou sobre R-sig-geo)
Spacedman
No Ubuntu, o driver é instalado por padrão. Esse não é o caso no MacOS X. Obrigado!
Nanounanue 02/07/2013
No seu código acima, no readOGRmétodo é possível usar um sql em vez de uma tabela completa?
Nanounanue 02/07/2013
Atualmente acho que não. Houve alguma conversa sobre r-sig-geo cerca de 2,5 anos atrás sobre isso, mas nada parece ter sido feito. Parece simples para adicionar uma wherecláusula e passar isso para OGR via setAttributeFilter, mas que todos tem que ser feito em código C e C ++ ...
Spacedman
8

Se você possui dados no Postgis, não os exporte para o shapefile. Do meu ponto de vista, é meio que um passo atrás.

Você pode consultar seu banco de dados postgis a partir do R usando instruções SQL, importando-os como quadros de dados e, como você está familiarizado com o R, faça todas as geoestatísticas necessárias a partir daí. Acredito que você também pode exportar seu resultado geoestatístico de volta para postgis.

Usando SQL com funções Postgis, você também pode fazer todo o tipo de análise espacial, como operações de sobreposição, distâncias e assim por diante.

Para a plotagem de mapas, eu usaria o QGIS , um software GIS de código aberto, que pode ler postgis diretamente (até onde eu sei que era o objetivo inicial do projeto), e a próxima versão 2.0 vem com muitos recursos para produzir mapas de ótima aparência .

Alexandre Neto
fonte
Ok, ótimo conselho, mas como eu quero automatizar tudo em R (incluindo parcelas) que vão para o QGis, interrompe o fluxo, não é?
Nanounanue 01/07/2013
Nesse caso, se você estiver confortável com isso, use R para plotar seus mapas. Mesmo assim, tendo os layouts dos projetos qgis preparados com base nos dados postgis (atualizados), eles também seriam atualizados. Eu acho que, no final, será uma escolha pessoal usar R ou QGIS.
Alexandre Neto
Obrigado pela sua resposta rápida, mas como posso fazer um gráfico usando R, de uma tabela no Postgis?
Nanounanue 01/07/2013
Não tenho muita experiência com R e não sei como você plotaria dados vetoriais usando-os (como eu disse que utilizo QGIS para isso), como plotamos shapfiles em R? Para conectar-se ao PostgresSQL da RI, usei o RPostgreSQL antes. Eu acho que rgdal ]. Boa sorte!
Alexandre Neto
5

O recém-introduzido pacote sf (sucessor de sp) fornece as funções st_read()e st_read_db(). Após este tutorial e pela minha experiência, é mais rápido do que as formas já mencionadas. Como o sf provavelmente substituirá sp um dia, também é uma boa ideia dar uma olhada agora;)

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

você também pode acessar o banco de dados usando o 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)

Com st_write()você pode fazer upload de dados.

andschar
fonte
1
Esta é a solução mais simples, existe uma vinheta cran.r-project.org/web/packages/sf/vignettes/sf2.html explicando como usar sf
Cedric
2

Você pode usar todas as ferramentas ao mesmo tempo com base em cada etapa da sua solução.

  • Se você deseja fazer análises geoestáticas, use os pacotes da R.R são mais robustos e permitem um resultado mais analítico. Você pode importar dados com base em consultas SQL.
  • Se você deseja agregar seus dados com base em uma lógica, você pode usar o PostGIS. Você pode responder a perguntas complexas, como quantos pontos estão dentro dos meus limites prescritos? Mas em grande escala.
  • Para o mapeamento, você pode usar R ou QGIS. O QGIS é mais direto, com R você pode lutar para alcançar o resultado desejado.

Poderíamos fornecer uma resposta mais específica se você nos desse mais detalhes sobre o seu problema

nickves
fonte
Você poderia fornecer um exemplo do último ponto, quero dizer, como posso fazer se eu quiser traçar um mapa com R a partir de uma tabela no Postgis?
Nanounanue 02/07/2013
@nanounanue sure: library ("rgdal") mydata = readOGR (dsn = "PG: dbname = <mydb>", layer = "schema.table") plot (mydata, axes = TRUE) title ("Meu gráfico").
nickves
também dê uma olhada nesta página: wiki.intamap.org/index.php/PostGIS
nickves
2

Eu também usaria uma combinação de rgdal e RPostgreSQL. Portanto, o mesmo código que o @Guillaume, exceto com um tryCatch que lida com mais linhas, um nome de tabela pseudoaleatório e o uso de uma tabela não registrada para obter melhor desempenho. (Observação: não podemos usar a tabela TEMP, porque não é visível no 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)
  })
}

Uso:

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

Mas isso ainda é dolorosamente lento:

Para um pequeno conjunto de polígonos (6 recursos, 22 campos):

parte postgis:

user  system elapsed
0.001   0.000   0.008

parte readOGR:

user  system elapsed
0.313   0.021   1.436
Fred Moser
fonte
2

Agora existe um pacote RPostGIS que pode importar geoms PostGIS para o R com consultas SQL.

Guillaume
fonte
1

Você também pode combinar rgdal e RPostreSQL. Esta função de exemplo cria uma tabela temporária com o RPostgreSQL e a envia ao readOGR para a saída de um objeto espacial. Isso é realmente ineficiente e feio, mas funciona muito bem. Observe que a consulta deve ser uma consulta SELECT e o usuário precisa ter acesso de gravação ao banco de dados.

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

Você pode chamá-lo com algo como:

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

Se você retornar uma consulta com 'ST_AsText (geom) como geomwkt' e buscar o resultado nos dados, poderá usar:

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

Ainda é dolorosamente lento .... 1 segundo para 100 geoms em um teste.

ideamotor
fonte