Creating and Manipulating Spatial Data

除了可视化和访问之外，GIS还必须能够创建和修改空间数据。 R的空间包为处理和创建空间数据提供了非常广泛和强大的功能套件。
复制、联结/剪裁是GIS的基础操作，因此在本节中我们将探讨如何在R中进行这些操作。
首先，我们将非空间数据加入到空间数据中，以便映射。 最后，我们将介绍空间连接，从而基于空间位置组合来自两个空间对象的信息。

• 创建空间数据
• 设置空间参考系
• 属性连接
• 剪切和空间连接

---

创建空间数据

创建矢量对象

vec <- vector("numeric", 3)
class(vec)

1.第一个参数：数据类型；第二个参数：数据长度
2.返回数据类型

创建数据帧对象

df <- data.frame(x = 1:3, y = c(1/2, 2/3, 3/4))
class(df)

Data Frame: 相当于R中的Excel表
以上均为非空间数据

创建空间点对象

sp1 <- SpatialPoints(coords = df) # input：a numeric matrix or data.frame
class(sp1)

SpatialPoints():来自SP安装包，输入必须为数值矩阵或数据帧。

增加数据

spdf <- SpatialPointsDataFrame(sp1, data = df)  # adding data from df
class(spdf)

1.将非空间数据(df)增加到空间点对象(sp1)中，并命名为spdf

在现实情况下，与上面创建的空间对象不同，您的空间数据将从外部创建的文件读取，例如，使用readOGR()。大多数空间数据都带有关联的CRS。

投影：设置和转换空间参考系

每个空间对象都应有一个坐标参考系CRS。
比较常见的：WGS84，是为GPS全球定位系统使用而建立的坐标系统。

设置CRS

proj4string(lnd) <- NA_character_ # remove CRS information from lnd
proj4string(lnd) <- CRS("+init=epsg:27700") # assign a new CRS

从lnd中移除坐标系
设置新的坐标系（此坐标系为已知坐标系）

转换坐标系

library(rgdal)
EPSG <- make_EPSG() # create data frame of available EPSG codes
EPSG[grepl("WGS 84$", EPSG$note), ] # search for WGS 84 code
lnd84 <- spTransform(lnd, CRS("+init=epsg:4326")) # reproject

保存坐标系

saveRDS(object = lnd84, file = "data/lnd84.Rds")

saveRDS():第一个参数：对象名称；第二个参数：对象存储位置。

移除中间对象

rm(lnd84) # remove the lnd object
# we will load it back in later with readRDS(file = "data/lnd84.Rds")

属性连接

读取矢量数据

library(rgdal) # ensure rgdal is loaded
# Create new object called "lnd" from "london_sport" shapefile
lnd <- readOGR(dsn = "data", "london_sport")
plot(lnd) # plot the lnd object (not shown)
nrow(lnd) # return the number of rows (not shown)

select theft

crime_data <- read.csv("data/mps-recordedcrime-borough.csv",stringsAsFactors = FALSE)
head(crime_data$CrimeType)
crime_theft <- crime_data[crime_data$CrimeType == "Theft & Handling", ]
head(crime_theft, 2)

.CSV格式的数据在R中读取后是数据帧格式

aggregate

crime_ag <- aggregate(CrimeCount ~ Borough, data = crime_theft, sum)
head(crime_ag, 2)

aggregate()三个参数:汇总规则，数据帧名称，函数
The ~ symbol means “by”: 按区名(Borough)汇总CrimeCount变量

高级查询 %in%

lnd$name %in% crime_ag$Borough
lnd$name[!lnd$name %in% crime_ag$Borough]
crime_ag$Borough[!crime_ag$Borough%in% lnd$name]
crime_ag[crime_ag$Borough == "NULL",]

a %in% b: 返回一串逻辑布尔值：若a1在b1中，返回TRUE；反之，返回False。以此类推…
a[!a %in% b]: 返回属于a但不属于b的具体值。

连接空间数据和非空间数据

library(dplyr)
lnd@data <- left_join(lnd@data, crime_ag, by = c('name' = 'Borough'))
lnd$CrimeCount

left_join():返回x中的所有行和x与y的所有列
属于x不属于y的行 返回NA

快速专题地图

library(tmap)
qtm(lnd, "CrimeCount") # input:shape object
lnd$Thefts  <- lnd$CrimeCount / 10000
qtm(lnd, "Thefts",style = "bw")
qtm(lnd, "Thefts",format = "World")

qtm():第一个参数输入为shape object, 包括：
SpatialPolygons(DataFrame)
SpatialPoints(DataFrame)
SpatialLines(DataFrame)
SpatialGrid(DataFrame)
SpatialPixels(DataFrame)
RasterLayer, RasterStack, or RasterBrick
第二个参数输入为：填充颜色
结果如图7所示。

tmap

library(tmap)
lnd_wgs = spTransform(lnd, CRS("+init=epsg:4326"))
osm_tiles = tmaptools::read_osm(bbox(lnd_wgs))
lnd_wgs$Thefts <- lnd$CrimeCount / 10000 
tm_shape(osm_tiles) +
  tm_raster() +
  tm_shape(lnd_wgs) +
  tm_fill("Thefts", fill.title = "Thefts\n(10000)", scale = 0.8, alpha = 0.5) +
  tm_layout(legend.position = c(0.89,0.02))

read_osm()读取开放的街道地图数据。1.读取OSM文件返回栅格数据。2.查询矢量化的OSM数据返回为空间多边形、线或点。

剪切和空间连接

除了属性（例如通过自治市镇名称）连接，还可以在R中进行空间连接，以伦敦自治市交通基础设施点为例。

bbox()

x <- c(1,1,1,2,2,2,3,3,3)
y <- c(1,2,3,1,2,3,1,2,3)
xy <- cbind(x,y)
S <- SpatialPoints(xy)
S@coords
bbox(S)

bbox():从空间数据中检索空间边界框

读取数据

library(rgdal)
stations <- readOGR(dsn = "data", layer = "lnd-stns")
# stations = read_shape("data/lnd-stns.shp") # from tmap
proj4string(stations) # this is the full geographical detail.
proj4string(lnd) # what's the coordinate reference system (CRS)
bbox(stations)
bbox(lnd)

2.读取车站数据
4.5. proj4string()的结果表明：车站的坐标参考系与lnd不同，因此需要转换坐标系。
6.7. bbox():返回空间数据的空间边界框，即经纬度的最大值与最小值。

重新投影

stations <- spTransform(stations, CRS(proj4string(lnd)))
plot(lnd) # plot London 
points(stations) # overlay the station points

车站的空间范围大于lnd，因此需要剪切。

剪切

stations11 <- stations[lnd, ] 
plot(lnd)
points(stations11)
# save(lnd, file="data/lnd.RData")
# save(stations, file="data/stations.RData")

1.选取lnd内部的车站空间点数据
2.绘图
3.4.保存数据
除了[]方法外，还有sp:over和gIntersects方法实现剪切。