我几乎一直都在使用R,现在我正在使用它进行空间数据挖掘。

我有一个PostGIS数据库,其中包含(显然)GIS数据。

如果我想进行统计空间分析和绘图是更好的方法:


将表导出为shapefile或;
直接工作到数据库?


#1 楼

如果您在rgdal软件包中具有PostGIS驱动程序功能,那么仅是创建连接字符串并使用它的问题。在这里,我使用默认凭据连接到本地数据库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)以执行任何操作。 rgeos中的readWKT函数可以帮助解决此问题。

要提防的事情通常是诸如无法映射到R数据类型的数据库列之类的东西。您可以在查询中包含SQL以进行转换,过滤或限制。不过,这应该可以帮助您入门。

评论


很好的答案,但是如何在rgadl中启用该功能(Postgis驱动程序)?我在Ubuntu 13.04中...

– Nanounanue
13年7月2日在20:15

你有吗? ogrDrivers()函数应该告诉您某个地方。如果不是这样的话,那是另一个问题(可能最好先用Google搜索,然后在R-sig-geo上询问)

– Spacedman
13年7月2日在22:47

在Ubuntu中,默认情况下已安装驱动程序。在MacOS X中不是这种情况。谢谢!

– Nanounanue
13年7月2日在23:28

在上面的代码中,是否可以在readOGR方法中使用sql而不是完整表?

– Nanounanue
13年7月2日在23:28

目前我认为不是。关于此事,大约在2.5年前,r-sig-geo上有一些chat不休,但似乎什么也没做。添加where子句并通过setAttributeFilter将其传递给OGR看起来很简单,但是所有这些都必须使用C和C ++代码来完成...

– Spacedman
13年7月3日在8:14

#2 楼

如果Postgis中有数据,请不要将其导出到shapefile。从我的角度来看,这有点后退。

您可以使用SQL语句从R查询Postgis数据库,并将其作为数据帧导入,并且由于您对R很熟悉,所以可以进行所有的地统计您需要从那里。我相信您也可以将地统计结果导出回postgis。

通过将SQL与Postgis函数配合使用,您还可以进行各种空间分析,例如叠加操作,距离等。

对于地图绘图,我将使用QGIS(一种开放源GIS软件),该软件可以直接读取postgis(据我所知,这是该项目的最初目标),而即将发布的2.0版将提供许多功能漂亮的地图。

评论


好的建议,但是由于我想自动执行R的所有工作(包括绘图),因此要中断流程,不是吗?

– Nanounanue
13年7月1日在23:26

在这种情况下,如果您愿意的话,只需使用R绘制地图即可。即使这样,让qgis项目布局基于postgis(更新的)数据来准备,它们也会随之更新。我想最终使用R还是QGIS将是一个个人选择。

–亚历山大·内托(Alexandre Neto)
2013年7月1日23:48



感谢您的快速回复,但是如何在Postgis中的表格中使用R作图?

– Nanounanue
13年7月1日在23:50

我对R不太熟悉,也不知道如何使用它绘制矢量数据(就像我说过的那样使用QGIS),如何在R中绘制shapfile?为了从R连接到PostgresSQL,我以前使用过RPostgreSQL。我认为rgdal]。祝好运!

–亚历山大·内托(Alexandre Neto)
13年7月2日在9:11

#3 楼

更新:
从v0.9开始,st_read_db()st_write_db()已失效,您可以简单地使用st_read()st_write()

新引入的sf-package(sp的后代)提供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_db(),您可以上传数据。

评论


这是最简单的解决方案,其中有一个小插图cran.r-project.org/web/packages/sf/vignettes/sf2.html解释了如何使用sf

–塞德里克
19年3月3日在15:17

#4 楼

您可以针对解决方案的每个步骤同时使用所有工具。


如果要进行地静力分析,请使用R。R的软件包更加可靠,可让您更具分析性的结果。您可以基于SQL查询导入数据。
如果要基于逻辑汇总数据,则可以使用PostGIS。您可以回答复杂的查询,例如哪些点在我规定的范围内?但是规模很大。
对于映射,可以使用R或QGIS。 QGIS更为直接,使用R可能会为获得所需的结果而苦苦挣扎。

如果您能为我们提供更多有关问题的详细信息,我们可以为您提供更具体的答案

评论


您能否提供最后一点的示例,我的意思是,如果我想从Postgis中的表格中绘制带有R的地图,该怎么办?

– Nanounanue
13年7月2日,下午3:55

@nanounanue确保:library(“ rgdal”)mydata = readOGR(dsn =“ PG:dbname = ”,layer =“ schema.table”)plot(mydata,axes = TRUE)title(“我的情节”)。

–尼克斯
13年7月2日在9:17



也请看此页面:wiki.intamap.org/index.php/PostGIS

–尼克斯
13年7月2日在9:20

#5 楼

我还将结合使用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


#6 楼

现在有一个RPostGIS软件包,可以使用SQL查询将PostGIS几何导入R中。

#7 楼

您也可以结合使用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")


#8 楼

如果以“ 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个几何。

#9 楼

Geotuple-https://github.com/rhansson/geotuple
是一个Web应用程序,用于连接R-Server和PostGIS(使用RPostgreSQL)