¿Trabajando con datos de PostGIS en R?

26

Trabajo con R casi todo el tiempo, y ahora lo estoy usando para hacer minería de datos espaciales.

Tengo una base de datos PostGIS con (obviamente) datos GIS.

Si quiero hacer un análisis estadístico espacial y trazar mapas es la mejor manera de:

  • exportar las tablas como shapefiles o;
  • ¿trabaja directamente en la base de datos?
pregunta nanounanue 01.07.2013 - 23:09

9 respuestas

32

Si tiene la capacidad del controlador PostGIS en el paquete rgdal, entonces solo es cuestión de crear una cadena de conexión y usarla. Aquí me estoy conectando a mi base de datos local gis usando credenciales predeterminadas, por lo que mi DSN es bastante simple. Es posible que deba agregar un host, nombre de usuario o contraseña. Ver gdal docs para información.

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

¿Qué tablas hay en esa base de datos?

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

Consigue uno:

> 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

¿Qué tengo?

> 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  

De lo contrario, puede usar la funcionalidad de la base de datos de R y consultar las tablas directamente.

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

Eso devuelve la geometría de la característica en df$geom , que deberá convertir a objetos de clase sp (SpatialPolygons, SpatialPoints, SpatialLines) para hacer cualquier cosa con. La función readWKT en rgeos puede ayudar con eso.

Las cosas a tener en cuenta son cosas como columnas de base de datos que no pueden asignarse a los tipos de datos R. Puede incluir SQL en la consulta para hacer conversiones, filtrado o limitación. Esto debería ayudarte a empezar, sin embargo.

    
respondido por el Spacedman 02.07.2013 - 10:01
8

Si tiene datos en Postgis, no los exporte a shapefile. Desde mi punto de vista, es un paso atrás.

Puede consultar su base de datos de Postgis desde R usando sentencias de SQL, importándolas como marcos de datos y, ya que está familiarizado con R, haga todas las geoestadísticas que necesite desde allí. Creo que también puedes exportar tu resultado geoestadístico a postgis.

Al utilizar SQL con las funciones de Postgis, también puede hacer todo tipo de análisis espacial, como operaciones de superposición, distancias, etc.

Para el trazado de mapas, usaría QGIS , un software OpenSource GIS, que puede leer postgis directamente (hasta donde sé que era el objetivo inicial del proyecto), y la próxima versión 2.0 viene con muchas características para producir mapas de gran apariencia .

    
respondido por el Alexandre Neto 02.07.2013 - 00:30
2

Puede usar todas las herramientas al mismo tiempo para cada paso de su solución.

  • Si desea realizar un análisis geoestático, los paquetes de R. R son más robustos y le permiten obtener un resultado más analítico. Puede importar datos en función de las dudas de SQL.
  • Si desea agregar sus datos basándose en una base lógica, puede usar PostGIS. ¿Puede responder preguntas complejas como cuáles son muchos puntos dentro de mis límites prescritos? Pero a gran escala.
  • Para la asignación, puede usar R o QGIS. QGIS es más sencillo, con R puede que luches por lograr el resultado deseado.

Podríamos proporcionarle una respuesta más específica si nos proporcionara más detalles de su problema

    
respondido por el nickves 02.07.2013 - 01:28
2

También iría en una combinación de rgdal y RPostgreSQL. Por lo tanto, el mismo código que @Guillaume, excepto con un tryCatch que maneja más líneas, un nombre de tabla pseudoaleatoria y el uso de una tabla sin registro para un mejor rendimiento. (Nota para mí: no podemos usar la tabla TEMP, porque no está visible desde 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")

Pero, esto sigue siendo dolorosamente lento:

Para un pequeño conjunto de polígonos (6 entidades, 22 campos):

parte postgis:

user  system elapsed
0.001   0.000   0.008

parte de readOGR:

user  system elapsed
0.313   0.021   1.436
    
respondido por el Fred Moser 03.06.2015 - 16:15
2

El recién introducido sf-package (succesor de sp) proporciona el% co_de Funciones% y st_read() . Después de este tutorial y desde mi experiencia, es más rápido que las formas ya mencionadas. Como sf probablemente reemplazará a sp algún día, también es una buena idea echarle un vistazo ahora;)

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

también puede acceder a la base de datos utilizando 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)

Con st_read_db() puedes subir datos.

    
respondido por el andrasz 20.04.2017 - 16:58
1

También puedes combinar rgdal y RPostreSQL. Esta función de ejemplo crea una tabla temporal con RPostgreSQL y la envía a readOGR para la salida de un objeto espacial. Esto es realmente ineficiente y feo, pero funciona bastante bien. Tenga en cuenta que la consulta debe ser una consulta SELECT y que el usuario debe tener acceso de escritura a la base de datos.

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

Puedes llamarlo con 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")
    
respondido por el Guillaume 27.04.2015 - 21:59
1

Ahora hay un paquete RPostGIS que puede importar geografías de PostGIS a R con SQL consultas

    
respondido por el Guillaume 08.03.2017 - 17:29
0

Si devuelve una consulta con 'ST_AsText (geom) como geomwkt' y obtiene el resultado en datos, puede 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))
}

Sigue siendo dolorosamente lento ... 1 segundo para 100 geoms en una prueba.

    
respondido por el ideamotor 21.04.2016 - 21:41
0

Geotuple - enlace es una aplicación web que conecta R-Server y PostGIS (mediante RPostgreSQL)

    
respondido por el user1493801 20.04.2017 - 17:30

Lea otras preguntas en las etiquetas