在R中如何将经纬度转换为LCC投影并提取NC文件时序数据

步骤1：加载依赖包 #

确保安装并加载raster和sf包（sf处理空间转换更直观）：

install.packages(c("raster", "sf"))
library(raster)
library(sf)

步骤2：创建目标点并转换投影 #

核心是严格复用RasterBrick的自定义LCC投影参数——你的LCC是自定义配置，没有对应标准EPSG代码，手动输入参数极易出错：

# 定义目标经纬度点（注意顺序：经度在前，纬度在后，WGS84坐标系）
target_wgs84 <- st_sfc(st_point(c(-158.844437, 18.241755)), crs = 4326)

# 直接提取RasterBrick的投影信息，避免手动编写出错
lcc_proj <- crs(brick.nc)

# 将WGS84点转换为LCC投影
target_lcc <- st_transform(target_wgs84, lcc_proj)

步骤3：提取时间序列 #

用raster包的extract函数直接提取，无需加载全量数据到内存：

# 提取该点的时间序列，可选bilinear插值（适合气温、降水等连续变量）
ts_result <- extract(brick.nc, target_lcc, method = "bilinear")

# 转换为标准时间序列格式（前提是栅格的时间元数据已正确读取）
ts_final <- ts(t(ts_result), start = start(brick.nc), frequency = 365)

排查错误的关键点 #

• 经纬度顺序：多数GIS工具采用经度（X）在前，纬度（Y）在后，之前转换错误大概率是顺序搞反了。
• 范围校验：提前确认转换后的点是否在栅格覆盖范围内，避免提取结果为NA：

# 检查LCC坐标是否在栅格范围内
point_coords <- st_coordinates(target_lcc)
ext_brick <- extent(brick.nc)
if(all(point_coords >= c(ext_brick@xmin, ext_brick@ymin)) && all(point_coords <= c(ext_brick@xmax, ext_brick@ymax))) {
  cat("点在栅格范围内，可正常提取\n")
} else {
  cat("点不在栅格覆盖区域，提取结果为NA\n")
}

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