我有一个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以进行转换,过滤或限制。不过,这应该可以帮助您入门。
#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 =
–尼克斯
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)
评论
很好的答案,但是如何在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