R语言中的地理/投影坐标系统（下）[翻译]
❝
「译者注」：在原文的本部分出现了许多世界地图，并且地图涉及到了我国领土边界的标识，但是这种标识是错误的。为了避免这一错误并尽可能地保留原文面貌，在确实需要标记区域边界时，删去了亚洲国家和地区；没有必要标识边界时则省略边界。所有改动处均有注明。
❞
为了与上篇衔接，可以先运行如下代码：
library(sf)
library(raster)
load(url("https://github.com/mgimond/Spatial/raw/main/Data/Sample1.RData"))7 「坐标系统转换」
上节展示的是如何定义（）或修改（）坐标系统。本节来展示如何转换坐标系统，即将空间对象的坐标映射到其他坐标系统上去。这个过程会产生新的坐标取值。
例如，将s.sf的坐标系统转换成WGS 1984地理（经纬度）坐标系统，我们将使用()函数：
s.sf.gcs <- st_transform(s.sf, "+proj=longlat +datum=WGS84")
st_crs(s.sf.gcs)## 译者注：输出结果略使用等价于上面proj4字符串的EPSG代码：
s.sf.gcs <- st_transform(s.sf, 4326)
st_crs(s.sf.gcs)## 译者注：输出结果略如果要转换栅格数据的坐标系统，使用的是()函数：
elev.r.gcs <- projectRaster(elev.r, crs="+proj=longlat +datum=WGS84")
crs(elev.r.gcs)## CRS arguments: +proj=longlat +datum=WGS84 +no_defs如果要使用EPSG代码，需要采用+init=EPSG:...语句：
elev.r.gcs <- projectRaster(elev.r, crs="+init=EPSG:4326")
crs(elev.r.gcs) ## CRS arguments: +proj=longlat +datum=WGS84 +no_defs当需要将图层叠加到基于、Bing或等地图服务平台的web图层时，通常就会使用地理坐标系统（即使这些地图平台最终也会投影到投影坐标系通——一般是网络墨卡托投影（Web ））。为了检查s.sf.gcs的投影是否转换正确，我们将使用工具包把它叠加到的顶部图层上：
library(leaflet)
leaflet(s.sf.gcs) %>% addPolygons() %>% addTiles()❝
译者注：关于工具包，可查看如下推文：
❞
下面我们将使用tmap工具包的世界地图数据探索其他转换过程。
library(tmap)
data(World) # 数据格式为sf# 检查当前坐标系统
st_crs(World)## 译者注：输出结果略下列代码将世界地图转换成一个自定义的等距方位 ( ) 投影，中心点为0度经度、0度纬度。这里我们使用的是proj4语法。
World.ae <- st_transform(World, "+proj=aeqd +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")让我们检查对象转换后的坐标系统：
st_crs(World.ae) ## 译者注：输出结果略绘制地图：
tm_shape(World.ae) + tm_fill()接下来的代码会将坐标系统转换成中心点在美国缅因州（西经69.8度，北纬44.5度）的等距方位投影：
World.aemaine <- st_transform(World, "+proj=aeqd +lat_0=44.5 +lon_0=-69.8 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")tm_shape(World.aemaine) + tm_fill()下面代码转换后是罗宾逊投影（ ）：
World.robin <- st_transform(World, "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")tm_shape(World.robin) + tm_fill()下面转换后是正弦投影（ ）：
World.sin <- st_transform(World, "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")tm_shape(World.sin) + tm_fill()转换后是墨卡托投影（ ）：
World.mercator <- st_transform(World, "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")tm_shape(World.mercator) + tm_fill()「一个转换失败的例子」
当坐标系统定义的切线或点使多边形面（）沿180度子午线割裂时，空间数据转换过程就会出问题。例如，将墨卡托投影的中心点设置西经69度时就会产生如下地图：
World.mercator2 <- st_transform(World, "+proj=merc +lon_0=-69 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")library(tidyverse)
World.mercator2_excluded_asia <- filter(World.mercator2, continent != "Asia")
tm_shape(World.mercator2_excluded_asia) + tm_borders()❝ ❞
多边形面被扯裂，但是R并不知道如何把它拼接起来。
一个解决方法是使用工具包中的s()函数。这个函数可以将多边形面沿指定经度分割，但它要求空间对象的结构必须为*形式（译者注：即sp格式），坐标系统为地理坐标系统（经纬度）。下面代码是一个流程示例：
library(maptools)# 转换成经纬度坐标
wld.ll <- st_transform(World, "+proj=longlat +datum=WGS84 +no_defs")
wld.ll_excluded_asia <- filter(wld.ll, continent != "Asia")# 转换成sp格式，再沿指定经度分割（本例为东经111度）
wld.sp <- nowrapSpatialPolygons(as(wld.ll, "Spatial"), offset = 111)
wld.sp_excluded_asia <- nowrapSpatialPolygons(as(wld.ll_excluded_asia, "Spatial"),offset = 111)# 重新转换成sf格式，再以西经69度为中心进行投影
# 然后绘图
wld.sf <- st_as_sf(wld.sp)
wld.sf_excluded_asia <- st_as_sf(wld.sp_excluded_asia)
wld.merc2.sf <- st_transform(wld.sf, "+proj=merc +lon_0=-69 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
wld.merc2.sf_excluded_asia <- st_transform(wld.sf_excluded_asia, "+proj=merc +lon_0=-69 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
tm_shape(wld.merc2.sf_excluded_asia) + tm_borders()❝ ❞
你可以注意到，南极大陆在转换中变形了。这种情况下，可以在使用s()函数处理前先从数据中删去南极大陆或者使用其他合适的投影。在下面的例子中，我们使用罗宾逊投影代替墨卡托投影：
wld.rob.sf <- st_transform(wld.sf, "+proj=robin +lon_0=-69 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")tm_shape(wld.rob.sf) + tm_fill()❝ ❞
使用s()函数的一个缺点是会丧失矢量数据的属性信息：
head(data.frame(wld.rob.sf), 4)一个（但不完美）的办法是在投影前先生成矢量数据的质心，然后再使用空间连接方法将质心的属性数据添加到投影后的数据上。
# 提取面的质心
pt <- st_centroid(World, of_largest_polygon = TRUE)# 将质心数据转换到罗宾逊投影
pt.rob <- st_transform(pt, "+proj=robin +lon_0=-69 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")# 进行空间连接（将质心的数据值添加到wld.rob.sf上去）
wld.rob.df.sf <- st_join(wld.rob.sf, pt.rob, join = st_contains)# 输出地图
wld.rob.df.sf_excluded_asia <- filter(wld.rob.df.sf, continent != "Asia" | is.na(continent))
tm_shape(wld.rob.df.sf_excluded_asia) + tm_polygons(col = "pop_est_dens", style = "quantile") +tm_legend(outside = TRUE)❝ ❞
虽然这种方法对大多数区域适用，但是少部分区域，如挪威显示是缺失值。为了探究原因，我们使用投影前的数据放大挪威并叠加质心数据。
library(dplyr)# 提取挪威和瑞典的范围
nor.bb <- World %>% filter(name == "Norway" | name == "Sweden") %>% st_bbox()# 绘制区域放大图，并添加质心作为参考
tm_shape(World, bbox = nor.bb) + tm_polygons(col = "pop_est_dens", style = "quantile") +tm_shape(pt) + tm_dots() +tm_text("name", just = "left", xmod = 0.5, size = 0.8) +tm_legend(outside = TRUE)你会注意到，挪威的质心落到了瑞典境内。这是因为挪威的形态呈“C”状，而()函数在计算质心时使用的是地理中心（ ）而不是“质量中心”（ of mass）。
另一个解决办法是使用()函数，它将会在每个多边形面内产生一个点。
pt <- st_point_on_surface(World)pt.rob <- st_transform(pt, "+proj=robin +lon_0=-69 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
wld.rob.df.sf <- st_join(wld.rob.sf, pt.rob, join = st_contains)
wld.rob.df.sf_excluded_asia <- filter(wld.rob.df.sf, continent != "Asia" | is.na(continent))
tm_shape(wld.rob.df.sf_excluded_asia) + tm_polygons(col = "pop_est_dens", style = "quantile") +tm_legend(outside = TRUE)❝ ❞ 8 「关于包含的说明」
从理论上讲，如果一个点位于一个闭合的区域内，那么不管如何进行投影变换，该点应该始终位于这个区域内，然而实际上并非总是如此。这是因为投影变换作用的是构成区域边界（线）的点，而非边界本身。如果一个点位于区域内但非常接近边界，那么在投影转换后它可能会“移”到边界的另一边从而位于区域外。在下列例子中，面和点都处在米勒坐标系统（ ）下，并且点位于面内。
# 定义投影
miller <- "+proj=mill +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
lambert <- "+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"# Subset the World data layer
wld.mil <- World %>% filter( iso_a3 == "CAN" | iso_a3 == "USA") %>% st_transform(miller)# Create polygon and point layers in the Miller projection
sf1 <- st_sfc( st_polygon(list(cbind(c(-13340256,-13340256,-6661069, -6661069, -13340256),c(7713751, 5326023, 5326023,7713751, 7713751 )))), crs = miller) pt1 <- st_sfc(st_multipoint(rbind(c(-11688500,7633570), c(-11688500,5375780),c(-10018800,7633570), c(-10018800,5375780),c(-8348960,7633570), c(-8348960,5375780))), crs = miller)
pt1 <- st_cast(pt1, "POINT") # Create single part points# Plot the data layers in their native projection
tm_shape(wld.mil) +tm_fill(col = "grey") + tm_graticules(x = c(-60,-80,-100,-120,-140), y = c(30,45, 60), labels.col = "white", col = "grey90") +tm_shape(sf1) + tm_polygons("red", alpha = 0.5, border.col = "yellow") +tm_shape(pt1) + tm_dots(size=0.2)
这些点非常接近边界，但还都在多边形面内。为了证明这一点，我们使用()函数：
st_contains(sf1, pt1)## Sparse geometry binary predicate list of length 1, where the predicate
## was `contains'
## 1: 1, 2, 3, 4, 5, 6如预期的那样，六个点均包含在内。
现在，让我们把数据转成兰勃特等角投影（ ）。
# Transform the data
wld.lam <- st_transform(wld.mil, lambert)
pt1.lam <- st_transform(pt1, lambert)
sf1.lam <- st_transform(sf1, lambert)# Plot the data in the Lambert coordinate system
tm_shape(wld.lam) +tm_fill(col = "grey") + tm_graticules(x = c(-60,-80,-100,-120,-140), y = c(30,45,60), labels.col = "white", col = "grey90") +tm_shape(sf1.lam) + tm_polygons("red", alpha = 0.5, border.col = "yellow") +tm_shape(pt1.lam) + tm_dots(size=0.2)只有三个点还在面内。我们可以再用()函数来验证：
st_contains(sf1.lam, pt1.lam)## Sparse geometry binary predicate list of length 1, where the predicate
## was `contains'
## 1: 1, 3, 5为了解决这个问题，我们应当让构成多边形边界的顶点更密集。顶点密度取决于保持这种包含关系所需要的分辨率大小，而最好的方法是试验。
我们使用()函数设置边界顶点的间隔距离为1 km（1000 m）：
# Add vertices every 1000 meters along the polygon's line segments
sf2 <- st_segmentize(sf1, 1000)# Transform the newly densified polygon layer
sf2.lam <- st_transform(sf2, lambert)# Plot the data
tm_shape(wld.lam) + tm_fill(col = "grey") + tm_graticules(x = c(-60,-80,-100,-120,-140), y = c(30,45, 60), labels.col = "white", col = "grey90") +tm_shape(sf2.lam) + tm_polygons("red", alpha = 0.5, border.col = "yellow") +tm_shape(pt1.lam) + tm_dots(size = 0.2)现在，投影转换后所有的点仍然都位于面内。我们可以通过下面方法检验：
st_contains(sf2.lam, pt1.lam)## Sparse geometry binary predicate list of length 1, where the predicate
## was `contains'
## 1: 1, 2, 3, 4, 5, 69 「创建天梭指示圆（ ）」
大多数投影都会扭曲一部分空间属性，特别是面积和形状。可视化这种扭曲性质的一个好方法就是创建大地圆（ ）。
首先，在经纬度坐标系统下定义作为圆心的点图层：
tissot.pt <- st_sfc(st_multipoint(rbind(c(-60,30), c(-60,45), c(-60,60),c(-80,30), c(-80,45), c(-80,60),c(-100,30), c(-100,45), c(-100,60),c(-120,30), c(-120,45), c(-120,60) )), crs = "+proj=longlat")
# Create single part points
tissot.pt <- st_cast(tissot.pt, "POINT")然后，我们使用工具包基于这些点构建大地圆：
library(geosphere)cr.pt <- list() # Create an empty list# Loop through each point in tissot.pt and generate 360 vertices at 300 km
# from each point in all directions at 1 degree increment. These vertices
# will be used to approximate the Tissot circles
for (i in 1:length(tissot.pt)){cr.pt[[i]] <- list(destPoint( as(tissot.pt[i], "Spatial"), b = seq(0,360,1), d = 300000) )
}# Create a closed polygon from the previously generated vertices
tissot.sfc <- st_cast(st_sfc(st_multipolygon(cr.pt ),crs = "+proj=longlat"), "POLYGON" )我们将通过计算这些圆的面积去验证它们是大地圆。使用的函数是sf工具包中的()，它在经纬度坐标系统下会计算地测面积（ area）。
tissot.sf <- st_sf(geoArea = st_area(tissot.sfc), tissot.sfc )在我们的例子中，真实面积应该是平方米。让我们计算输出结果的误差。数值将以小数的形式展示：
((pi*300000^2) - as.vector(tissot.sf$geoArea))/(pi*300000^2)## [1] -0.0008937164 0.0024530577 0.0057943110 -0.0008937164 0.0024530577
## [6] 0.0057943110 -0.0008937164 0.0024530577 0.0057943110 -0.0008937164
## [11] 0.0024530577 0.0057943110在所有情况下，误差均小于0.1%。这个误差主要是由圆参数的离散化引起的（「译者注」：该处的意思应该是指矢量数据中的圆实际是由离散的点构成的）。
让我们看看一些常见的坐标系统所带来的空间扭曲。我们从墨卡托投影开始：
# Transform geodesic circles and compute area error as a percentage
tissot.merc <- st_transform(tissot.sf, "+proj=merc +ellps=WGS84")
tissot.merc$area_err <- round((st_area(tissot.merc, tissot.merc$geoArea)) / tissot.merc$geoArea * 100, 2)# Plot the map
tm_shape(World, bbox = st_bbox(tissot.merc), projection = st_crs(tissot.merc)) + tm_borders() + tm_shape(tissot.merc) + tm_polygons(col = "grey", border.col = "red", alpha = 0.3) + tm_graticules(x = c(-60,-80,-100,-120,-140), y = c(30,45,60),labels.col = "white", col = "grey80") +tm_text("area_err", size =.8, alpha = 0.8, col = "blue")墨卡托投影很好地保留了对象的形状，但是面积扭曲较严重。
接下来，我们探究以北纬45度、西经100度为中心的兰勃特等面积投影（ equal area ）：
# Transform geodesic circles and compute area error as a percentage
tissot.laea <- st_transform(tissot.sf, "+proj=laea +lat_0=45 +lon_0=-100 +ellps=WGS84")
tissot.laea$area_err <- round((st_area(tissot.laea ) - tissot.laea$geoArea) /tissot.laea$geoArea * 100, 2)# Plot the map
tm_shape(World, bbox = st_bbox(tissot.laea), projection = st_crs(tissot.laea)) + tm_borders() + tm_shape(tissot.laea) + tm_polygons(col="grey", border.col = "red", alpha = 0.3) + tm_graticules(x=c(-60,-80,-100, -120, -140), y = c(30,45, 60),labels.col = "white", col="grey80") +tm_text("area_err", size=.8, alpha=0.8, col="blue")遍及48个州，面积误差都几乎为0。但是注意，从投影中心到外围圆的形状会变得扭曲。
再接下来，我们将探究罗宾森投影：
# Transform geodesic circles and compute area error as a percentage
tissot.robin <- st_transform(tissot.sf, "+proj=robin +ellps=WGS84")
tissot.robin$area_err <- round((st_area(tissot.robin ) - tissot.robin$geoArea) /tissot.robin$geoArea * 100, 2)# Plot the map
tm_shape(World, bbox = st_bbox(tissot.robin), projection = st_crs(tissot.robin)) + tm_borders() + tm_shape(tissot.robin) + tm_polygons(col = "grey", border.col = "red", alpha = 0.3) + tm_graticules(x =c(-60,-80,-100, -120, -140), y = c(30,45, 60),labels.col = "white", col = "grey80") +tm_text("area_err", size=.8, alpha = 0.8, col = "blue")北美大陆的形状和面积都发生了肉眼可见的扭曲。
