如何基于经纬度值确认或转换数据框的空间投影？

我正在开展美国地图绘制项目，原始数据为带有Lambert Conformal Conic投影的ncdf文件。将其转换为栅格时经纬度变量出现错误，因此我将各变量单独处理后存入数据框。

将数据绘图后，在采用"+proj=longlat +datum=NAD83 +no_defs"投影的美国平面shapefile上显示为弯曲状态。我尝试将数据框中的经纬度转换为坐标点并重新投影至shapefile的投影，但未取得效果。

相关文件：

• 数据框：Data1deg.xlsx
• Shapefile：cb_2018_us_state_20m.shp

我的尝试代码：

library(sf)
library(sp)
library(rgdal)
library(readxl)

#Processed data
datafile <- readxl::read_xlsx("Data1deg.xlsx")

#Shape File
Shape <- readOGR("cb_2018_us_state_20m.shp")
Shapecrs <- Shape@proj4string@projargs
statenums <- c("0", "1", "2", "3", "4", "5", "6", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35", "37", "38", "39", "40", "41", "42", "43", "44", "45", "46", "47", "49", "50", "51")
Shape1 <- Shape[statenums, ]

#conrvrting dataframe long and lat to coordinate points
coordinates <-  st_as_sf(datafile, coords = c("xv", "yv"))
st_is_longlat(coordinates)

#original data projection
coordinates_geo <- st_set_crs(coordinates, "+proj=lcc +lat_0=40.0000076293945 +lon_0=-97 +lat_1=30 +lat_2=45 +x_0=0 +y_0=0 +R=6370000 +units=m +no_defs")
plot(coordinates_geo)

#Attempts to reproject data
reprojcoord <- st_transform(coordinates_geo, Shapecrs)
plot(reprojcoord)
reprojcoord1 <- st_transform_proj(coordinates_geo, Shapecrs)
plot(reprojcoord1)

请问是否可以仅通过经纬度值来确认或更改数据框的投影？

核心结论 #

不能仅通过经纬度值自动确认或更改投影，但可以结合已知的原始投影信息手动修正，这也是解决你当前问题的关键。

问题根源分析 #

你当前的错误在于：将数据框中的xv/yv直接当作经纬度处理，但实际上这两个值是Lambert Conformal Conic投影下的平面坐标（米为单位），不是WGS84或NAD83的经纬度。错误地给平面坐标赋予LCC投影后再转地理坐标，导致坐标混乱，最终绘图弯曲。

修正步骤 #

1. 正确定义原始数据的投影
   原始ncdf是LCC投影，xv/yv是该投影下的平面坐标，需先将数据框转为sf对象并赋予正确的LCC投影：

   # 正确创建sf对象：xv/yv是LCC投影的平面坐标，不是经纬度
   data_sf <- st_as_sf(datafile, coords = c("xv", "yv"), crs = "+proj=lcc +lat_0=40.0000076293945 +lon_0=-97 +lat_1=30 +lat_2=45 +x_0=0 +y_0=0 +R=6370000 +units=m +no_defs")
   

2. 转换到shapefile的投影
   直接将上述sf对象转换为shapefile的NAD83地理坐标：

   # 转换到shapefile的投影
   data_nad83 <- st_transform(data_sf, "+proj=longlat +datum=NAD83 +no_defs")
   

3. 合并绘图验证
   用sf包的st_read替代老旧的readOGR读取shapefile，然后一起绘图：

   # 读取并筛选shapefile
   shape_sf <- st_read("cb_2018_us_state_20m.shp") %>%
     filter(STATEFP %in% statenums)
   
   # 绘图验证
   plot(shape_sf$geometry, col = "lightgray")
   plot(data_nad83$geometry, add = TRUE, pch = 16, col = "red")
   

关于投影确认的补充 #

如果完全不知道原始投影，仅靠坐标值无法准确判断（不同投影可能有相似的坐标范围）。但你已知原始ncdf的投影是LCC，直接用该投影参数定义即可。后续遇到类似问题，可通过ncdf4包读取ncdf文件的元数据提取投影信息，避免手动输入出错。

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