R语言矢量数据空间分析一：入门及rgeo包简介

首先确定一个概念，虾神此次在这里讲的“矢量数据空间分析”并非指的是那些高大上的诸如什么空间聚类、密度分析、插值神马的分析，而且专值的“矢量数据”的一些基础的图形操作，更准确的说，应该是指矢量数据相关的几何特征行为的空间关系和空间操作，也就是下面这些OGC规范里面的东西：



实际上这些都是空间分析里面最基础的一些内容，几乎可以归类到“查询”里面，属性的查询很容易被理解，空间上的查询比属性的稍微复杂那么一点点，所以专门来讲讲。

当然，肯定也是不讲什么射线法、转角法这种基础空间算法，这里专门就讲如何用R语言来实现这些。

实际上如果熟悉诸如ArcGIS一类的软件的话，这些工具都是现成的，可能大家都不觉得需要专门来讲讲怎么去做，Toolbox里面的工具双击就完成了。为什么虾神要专门写一个系列（起码五六章+）来讲这些基础操作呢？

原因就是近期在用R语言的时候，发现属性筛选（查询、处理）的包和方法，在R语言里面俯拾皆是，比如虾神最喜欢的sqldf包，以及使用最广泛的dplyr包，这些包都提供了完善的属性数据处理方法。但是作为空间统计专业户，不知道大家会不会有和虾神一样的疑问：R里面有没有空间筛选的包呢？比如我就需要某条道路左右100米范围内的数据，超过这个范围的，统统不要。

上面那个问题，在ArcGIS里面，那是简单到惨绝人寰……但是要在R语言里面做，怎么办？

所以基于上面的问题，虾神准备专门出一个系列，来讲讲如何在R语言里面做这类分析。

在R语言里面，我们要使用的矢量数据空间分析的包，就是大名鼎鼎的rgeo包。全称是：R for Geometry Engine - Open Source (R语言版开源几何引擎接口）熟悉OGC和GDAL的同学是不是有似曾相识的感觉？没错，这个就是OGC里面那套大名鼎鼎的几何引擎的R语言版本：


官方地址是：http://trac.osgeo.org/geos

具体的说明，大家自行去查询，我这里就不多说了，下面我们来看看R语言怎么来用这个包。

首先当然是安装：

R语言的软件安装是所有语言里面最简单的，只要敲一个命令就是：

install.packages("rgeo")



当然，我在北京，选择的是清华大学的镜像，南方的同学可以选择合肥中科大的镜像。

只要不报错，就OK了，如果安装的时候报错，就得具体问具体分析了，一般来说R语言安装的错误，大部分都与rtools有关，这个大家自己去百度。

当然，你在安装的时候，肯定还会有一堆的依赖包被安装，比如sp包、maptools包啥的，大家安装的时候安按照默认方式来安装就行。

下面做一个测试：

首先在北京范围内，生成一个1000个随机点，并且生成一个点图成：

library(rgeos)
library(sp)
library(maptools)

path <- "E:/data/R/"
crs <-  CRS("+init=epsg:4326")

#加载并且绘制出北京行政区划
beijing <- readShapePoly(paste(path,"beijing.shp",sep =""))
proj4string(beijing) <- crs
plot(beijing,axes=TRUE,col="#FFFFBE")

#生成随机点
pntdf = data.frame(
  lng = runif(1000,min=115.5,max=117.5),
  lat = runif(1000,min=39.5,max=41),
  id = 1:1000,
  flag = rep(0,1000)
)
pnt <- SpatialPointsDataFrame(pntdf,data = pntdf,proj4string=crs)

#把点绘制上去
points(pnt,pch=20,cex=0.8)

运行完了，结果是这样的：

然后再以北京中心tianan -door 生成一个中心点：

ctpnt = SpatialPoints(data.frame(lon=c(116.391),lat=c(39.906)),
proj4string= crs)
points(ctpnt,pch=18,cex=2,col="red")

以这个中心点，绘制一个0.2度缓冲区：
因为这里用的数据是wgs 84的空间参考，单位就是度，所以缓冲区生成的单位只能是度了。

ctbf <- gBuffer(ctpnt,id=1,width = 0.2,quadsegs=10)
ctbfxy <- slot(slot(ctbf@polygons[[1]],"Polygons")[[1]],"coords")
lines(x=ctbfxy[,1],y=ctbfxy[,2],col="blue",lwd=3)

最后来查询一下，这个缓冲区内有哪些点：

（这样写循环，是为了了解方法后面的原理，实际遍历是效率最低的方式，正确的方法应该用R语言推荐的矩阵运算，请关注后续的内容，这种方法在实际工作中，不建议使用，下面的雷同）

for (n in 1:length(pnt)){
  temppnt = SpatialPoints(data.frame(x=pnt@coords[n,1],
                                     y=pnt@coords[n,2]),
                          proj4string= crs)
  if(gIntersects(temppnt,ctbf)){
    pnt@data$flag[n]<- 1
  }
}
selPnt <- subset(pnt,pnt@data$flag==1)
points(selPnt,pch=3,col="red")

打完收工。

从上面这个例子可以看见，R语言也可以做各种空间查询，但是相对来说，不如专业的GIS软件那么方便罢了，至于上面这些语言和方法是啥意思，如何用的，请听下回分解。

最后，扩展一个，来个线要素缓冲，结果如下：


代码如下：

bjline = Lines(list(Line(cbind(c(115.8,116,116.3,117.2),
                               c(41,39.7,40.3,39.8)))),ID = "1")
sline = SpatialLines(list(bjline),proj4string=crs)
lineBF = gBuffer(sline,width = 0.2)
pnt@data$flag[] <-0
for (n in 1:length(pnt)){
  temppnt = SpatialPoints(data.frame(x=pnt@coords[n,1],
                                     y=pnt@coords[n,2]),
                          proj4string= crs)
  if(gIntersects(temppnt,lineBF)){
    pnt@data$flag[n]<- 1
  }
}
selPnt <- subset(pnt,pnt@data$flag==1)


plot(beijing,axes=TRUE,col="#FFFFBE")
points(pnt,pch=20,cex=0.8)
lines(sline,lwd=3,col="blue")
bl = slot(slot(lineBF@polygons[[1]],"Polygons")[[1]],"coords")
lines(x=bl[,1],y=bl[,2],col="red",lwd=2)
points(selPnt,pch=3,col="red")

待续未完。