Spatial data in R

• 数据结构
• 基础作图
• 构建缓冲区
• 象限分颜色表示

---

加载数据

第一部分我们已经看了R的语法并安装了必要的包，我们可以开始查看一些真实的空间数据。

下载

本教程的数据下载链接：Spatial data，打开链接后，依次点击右手边的Clone or download、Download ZIP，下载完成后解压即可。

项目文件

在RStudio中的所有操作建议在项目文件中进行，刚才下载的文件里已经包含了本教程的下载文件。打开步骤：
打开RStudio，依次点击左上角的File -> Open File…，找到上一步下载并解压后的文件夹，双击打开Creating-maps-in-R.Rproj项目文件。如图2所示：

数据加载

数据加载之前最好通过其他软件如LibreOffice(免费的办公软件)预览一下数据结构。然后在RStudio输入以下代码：

library(rgdal)  # load package ：Geospatial Data Abstraction Library
lnd <- readOGR(dsn = "data", layer = "london_sport") # read data

**readOGR():**读取shapfile文件,并在R中赋予它一个新的对象名：lnd
第一个参数，dsn = shapfile文件所在的文件夹名
第二个参数，layer = shapfile文件名，省略扩展名
文件名均用字符串形式表示(即:”文件名”)
R有一个默认参数执行顺序：readOGR(“data”, “london_sport”)效果同上
至于如何加载其他不同类型的空间数据，可以通过?readOGR查看帮助文档。

数据结构

属性槽

R中的空间对象由不同的**属性槽**组成，比较重要的：
@: 获取对象的属性槽
@data: 数据槽，也就是数据属性表。
@polygons：多边形槽，顶点的空间位置。

基本函数

以2001年的伦敦自治市的体育数据为例:

head函数

head(lnd@data, n = 2)  # show the first two lines of data
head(lnd@polygons[[1]]@Polygons[[1]]@coords, 3)
mean(lnd$Partic_Per)

1.head(lnd@data):默认输出数据属性表的前六行
1.head(lnd@data, n = 2):输出数据属性表的前两行
2.lnd@polygons[[1]]: 第一个多边形
2.lnd@polygons[[1]]@Polygons[[1]]: 第一个多边形的**第一个多边形形状(通常也只有一个)**
2.lnd@polygons[[1]]@Polygons[[1]]@coords: 第一个多边形的第一个多边形形状的**坐标**
3.$: 获取数据属性表中的变量
3.Partic_Per :sports participation per 100 people

sapply函数

a <- 1:10
sqrt_a <- sapply(a,sqrt)
double_a <- sapply(a,function(x) 2*x)
sapply(lnd@data, class) # check the classes of all the variables
lnd$Pop_2001 <- as.numeric(as.character(lnd$Pop_2001))

1.生成整数：1-10
2.求a的平方根，sqrt_a的长度与a相等
3.求2a，sqrt_a的长度与a相等
4.查看变量属性
5.将因子变量转化为数值型变量

基础作图

底图

plot(lnd)

条件查询

lnd@data[1:2, 1:3]
lnd@data[lnd$Partic_Per < 15, ]

[1:2, 1:3]：显示数据的前两行，前三列
[lnd$Partic_Per < 15, ]：显示数据中变量Partic_Per小于15的所有行，**逗号**后为空：默认显示：所有列

作图

sel <- lnd$Partic_Per > 20 & lnd$Partic_Per < 25
head(sel)
plot(lnd[sel, ])

plot(lnd, col = "lightgrey") # plot the london_sport object
sel <- lnd$Partic_Per > 25
plot(lnd[ sel, ], col = "turquoise", add = TRUE) # add selected zones to map

col = “lightgrey”: 浅灰色
col = “turquoise”: 绿宝石色
add = TRUE: 在底图**plot(lnd, col = "lightgrey")**上增加新的图层

构建缓冲区

选择地理中心位于伦敦内部地理中心10公里内的所有区域

质心和缓冲区

library(rgeos)  # Geometry Engine - Open Source (GEOS) 
plot(lnd, col = "grey")
cent_lnd <- gCentroid(lnd[lnd$name == "City of London",]) 
points(cent_lnd, cex = 3)
lnd_buffer <- gBuffer(spgeom = cent_lnd, width = 10000) 

2:绘制底图
3:寻找质心
4:查看质心位置
5:构建10km的缓冲区

方法一：相交区域

lnd_central <- lnd[lnd_buffer,] # the selection is too big!
plot(lnd_central, col = "lightblue", add = T)
plot(lnd_buffer, add = T) # some areas just touch the buffer

1:选取相交区域
2:绘图
3:绘制缓冲线

方法二：内部缓冲点

plot(lnd, col = "grey")
lnd_cents <- SpatialPoints(coordinates(lnd),proj4string = CRS(proj4string(lnd))) # create spatialpoints
sel <- lnd_cents[lnd_buffer,] # select points inside buffer
points(sel) # show where the points are located
lnd_central <- lnd[sel,] # select zones intersecting w. sel
plot(lnd_central, add = T, col = "lightslateblue", border = "grey")
plot(lnd_buffer, add = T, border = "green", lwd = 2)

1:由于方法一修改了底图，在此重新绘制底图
2:生成空间点
3:创建select:选取缓冲区内部的点
4:查询符合条件的点
5:选取相交区域
6:绘图
7:绘制缓冲线

图上增加文本

text(coordinates(cent_lnd), "Central\nLondon")

结果如图4所示：

象限分颜色表示

获取经纬度

lng <- coordinates(gCentroid(lnd))[[1]]
lat <- coordinates(gCentroid(lnd))[[2]]

coordinates(lnd): 返回两列数：第一列**经度，第二列纬度**。
2.lng: 获取经度
3.lat: 获取纬度
这个教程来自github上的一个开源项目，原项目的经纬度写反了，我发现了这个错误。

方向布尔值(TRUE/FALSE)

east <- sapply(coordinates(lnd)[,1], function(x) x > lng)
west <- sapply(coordinates(lnd)[,1], function(x) x < lng)
north <- sapply(coordinates(lnd)[,2], function(x) x > lat)
south <- sapply(coordinates(lnd)[,2], function(x) x < lat)

挑战

1.以质心为标准，分颜色表示四个象限

• 提示1：

plot(lnd)
plot(lnd[east & north,],add = TRUE, col = "red" )
llgridlines(lnd, lty= 3, side ="EN", offset = -0.5)

2.分别将四个象限的多边形合并，即只留下四个多边形

• 提示2：

london = gUnaryUnion(lnd) # rgeos
plot(london)

参考结果如图6所示：