引言
过去一次推送中，记录了怎么在地图上标注一个长方形采样区域。
（点击 👉 即可跳转。）
1 需要的包
需要从 GitHub 安装
rnaturalearthhires。
remotes::install_github("ropensci/rnaturalearthhires")加载包。
library(tidyverse)
library(sf)
library(rnaturalearth)
library(rnaturalearthhires)2 获取北美洲地图
Map1 <- ne_countries(
continent = "north america",
scale = "small", # 低分辨率
returnclass = "sf"
)3 制作地图轮廓
这段代码定义了一个名为 MakeOutline 的函数，用于生成一个给定区域的轮廓坐标。
函数接收经纬度的最小和最大值作为参数，以及生成轮廓的点数 N。函数首先创建经度 Longs 和纬度 Lats 序列，然后合并这两个序列以形成轮廓的坐标。
最后，通过调用 plot 函数，我们可以将生成的轮廓坐标以图形的形式展示出来。这里使用 5 个点作为示例进行测试。
# Map outline coordinates function
MakeOutline <- function(LonMin = 0, LonMax = 90, LatMin = 0, LatMax = 90, N = 10) {
# Longitudes
Longs <- c(
rep(LonMin, N), # 左边
seq(LonMin, LonMax, len = N), # 上边
rep(LonMax, N), # 右边
seq(LonMax, LonMin, len = N), # 下边
LonMin # 闭合
)
# Latitudes
Lats <- c(
seq(LatMin, LatMax, len = N), # 左边
rep(LatMax, N), # 上边
seq(LatMax, LatMin, len = N), # 右边
rep(LatMin, N), # 下边
LatMin # 闭合
)
# Return values
Outline <- cbind(Longs, Lats)
return(Outline)
}
# 测试函数
plot(MakeOutline(N = 5), pch = 19, axes = FALSE, asp = 1)4 设定地图范围
指定了一个具体的地理区域，其经度范围从 -120 到 -60，纬度范围从 20 到 70。
函数调用时设置了 N = 10，这意味着每条边将有 10 个点，从而生成一个相对精细的轮廓。
使用 head(10) 函数，我们仅展示了生成的轮廓坐标中的前 10 个点，作为预览。
LonMin <- -120
LonMax <- -60
LatMin <- 20
LatMax <- 70
# Check out
MakeOutline(
LonMin = LonMin,
LonMax = LonMax,
LatMin = LatMin,
LatMax = LatMax,
N = 10 # 轮廓每条边有多少个点
) |>
head(10) Longs Lats
[1,] -120 20.00000
[2,] -120 25.55556
[3,] -120 31.11111
[4,] -120 36.66667
[5,] -120 42.22222
[6,] -120 47.77778
[7,] -120 53.33333
[8,] -120 58.88889
[9,] -120 64.44444
[10,] -120 70.000005 构造多边形
使用 MakeOutline 函数创建一个多边形，并将其投影到 WGS84 坐标系中。
首先，我们定义了一个 CRS（Coordinate Reference System，坐标参考系统）Proj1，它代表 WGS84 经纬度坐标系。
然后，把轮廓数据封装在一个列表中，以便使用 sf 包中的 st_polygon 函数转换为一个地理多边形。
最后，我们使用 st_sfc 函数创建一个简单特征集合，指定其 CRS 为 Proj1。
Proj1 <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
Outline1 <- list(
MakeOutline(
LonMin = LonMin,
LonMax = LonMax,
LatMin = LatMin,
LatMax = LatMax,
N = 10
)
) |>
st_polygon() |>
st_sfc(crs = Proj1)7 最终地图
首先，将 Outline1 从 WGS84 坐标系转换到另一个坐标系 Proj2，并计算转换后多边形的边界框 Lim1。
接着，使用 ggplot2 包创建一个地图，其中包含了三个图层：底层是一个灰色的地图 Map1，中间是之前创建的橙色轮廓 Outline1，最上层是 60 个随机分布的点。这些点在指定的经纬度范围内生成，并转换为简单的地理特征。
最后，通过 coord_sf 函数设置坐标系统为 Proj2，并调整x轴和y轴的范围以及比例，使得地图展示更为合适。
clip = "off"确保了超出设定边界的部分不会被裁剪掉。c(1, 1.2)是一个乘数因子，它被应用到ylim的上下界上。这意味着 y 轴的最小值将被设置为Lim1[2]（最小纬度）乘以1，而最大值将被设置为Lim1[4]（最大纬度）乘以1.2。给 y 轴上界的显示提供了额外的空间。
# 新投影
Proj2 <- "+proj=laea +lon_0=-100 +lat_0=30"
# 轮廓应用新投影后的范围
Outline1Trans <- Outline1 |>
st_transform(crs = Proj2)
Lim1 <- st_bbox(Outline1Trans)
print(Lim1) xmin ymin xmax ymax
-2084150 -1106436 4068494 4722623# 设置绘图主题和字体
FONT <- "sans"
theme_set(
theme_grey(
base_family = FONT,
base_size = 16,
base_line_size = 0.5
) +
theme(axis.ticks = element_blank())
)
# 绘图
ggplot() +
geom_sf(
data = Map1,
color = "white",
fill = "grey50",
linewidth = 0,
alpha = 0.5
) +
geom_sf(
data = Outline1,
color = "orange",
fill = NA,
linewidth = 0.7
) +
geom_sf(
data = data.frame(
x = runif(60, LonMin, LonMax),
y = runif(60, LatMin, LatMax)
) |>
st_as_sf(
coords = c("x", "y"),
crs = "+proj=longlat"
)
) +
scale_fill_brewer(guide = "none", palette = "Paired") +
scale_x_continuous(breaks = seq(LonMin, LonMax, len = 3)) +
scale_y_continuous(breaks = seq(LatMin, LatMax, len = 3)) +
coord_sf(
crs = Proj2,
xlim = c(Lim1[1], Lim1[3]),
ylim = c(Lim1[2], Lim1[4]) * c(1, 1.2),
clip = "off" # off 则不裁剪超出 lim 部分
)8 裁剪地图
使用 sf 包来裁剪地图 Map1，使其只显示在之前创建的轮廓 Outline1 范围内。
首先，通过 sf_use_s2(FALSE) 禁用 S2 库的地理计算。
然后，使用管道操作符 |> 将 Map1 传递给 st_intersection 函数，该函数与 Outline1 做交集运算，从而得到新的地图 Map2，它仅包含轮廓内的部分。
这样，我们就能够得到一个与指定轮廓完全匹配的地图区域。
# 地图裁剪到轮廓范围内
sf_use_s2(FALSE)
Map2 <- Map1 |>
st_intersection(Outline1)
# 地图
set.seed(1)
P0 <- ggplot() +
geom_sf(
data = Map2,
color = "white",
fill = "grey50",
linewidth = 0,
alpha = 0.5
) +
geom_sf(
data = Outline1,
color = "grey0",
fill = NA,
linewidth = 0.3,
linetype = 2
) +
geom_sf(
data = data.frame(
x = runif(60, LonMin, LonMax),
y = runif(60, LatMin, LatMax),
g = rep(c("T1", "T2", "Ck"), 20)
) |>
st_as_sf(
coords = c("x", "y"),
crs = "+proj=longlat"
),
aes(fill = g, shape = g),
size = 3
) +
scale_fill_brewer(name = NULL, palette = "Paired") +
scale_shape_manual(name = NULL, values = c(21, 22, 24)) +
scale_x_continuous(breaks = seq(LonMin, LonMax, len = 3)) +
scale_y_continuous(breaks = seq(LatMin, LatMax, len = 3)) +
labs(
x = NULL,
y = NULL
) +
coord_sf(
crs = Proj2,
xlim = c(Lim1[1], Lim1[3]),
ylim = c(Lim1[2], Lim1[4]) * c(1, 1.2),
clip = "off" # off 则不裁剪超出 lim 部分
)
P0ggsave(
"map-add-polygon-0.png",
plot = P0,
width = 6,
height = 4,
dpi = 700,
bg = "white"
)9 直接在地图边缘添加坐标
如题，也可以不要坐标轴，直接在地图边缘添加坐标。
Outline2Df <- MakeOutline(
LonMin = LonMin,
LonMax = LonMax,
LatMin = LatMin,
LatMax = LatMax,
N = 2
) |>
as.data.frame()
print(Outline2Df) Longs Lats
1 -120 20
2 -120 70
3 -120 70
4 -60 70
5 -60 70
6 -60 20
7 -60 20
8 -120 20
9 -120 20Outline2Trans <- Outline2Df |>
st_as_sf(
coords = c("Longs", "Lats"),
crs = Proj1
) |>
st_transform(crs = Proj2) |>
st_coordinates() |>
data.frame(
Longs = paste0(Outline2Df$Longs, "°W"),
Lats = paste0(Outline2Df$Lats, "°N")
)
print(Outline2Trans) X Y Longs Lats
1 -2084150.1 -937575.1 -120°W 20°N
2 -801059.8 4449924.8 -120°W 70°N
3 -801059.8 4449924.8 -120°W 70°N
4 1528271.3 4722623.3 -60°W 70°N
5 1528271.3 4722623.3 -60°W 70°N
6 4068493.8 -427071.9 -60°W 20°N
7 4068493.8 -427071.9 -60°W 20°N
8 -2084150.1 -937575.1 -120°W 20°N
9 -2084150.1 -937575.1 -120°W 20°N# 直接在地图边缘添加坐标
P1 <- P0 +
geom_text(
data = Outline2Trans,
aes(X, Y, label = paste(Longs, Lats, sep = ", ")),
size = 3,
color = "grey0",
family = FONT
) +
theme_void()
P1ggsave(
"map-add-polygon-1.png",
plot = P1,
width = 6,
height = 4,
dpi = 700,
bg = "white"
)
P2 <- P0 +
geom_label(
data = Outline2Trans,
aes(X, Y, label = paste(Longs, Lats, sep = ", ")),
size = 3,
color = "grey0",
label.r = unit(0, "pt"),
family = FONT
) +
theme_void()
P2ggsave(
"map-add-polygon-2.png",
plot = P2,
width = 6,
height = 4,
dpi = 700,
bg = "white"
)