在R中转换Behrmann投影至WGS84经纬度投影遇错求助

我之前处理全球网格shapefile时也碰到过一模一样的问题！这种非有限点错误大多是因为个别多边形的顶点在投影转换时超出了有效范围，或者原始几何本身有小问题（比如自相交、极点附近的异常点）。咱们一步步来搞定它：

先定位问题根源 #

首先，咱们找出那个出问题的Polygon 1106，看看它的情况：

# 加载必要包
library(raster)
library(rgdal)
library(maptools)

data("wrld_simpl")
tmp <- tempfile()
download.file("https://github.com/darunabas/extras/blob/master/temp_shapefile.zip?raw=true", destfile = tmp)
unzip(tmp, exdir = ".")
s <- rgdal::readOGR("temp_shapefile")
proj4string(s) = CRS("+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +ellps=WGS84 +units=m +no_defs")

# 定位并查看问题要素
problem_poly <- s[1106, ]
plot(problem_poly)

你会发现这个多边形大概率在南极/北极附近——Behrmann投影对极点附近的点转换时容易出现数值溢出，导致非有限值。

方案1：用SF包（推荐，更现代稳定） #

现在R的空间数据处理已经逐渐转向sf包了，它处理几何错误和投影转换的能力比旧的sp/rgdal更健壮。试试这个：

rm(list = ls())
library(sf)
library(maptools)
data("wrld_simpl")

# 下载解压数据
tmp <- tempfile()
download.file("https://github.com/darunabas/extras/blob/master/temp_shapefile.zip?raw=true", destfile = tmp)
unzip(tmp, exdir = ".")

# 读取并设置原始投影
s_sf <- st_read("temp_shapefile")
st_crs(s_sf) <- "+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +ellps=WGS84 +units=m +no_defs"

# 先修复几何，再转换投影
s_sf_longlat <- s_sf %>%
  st_make_valid() %>%  # 自动修复自相交、无效几何等问题
  st_transform(crs = "+proj=longlat +datum=WGS84")

# 如果还有个别无效要素，过滤掉
valid_check <- st_is_valid(s_sf_longlat)
s_sf_clean <- s_sf_longlat[valid_check, ]

# 验证和wrld_simpl匹配
plot(st_geometry(wrld_simpl))
plot(st_geometry(s_sf_clean), add = TRUE, col = "orange", lwd = 0.5)

st_make_valid会帮你自动修复大多数几何小问题，st_transform处理边缘点时也会更智能，基本能解决这个错误。

方案2：继续用SP/RGDAL包 #

如果你习惯用旧的sp体系，也可以先修复几何再转换：

# 用缓冲0的方法修复几何错误
s_fixed <- gBuffer(s, byid = TRUE, width = 0)
# 再尝试转换
p <- spTransform(s_fixed, CRS("+proj=longlat +datum=WGS84"))

# 如果还是不行，直接移除出问题的要素
s_clean <- s[-1106, ]
p <- spTransform(s_clean, CRS("+proj=longlat +datum=WGS84"))

缓冲0的操作会重新生成几何，修复自相交等问题；如果还是不行，直接删掉那个有问题的多边形，不影响整体网格的使用。

总结一下 #

• 优先用sf包，它是现在R空间处理的主流，坑更少；
• 转换前一定要确保几何有效，这是解决这类投影错误的关键；
• 极点附近的网格要素是这类问题的高发区，针对性处理即可。

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