sf | 创建空间矢量对象及其投影设置

R语言的空间数据主要有sp和sf两种格式。本篇介绍使用sf工具包中的相关函数创建sf格式的空间矢量对象，以及设置和转换矢量对象的投影坐标信息。

1 创建空间矢量对象

1.1 sf对象的组成

sf对象是由两部分信息构成的特殊数据框：

• 一部分为要素属性信息，是一个普通的数据框对象；

• 另一部分为几何信息，是一个list对象，sf工具包称其为sfc (simple feature geometry list-column)，它储存矢量数据中所有要素的几何信息，其中的任意一个元素都代表着一个几何要素，也就是sfg对象 (simple feature geometry)。

使用st_geometry()可以提取sf对象的几何信息

library(sf)
nc <- st_read(system.file("shape/nc.shp", package="sf"))  # 示例数据

# 提取几何信息
# 方法1——使用基本操作符
geom_nc <- nc$geom
# 方法2——使用sf包中的st_geometry函数
geom_nc2 <- st_geometry(nc)

# 删去几何信息
df_nc <- st_drop_geometry(nc)

# 查看数据类型
class(nc)
class(geom_nc)
identical(geom_nc, geom_nc2)
class(geom_nc[[1]]) # sfc对象的任意一个元素是sfg对象
class(df_nc) # 删去几何信息后变成了普通数据框

# 绘图
plot(nc) # sf对象
plot(geom_nc) # sfc对象

# 部分输出结果
> # 查看数据类型
> class(nc)
[1] "sf"         "tbl_df"     "tbl"        "data.frame"
> class(geom_nc)
[1] "sfc_MULTIPOLYGON" "sfc"             
> identical(geom_nc, geom_nc2)
[1] TRUE
> class(geom_nc[[1]])
[1] "XY"           "MULTIPOLYGON" "sfg"  
> class(df_nc) # 删去几何信息后变成了普通数据框
[1] "tbl_df"     "tbl"        "data.frame"

• 比较两幅图可以发现，plot(nc)是根据nc的各个属性进行分类涂色，而plot(geom_nc)则仅仅绘制nc的形状，这是因为nc本身包含两部分，而geom_nc作为sfc对象只包含它的几何信息部分。

1.2 创建单个空间要素 (sfg)

sfc中的文本信息可以使用st_as_text()进行提取，sf称这些文本信息为well-known text (WKT)。

text_geom_nc <- st_as_text(geom_nc)

创建sfg需要先写出WKT，再根据要素类型选择相应函数转化为sfg对象

• st_point

• st_multipoint

• st_linestring

• st_multilinestring

• st_polygon

• st_multipolygon

• st_geometrycollection

par(mfrow = c(2, 4))

# 数据准备：定义若干个点
p1 <- c(1, 1)
p2 <- c(3, 1)
p3 <- c(2, 3)
p4 <- c(1, 4)
p5 <- c(3, 4)
p6 <- c(2, 2)

# 创建单个点要素
point <- st_point(p1)
plot(point)

# 创建单个要素含有多个点，需要将包含点的坐标合并成矩阵
mp <- rbind(p1, p2, p3)

mpoint <- st_multipoint(mp)
plot(mpoint)

# 创建单个线要素，需要将包含点的坐标合并成矩阵
mp <- rbind(p1, p2, p3)
line <- st_linestring(mp)
plot(line)

# 创建单个要素含有多条线，需要将包含线的矩阵合并成list
mp <- rbind(p1, p2, p3)
mp2 <- rbind(p4, p5, p6)
ml <- list(mp, mp2)

mline <- st_multilinestring(ml)
plot(mline)

# 创建单个面要素，首尾点坐标必须相同
mp3 <- rbind(p1, p2, p3, p1)
polygon <- st_polygon(list(mp3))
plot(polygon, col = "grey", border = "red")

mp4 <- rbind(p4, p5, p6, p4)
polygon2 <- st_polygon(list(mp3, mp4))
plot(polygon2, col = "grey", border = "red")

# 创建单个要素包含多个面，使用多个list组成一个list
mpolygon <- st_multipolygon(list(list(mp3), list(mp4)))
plot(mpolygon, col = "grey", border = "red")

# 创建单个要素包含点、线、面中的两类及以上，可以分别先创建单类型要素，再合并成list
mixsf <- st_geometrycollection(list(mpoint, mline, polygon))
plot(mixsf)

1.3 创建多要素几何对象 (sfc)

从集合与元素关系的角度去看，sfg是一个几何元素，而sfc是由若干个sfg组成的集合，因此创建sfc对象需要先创建若干个sfg对象，再使用st_sfc()进行合并。

以点要素为例：

# 数据准备：定义若干个点
p1 <- c(1, 1)
p2 <- c(3, 1)
point1 <- st_point(p1)
point2 <- st_point(p2)

point <- st_sfc(point1, point2)
plot(point)

1.4 创建sf点对象

手动创建sfg和sfc在大多数情况下并无太大实际意义。现实操作常常是从外部导入含有属性和地理坐标信息的表格 (excel文件、 csv等)，然后将其直接转为矢量数据，使用的函数是st_sf()或st_as_sf()。

library(tidyverse)

# 生成示例数据，数据中包含地理坐标信息
dta01 <- data.frame(
  ID = 1 : 10,
  x = runif(10, min = 116, max = 117),
  y = runif(10, min = 30, max = 33)
)

# 指定地理坐标信息，x表示经度，y表示纬度
dta01<- st_as_sf(dta01, coords = c(x = "x", y = "y"))

当需要使用已有空间数据的地理信息时，可以先提取出WKT信息，再调用st_as_sf()中的wkt参数

# 使用已有文件的空间信息
dta02 <- data.frame(ID = 1 : 10) # 数据中没有地理坐标信息
dta02$geom <- st_as_text(st_geometry(dta01)) # 提取dta01的WKT信息
dta02 <- st_as_sf(dta02, wkt = "geom")

## 以上代码也可以使用管道操作符书写
data.frame(ID = 1 : 10) %>%
  mutate(geom = st_as_text(st_geometry(dta01))) %>%
  st_as_sf(wkt = "geom") -> dta02

2 设置投影坐标系

当一个矢量数据本身存在地理坐标或投影坐标时，通过st_read()读入R空间后仍然会保留原先的坐标系统。

使用st_crs()可以读取矢量数据的坐标系统

crs_nc <- st_crs(nc)
crs_nc <- st_crs(nc)

使用st_is_longlat()查看对象是否为经纬度坐标

st_is_longlat(nc)

在不改变坐标值的前提下，修改坐标系统，可以使用st_crs<-或st_set_crs()

# 方法1
st_crs(nc) <- st_crs(nc) 
# 方法2
nc <- st_set_crs(nc, st_crs(nc))

• 以上代码是将nc的坐标系统强加给nc，但是由于二者的坐标系统本身并不一致，这一操作实际上毁坏了nc原有的坐标系统。更多地，我们是在为从外部导入含有坐标的表格设置坐标系统时才使用这两个函数。

R中的坐标系统采用的是EPSG代码，比如EPSG等于4326时代表的是GCS_WGS_1984地理坐标系，4214代表的是GCS_Beijing_1954地理坐标系。网站http://epsg.io/可以通过输入EPSG代码查询对应的坐标系统，网址https://blog.csdn.net/qq_41441896/article/details/104525296列举了一些常用的地理和投影坐标系。除了使用固定的EPSG代码外，还可以通过编辑WKT来自定义坐标系统。

例如，已有表格使用的是GCS_WGS_1984地理坐标，在导入R空间后，可以做如下设置：

data.frame(
  ID = 1 : 10,
  x = runif(10, min = 116, max = 117),
  y = runif(10, min = 30, max = 33)
) %>% 
  st_as_sf(coords = c(x = "x", y = "y")) %>%
  st_set_crs(4326) -> dta

• 这样原先普通的数据框就变成了拥有坐标系统的空间矢量对象了。

使用st_transform()可以进行坐标变换

# 使dta的坐标系统与nc保持一致
dta %>%
  st_transform(st_crs(nc))-> dta

• st_transform()与st_set_crs()的区别在于，它是通过改变对象的坐标值去匹配目标坐标系统从而改变坐标系统的。

 

 

---