如何用R空间包（如SF）将爱尔兰格网参考转换为经纬度？

首先明确：爱尔兰格网参考（针对北爱尔兰）对应的CRS确实是EPSG:29903（Irish Transverse Mercator），你之前的搜索是正确的。可以用R的sf包或Python的geopandas完成批量转换，核心步骤是先将格网字符串解析为东距（Easting）和北距（Northing），再通过空间包的坐标转换功能转成WGS84经纬度（EPSG:4326）。

方案一：R + sf包 #

1. 解析格网参考为东距/北距 #

先写一个解析函数，处理格网字符串的字母区偏移和数字精度：

parse_irish_grid <- function(grid_ref) {
  # 清理空格，统一格式
  grid_ref <- gsub(" ", "", grid_ref)
  letter <- substr(grid_ref, 1, 1)
  nums <- substr(grid_ref, 2, nchar(grid_ref))
  
  # 拆分东距、北距数字部分
  n <- nchar(nums)
  east_nums <- substr(nums, 1, n/2)
  north_nums <- substr(nums, n/2 + 1, n)
  
  # 爱尔兰格网字母区对应的坐标偏移（单位：米）
  letter_offsets <- list(
    A = c(E=0, N=400000), B = c(E=100000, N=400000),
    C = c(E=200000, N=400000), D = c(E=300000, N=400000),
    E = c(E=0, N=300000), F = c(E=100000, N=300000),
    G = c(E=200000, N=300000), H = c(E=300000, N=300000),
    J = c(E=0, N=200000), K = c(E=100000, N=200000),
    L = c(E=200000, N=200000), M = c(E=300000, N=200000),
    N = c(E=0, N=100000), O = c(E=100000, N=100000),
    P = c(E=200000, N=100000), Q = c(E=300000, N=100000),
    R = c(E=0, N=0), S = c(E=100000, N=0),
    T = c(E=200000, N=0), U = c(E=300000, N=0)
  )
  offset <- letter_offsets[[letter]]
  
  # 计算精度：8位数字对应10米，6位对应100米，4位对应1公里
  precision <- 10^(5 - nchar(east_nums))
  
  # 计算最终东距、北距
  easting <- offset["E"] + as.integer(east_nums) * precision
  northing <- offset["N"] + as.integer(north_nums) * precision
  
  return(data.frame(Easting = easting, Northing = northing))
}

2. 批量转换为经纬度 #

假设你的DataFrame中格网参考列名为grid_ref，执行以下代码：

library(sf)
library(dplyr)

# 解析格网，添加东距/北距列
NI <- NI %>% bind_cols(parse_irish_grid(.$grid_ref))

# 转换为sf空间对象，指定原CRS为29903
NI_sf <- st_as_sf(NI, coords = c("Easting", "Northing"), crs = 29903)

# 转换为WGS84经纬度（EPSG:4326）
NI_sf_wgs84 <- st_transform(NI_sf, crs = 4326)

# 提取经纬度到原DataFrame，移除空间几何信息
NI <- NI_sf_wgs84 %>%
  mutate(
    Longitude = st_coordinates(.)[,1],
    Latitude = st_coordinates(.)[,2]
  ) %>%
  st_drop_geometry()

方案二：Python + geopandas #

1. 解析格网参考为东距/北距 #

同样先写解析函数：

import pandas as pd
import geopandas as gpd

def parse_irish_grid(grid_ref):
    grid_ref = grid_ref.replace(" ", "")
    letter = grid_ref[0]
    nums = grid_ref[1:]
    n = len(nums)
    
    # 拆分东距、北距数字
    east_nums = nums[:n//2]
    north_nums = nums[n//2:]
    
    # 字母区坐标偏移
    letter_offsets = {
        'A': (0, 400000), 'B': (100000, 400000),
        'C': (200000, 400000), 'D': (300000, 400000),
        'E': (0, 300000), 'F': (100000, 300000),
        'G': (200000, 300000), 'H': (300000, 300000),
        'J': (0, 200000), 'K': (100000, 200000),
        'L': (200000, 200000), 'M': (300000, 200000),
        'N': (0, 100000), 'O': (100000, 100000),
        'P': (200000, 100000), 'Q': (300000, 100000),
        'R': (0, 0), 'S': (100000, 0),
        'T': (200000, 0), 'U': (300000, 0)
    }
    e_offset, n_offset = letter_offsets[letter]
    
    # 计算精度
    precision = 10 ** (5 - len(east_nums))
    
    # 计算东距、北距
    easting = e_offset + int(east_nums) * precision
    northing = n_offset + int(north_nums) * precision
    
    return easting, northing

2. 批量转换为经纬度 #

# 解析格网，添加东距/北距列
NI[['Easting', 'Northing']] = NI['grid_ref'].apply(lambda x: pd.Series(parse_irish_grid(x)))

# 创建GeoDataFrame，指定原CRS
gdf = gpd.GeoDataFrame(
    NI,
    geometry=gpd.points_from_xy(NI.Easting, NI.Northing),
    crs="EPSG:29903"
)

# 转换为WGS84经纬度
gdf_wgs84 = gdf.to_crs("EPSG:4326")

# 提取经纬度到原DataFrame
NI['Longitude'] = gdf_wgs84.geometry.x
NI['Latitude'] = gdf_wgs84.geometry.y

注意事项 #

• 解析函数支持4位、6位、8位数字的格网参考（对应不同精度）
• 如果数据中有小写字母、多余空格等不规范格式，需先做清洗（比如统一转为大写、去除非必要字符）
• 北爱尔兰的格网参考确实使用EPSG:29903，爱尔兰共和国则常用EPSG:29902，根据你的数据来源确认即可

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