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

待续未完。

转载自:https://blog.csdn.net/allenlu2008/article/details/70338447

You may also like...