R语言地图图层不匹配求助：土耳其断层线地图绘制

关键问题修正 #

1. 停止错误覆盖坐标系 #

读取GEM断层数据时，不要直接用st_set_crs强制设定CRS——如果原始shp文件自带正确的坐标信息，这会强制覆盖导致坐标错位。正确流程是：

# 先读取原始数据并查看其CRS
shapename <- read_sf('gem_active_faults_harmonized.shp')
print(st_crs(shapename)) # 确认原始坐标系

# 仅当原始CRS不是WGS84（EPSG:4326）时，执行转换
if (st_crs(shapename)$epsg != 4326) {
  shapename <- shapename %>% st_transform(crs = "EPSG:4326")
}

2. 移除冗余坐标转换 #

gisco_get_countries已经指定epsg="4326"，和你定义的crsLONGLAT完全一致，无需再次转换。简化获取土耳其边界的函数：

get_sf <- function() {
  giscoR::gisco_get_countries(
    year = "2016",
    epsg = "4326",
    resolution = "10",
    country = "Turkey"
  )
}
country_transformed <- get_sf()

3. 裁剪断层数据到土耳其范围 #

将断层数据裁剪至土耳其边界内，既减少无关数据干扰，也能确保图层精准对齐：

# 裁剪断层数据
shapename_cropped <- st_crop(shapename, country_transformed)

修正后的完整绘图代码 #

library(elevatr)
library(terra)
library(tidyverse)
library(sf)
library(giscoR)
library(marmap)

# 获取土耳其地理数据
get_sf <- function() {
  giscoR::gisco_get_countries(
    year = "2016",
    epsg = "4326",
    resolution = "10",
    country = "Turkey"
  )
}
country_transformed <- get_sf()

# 读取并处理断层线数据
shapename <- read_sf('gem_active_faults_harmonized.shp')
# 转换坐标系（仅当原始CRS不符时）
if (st_crs(shapename)$epsg != 4326) {
  shapename <- shapename %>% st_transform(crs = "EPSG:4326")
}
# 裁剪到土耳其范围
shapename_cropped <- st_crop(shapename, country_transformed)

# 绘制地图
ggplot() +
  geom_sf(data = country_transformed) +
  geom_sf(data = shapename_cropped, color = "orange", size = 1) +
  coord_sf(xlim = c(26, 45), ylim = c(36, 42)) +
  theme(legend.position = "none")

额外验证提示 #

• 确认下载的GEM断层数据文件完整，需包含.shp、.shx、.dbf等配套文件，缺失会导致读取异常。
• 检查coord_sf的经纬度范围是否符合需求，当前设置的东经26°-45°、北纬36°-42°覆盖了土耳其主要区域。

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