如何将sf类多边形t1投影为UTM并转换为owin对象？

核心逻辑 #

spatstat的as.owin()仅支持平面投影坐标（如UTM），不支持经纬度类型的地理坐标，因此需要先将ETRS89经纬度的sf对象转换至对应区域的UTM投影，再执行转换操作。

具体步骤与代码示例 #

• 步骤1：计算多边形所在的UTM带号
  ETRS89与WGS84经纬度坐标兼容度极高，可通过多边形中心经度计算UTM带号：

  # 获取多边形中心坐标
  centroid <- st_centroid(t1)
  # 提取中心经度（单位：度）
  lon <- st_coordinates(centroid)[, "X"]
  # 计算UTM带号
  utm_zone <- floor((lon + 180) / 6) + 1
  

• 步骤2：构建UTM投影的EPSG代码
  根据所在半球选择对应EPSG编码前缀：北半球用326XX，南半球用327XX（XX为计算出的带号），示例为北半球情况：

  utm_crs <- paste0("EPSG:326", utm_zone)
  

• 步骤3：将sf对象转换至UTM投影
  借助sf包的st_transform()完成投影转换：

  t1_utm <- st_transform(t1, crs = utm_crs)
  

• 步骤4：转换为owin对象
  此时可正常使用as.owin()完成转换：

  library(spatstat.geom)
  t1_owin <- as.owin(t1_utm)
  

额外说明 #

• 若不确定半球，可通过中心纬度判断：纬度>0为北半球，<0为南半球，对应调整EPSG编码前缀（326/327）。
• 若多边形跨多个UTM带，可选择覆盖范围最广的带号，或考虑使用LAEA等其他平面投影，UTM仍是区域平面投影的首选方案。

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