如何选择正确CRS以绘制从其他地图提取的坐标

我用juicr包从一张历史物种分布图中提取了点的坐标，但在R中把这些点绘制到世界地图上时，它们显示在错误的位置。我猜测这是CRS（坐标参考系统）的问题，但不知道该如何确定正确的投影，也不清楚怎么把提取的点投影到正确位置。

以下是我用来绘制坐标的代码：

# 设置主题
theme_set(theme_bw())

# 创建基础地图
world <- ne_countries(scale = "medium", returnclass = "sf")
class(world)

base.map <- ggplot(data = world) + 
  geom_sf(color = "dark gray", fill = "white") +
  coord_sf(xlim = c(20, 62), ylim = c(22, 48), expand = FALSE) +
  theme(panel.background = element_rect(fill = "#c6dbef")) +
  theme(plot.title = element_text(lineheight=.8, face="bold"), legend.position = "none")


# 创建地中海柏木分布图
# 从原始地图提取的坐标样本

csem$x <- c(21.06574, 21.45329, 21.64706, 21.84083, 22.95502, 23.68166, 24.02076, 23.87543, 24.0692) 
csem$y <- c(31.85996, 31.85996, 31.62298, 31.82047, 34.50628, 34.62478, 34.46679, 34.38779, 34.34829)
 


# 绘制地中海柏木种群分布
csem.map <- base.map +
  geom_point(data = csem, aes(x = x, y = y), size = 0.75)+
  scale_color_manual(values=c("#56B4E9", "black")) +
  labs(title = "Browicz记录的地中海柏木分布") + 
  theme(plot.title = element_text(lineheight=.8, face="bold"), 
        plot.caption = element_text(face = "italic"), 
        legend.position = "none") + 
  theme(axis.title  = element_blank()) +
  theme(plot.caption = element_text(hjust = 0))

csem.map

核心问题是：你用juicr提取的是图像的像素坐标，而非地理坐标系（如WGS84）下的经纬度，直接和ne_countries返回的WGS84坐标（EPSG:4326）叠加自然会错位。你需要对提取的像素坐标做地理配准，步骤如下：

1. 获取地理控制点 #

从原始历史地图上找到至少3个已知真实地理坐标的点：

• 比如地图上标注的经纬线交点、知名城市（如开罗、雅典）、海岸线的明显拐点
• 用juicr提取这些点的像素坐标（x,y），同时查找它们对应的真实经纬度（lon, lat，WGS84坐标系）

2. 在R中执行地理配准 #

用sf和gdalUtils包完成坐标转换，示例代码如下：

首先安装并加载所需包：

install.packages(c("sf", "gdalUtils", "ggplot2", "rnaturalearth"))
library(sf)
library(gdalUtils)
library(ggplot2)
library(rnaturalearth)

步骤1：定义控制点 #

假设你找到了3组匹配的控制点：

# 像素坐标（juicr提取的x,y）
pixel_pts <- data.frame(x = c(px1, px2, px3), y = c(py1, py2, py3))
# 对应的真实地理坐标（经纬度，WGS84）
geo_pts <- data.frame(lon = c(lon1, lon2, lon3), lat = c(lat1, lat2, lat3))

步骤2：生成配准后的坐标 #

将控制点用于推导变换规则，转换所有提取的点：

# 将控制点转为sf对象
pixel_sf <- st_as_sf(pixel_pts, coords = c("x", "y"), crs = NA)
geo_sf <- st_as_sf(geo_pts, coords = c("lon", "lat"), crs = 4326)

# 计算仿射变换参数并转换所有提取点
csem_sf <- st_as_sf(csem, coords = c("x", "y"), crs = NA)
csem_geo <- st_transform(csem_sf, crs = st_crs(geo_sf), affine = st_affine(pixel_sf, geo_sf))

步骤3：加载配准后的点并绘图 #

# 重新绘图
csem.map <- base.map +
  geom_sf(data = csem_geo, size = 0.75, color = "#56B4E9") +
  labs(title = "Browicz记录的地中海柏木分布") + 
  theme(plot.title = element_text(lineheight=.8, face="bold"), 
        plot.caption = element_text(face = "italic"), 
        legend.position = "none") + 
  theme(axis.title  = element_blank()) +
  theme(plot.caption = element_text(hjust = 0))

csem.map

3. 确定原始地图的CRS #

如果原始历史地图有标注投影信息，直接使用对应的EPSG代码；如果没有，通过以下方式推导：

• 查看地图的出版信息，通常会标注投影类型（如UTM、兰伯特投影等）
• 用多个控制点拟合投影参数，通过sf::st_crs验证转换后的点是否匹配地图上的位置

内容的提问来源于stack exchange，提问作者Harry