R语言中的地理/投影坐标系统(下)[翻译]
原文链接:https://mgimond.github.io/Spatial/coordinate-systems-in-r.html。
译文分上、下两篇,这里为下篇。
❝
「译者注」:在原文的本部分出现了许多世界地图,并且地图涉及到了我国领土边界的标识,但是这种标识是错误的。为了避免这一错误并尽可能地保留原文面貌,在确实需要标记区域边界时,删去了亚洲国家和地区;没有必要标识边界时则省略边界。所有改动处均有注明。
❞
为了与上篇衔接,可以先运行如下代码:
library(sf)
library(raster)
load(url("https://github.com/mgimond/Spatial/raw/main/Data/Sample1.RData"))
7 「坐标系统转换」
上节展示的是如何定义(define)或修改(modify)坐标系统。本节来展示如何转换坐标系统,即将空间对象的坐标映射到其他坐标系统上去。这个过程会产生新的坐标取值。
例如,将s.sf
的坐标系统转换成WGS 1984地理(经纬度)坐标系统,我们将使用st_transform()
函数:
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)## 译者注:输出结果略
如果要转换栅格数据的坐标系统,使用的是projectRaster()
函数:
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
当需要将图层叠加到基于Google、Bing或OpenStreetMap等地图服务平台的web图层时,通常就会使用地理坐标系统(即使这些地图平台最终也会投影到投影坐标系通——一般是网络墨卡托投影(Web Mercator projection))。为了检查s.sf.gcs
的投影是否转换正确,我们将使用leaflet
工具包把它叠加到OpenStreetMap的顶部图层上:
library(leaflet)
leaflet(s.sf.gcs) %>% addPolygons() %>% addTiles()
❝
译者注:关于
leaflet
工具包,可查看如下推文:
leaflet |(1)在R语言中导入高德地图
leaflet |(2)如何往地图上添加一只冰墩墩(添加各类要素)
❞
下面我们将使用tmap
工具包的世界地图数据探索其他转换过程。
library(tmap)
data(World) # 数据格式为sf# 检查当前坐标系统
st_crs(World)## 译者注:输出结果略
下列代码将世界地图转换成一个自定义的等距方位 (Azimuthal equidistant) 投影,中心点为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()
下面代码转换后是罗宾逊投影(Robinson projection):
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()
下面转换后是正弦投影(sinusoidal projection):
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()
转换后是墨卡托投影(Mercator projection):
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()
「一个转换失败的例子」
当坐标系统定义的切线或点使多边形面(polygon)沿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并不知道如何把它拼接起来。
一个解决方法是使用maptools
工具包中的nowrapSpatialPolygons()
函数。这个函数可以将多边形面沿指定经度分割,但它要求空间对象的结构必须为Spatial*
形式(译者注:即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()
❝
译者注:此处删去了亚洲国家和地区。
❞
你可以注意到,南极大陆在转换中变形了。这种情况下,可以在使用nowrapSpatialPolygons()
函数处理前先从数据中删去南极大陆或者使用其他合适的投影。在下面的例子中,我们使用罗宾逊投影代替墨卡托投影:
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()
❝
译者注:原文此处显示了国家和地区的边界。这里使用
tm_fill()
代替原文的tm_borders()
函数以隐藏边界。❞
使用nowrapSpatialPolygons()
函数的一个缺点是会丧失矢量数据的属性信息:
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”状,而st_centroid()
函数在计算质心时使用的是地理中心(geometric center)而不是“质量中心”(center of mass)。
另一个解决办法是使用st_point_on_surface()
函数,它将会在每个多边形面内产生一个点。
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 coordinate system)下,并且点位于面内。
# 定义投影
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()
函数:
st_contains(sf1, pt1)## Sparse geometry binary predicate list of length 1, where the predicate
## was `contains'
## 1: 1, 2, 3, 4, 5, 6
如预期的那样,六个点均包含在内。
现在,让我们把数据转成兰勃特等角投影(Lambert conformal projection)。
# 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()
函数来验证:
st_contains(sf1.lam, pt1.lam)## Sparse geometry binary predicate list of length 1, where the predicate
## was `contains'
## 1: 1, 3, 5
为了解决这个问题,我们应当让构成多边形边界的顶点更密集。顶点密度取决于保持这种包含关系所需要的分辨率大小,而最好的方法是试验。
我们使用st_segmentize()
函数设置边界顶点的间隔距离为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, 6
9 「创建天梭指示圆(Tissot indicatrix circles)」
大多数投影都会扭曲一部分空间属性,特别是面积和形状。可视化这种扭曲性质的一个好方法就是创建大地圆(geodesic circles)。
首先,在经纬度坐标系统下定义作为圆心的点图层:
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")
然后,我们使用geosphere
工具包基于这些点构建大地圆:
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
工具包中的st_area()
,它在经纬度坐标系统下会计算地测面积(geodesic 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度为中心的兰勃特等面积投影(Lambert azimuthal equal area projection):
# 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")
北美大陆的形状和面积都发生了肉眼可见的扭曲。
R语言中的地理/投影坐标系统(下)[翻译]相关推荐
- R语言中的地理/投影坐标系统(上)[翻译]
原文链接:https://mgimond.github.io/Spatial/coordinate-systems-in-r.html. 译文分上.下两篇,这里为上篇. 1 「关于PROJ环境更新的说 ...
- 技巧 | 在R语言中使用高德地图的API进行地理/逆地理编码(地址与经纬度的相互转换)...
高德地图和百度地图都提供了坐标拾取系统,通过坐标查询或坐标反查操作可以查询一个地址对应的经纬度或经纬度对应的地址名称.但是,手动查询的方式效率很低,也不能进行批量查询. 本篇就来介绍在R语言中调用高德 ...
- r语言 rgl 强制过程中_一个R语言中操纵矢量空间数据的标准化工具—sf
注: 本文是R语言sf包的核心开发者和维护者--来自德国明斯特大学的地理信息学教授:Edzer Pebesma 的一篇关于sf包的简介,发表于2018年7月的R语言期刊,主要讲述了sf的定位.功能. ...
- raster | R语言中的空间栅格对象及其基本处理方法(Ⅰ)
前面的系列推文已经完成了对R语言中的两个管理空间矢量数据的工具包(sf和sp)的介绍,以及空间自相关.空间插值等空间分析方法. 这里小编再推出一个系列来介绍R语言中管理空间栅格数据的工具包:raste ...
- R语言中如何进行PCA分析?利用ggplot和prcomp绘制基因表达量分析图
学习笔记的主要内容是在R语言中利用ggplot2进行PCA分析和绘图,包括简单分析与操作流程,对比不同方式得到的结果差异,提供脚本代码供练习. PCA分析的原理 在处理基因差异表达数据时,有时候需要分 ...
- R语言GD包地理探测器分析时报错、得不到结果等情况的解决方案
本文介绍在利用R语言的GD包,实现自变量最优离散化方法选取与执行.地理探测器(Geodetector)操作时,出现各类报错信息.长时间得不到结果等情况的解决方案. 在之前的文章R语言GD包基于 ...
- R开发(part8)--应用R语言中的函数环境空间
学习笔记,仅供参考,有错必纠 文章目录 R开发 应用R语言中的函数环境空间 R语言的函数环境空间 封闭环境 绑定环境 运行环境 调用环境 函数环境空间图示 R开发 应用R语言中的函数环境空间 R语言的 ...
- r语言中正定矩阵由于误差不正定_R语言之数据处理(一)
在上一篇小文中,提到了关于R语言导入数据的一些方法,之后的重点就转向了数据的处理上.数据处理其实在整个数据分析项目中所占用的时间是比较多的,所以根据处理的目的不同,也有不同的处理方法.在R语言中,我通 ...
- r语言中的或怎么表示什么不同_R经典入门 之 R语言的基本原理与概念 -- 200430
一.基本原理 R是一种解释型语言,输入的命令可以直接被执行,不同于C等编译语言需要构成完整的程序才能运行. R的语法非常简单和直观.合法的R函数总是带有圆括号的形式,即使括号内没有内容(如,ls()) ...
最新文章
- 关于js中的时间——计算时间差等
- boost::heap模块实现可变堆的测试程序
- [C++11]initializer_lisr模板类的使用
- ‘MIX_INIT_MP3’ was not declared in this scope,这是什么情况?
- SUSE Linux Enterprise Server 12 SP5 Install
- Centos重置root密码
- 如何查看虚拟机服务器ftp,如何通过FTP工具查看虚拟空间使用了多少?
- UE4 虚幻引擎,处理PBR材质
- 通信原理、模电——部分英文术语对照表
- 关于 IE 浏览器打开时速度过慢的问题
- 大学英语精读第三版(第五册)复习笔记——文章内容摘要
- 关于locahost:8080一直在等待却不报错
- 获取crumbIssuer
- ios 扫码枪外设 键盘模式_想把 iPad 当笔记本电脑用?可以试试这款外接键盘
- 关于华为云短信接口对接问题
- 面经九2023.2.3上午笔试和群面
- Nginx部署Vue项目动态路由刷新404
- java执行sql列名无效_列名无效!java代码里的SQL语句!数据库里可以得到正确为什么放java里出错了?...
- 软件测试Bug评测 之Serverity(严重程度)、Priority(优先级)
- windows2008r2 sp1 官方杀毒软件Microsoft Security Essentials(mse)_php_sir_新浪博客