如何计算Albers投影下经纬度网格与坐标轴的交点坐标？

嘿，我看你现在的问题是要找到Albers投影下绘图边界和经纬度网格的交点坐标，用来给ggplot的坐标轴标注经纬度度数对吧？其实你不用纠结反向求解Y的思路，换个更直接的方式就行——我们可以先找到这些交点在WGS84下的位置，再转换到Albers投影得到对应的坐标，这样就能轻松设置刻度了。下面是具体的实现步骤和代码，结合你现有的代码修改：

你需要的是绘图边界（x轴/ y轴）与目标经纬度网格线的交点，这些交点的经纬度是确定的（比如y轴左边界与北纬40°的交点，本质是北纬40°纬线和绘图左边界的交叉点）。我们可以先构造这些经纬线的WGS84坐标，转换到Albers投影后，找到它们和绘图边界的交点，就能得到需要的刻度坐标。

library(raster)
library(sp)
library(ggplot2)
library(plyr)

# 你的投影定义（修正原代码中笔误的aea_China为aea）
aea <- CRS('+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +x_0=4000000 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 ')
wgs84 <- CRS('+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0')

# 绘图范围
ext <- c(3640676 ,5672676 ,3683208 ,5619208 )

# --- 生成y轴刻度：左边界与目标纬度的交点 ---
# 自定义需要标注的纬度（可根据需求调整）
target_lats <- seq(35, 55, 5)
lat_ticks <- lapply(target_lats, function(lat) {
  # 创建横跨中国区域的北纬lat纬线（WGS84）
  lon_seq <- seq(70, 140, 0.1)
  line_df <- data.frame(lon = lon_seq, lat = lat)
  line_sp <- SpatialLines(list(Lines(list(Line(line_df)), ID = "1")), proj4string = wgs84)
  
  # 转换到Albers投影
  line_aea <- spTransform(line_sp, aea)
  line_aea_df <- as.data.frame(line_aea)
  
  # 找到最接近绘图左边界（x=ext[1]）的点，即为交点
  closest_idx <- which.min(abs(line_aea_df$x - ext[1]))
  data.frame(aea_y = line_aea_df$y[closest_idx], lat_label = paste0(lat, "°N"))
}) %>% do.call(rbind, .)

# --- 生成x轴刻度：下边界与目标经度的交点 ---
# 自定义需要标注的经度（可根据需求调整）
target_lons <- seq(90, 120, 5)
lon_ticks <- lapply(target_lons, function(lon) {
  # 创建横跨中国区域的东经lon经线（WGS84）
  lat_seq <- seq(10, 60, 0.1)
  line_df <- data.frame(lon = lon, lat = lat_seq)
  line_sp <- SpatialLines(list(Lines(list(Line(line_df)), ID = "1")), proj4string = wgs84)
  
  # 转换到Albers投影
  line_aea <- spTransform(line_sp, aea)
  line_aea_df <- as.data.frame(line_aea)
  
  # 找到最接近绘图下边界（y=ext[3]）的点，即为交点
  closest_idx <- which.min(abs(line_aea_df$y - ext[3]))
  data.frame(aea_x = line_aea_df$x[closest_idx], lon_label = paste0(lon, "°E"))
}) %>% do.call(rbind, .)

# --- 绘制经纬度网格（沿用你的代码并修正）---
lon_tick <- seq(-180,180,5)
lat_tick <- seq(-90,90,5)
j_line <- plyr::alply(lon_tick,1,function(x) cbind(x,lat_tick))
i_line <- plyr::alply(lat_tick,1,function(x) cbind(lon_tick,x) )
lines <- spLines(c(i_line,j_line),crs=wgs84)
lonlat_line <- spTransform(lines, aea)
lonlat_shp <- SpatialLinesDataFrame(lonlat_line, data = data.frame(ID=1))

# --- 最终绘图，替换刻度为经纬度 ---
p <- ggplot() +
  geom_polygon(data=lonlat_shp,aes(x=long,y=lat,group=group), fill=NA,color = "black",linetype=2) +
  coord_cartesian(xlim=ext[1:2],ylim=ext[3:4]) +
  # x轴替换为经度刻度
  scale_x_continuous(breaks = lon_ticks$aea_x, labels = lon_ticks$lon_label) +
  # y轴替换为纬度刻度
  scale_y_continuous(breaks = lat_ticks$aea_y, labels = lat_ticks$lat_label) +
  labs(x = "经度", y = "纬度") +
  theme_bw()

print(p)

1. 为什么不用反向求解？因为投影转换是双向的，但直接构造经纬线找交点的方式更直观，也避免了复杂的逆运算；
2. 你可以根据自己的绘图范围，灵活调整target_lats、target_lons，以及构造经纬线时的lon_seq、lat_seq范围，确保覆盖你的绘图区域；
3. 原代码中的aea_China是笔误，我统一替换成了你定义的aea投影。

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