如何通过gArea函数将UTM坐标系家域计算结果转换为km²单位？

问题描述 #

使用如下UTM投影坐标样例数据开展动物家域计算：

x        y
437850.3 7220701
465101.3 7210903
489314.6 7159065
513795.7 7114472
532871.0 7075753

使用如下R代码完成95%核密度家域估算流程：

# 加载依赖包
library(rgdal)
library(ks)
library(rgeos)
library(maptools)

# 为点数据设置UTM投影
coordinates(data) <- c("x","y")
proj4string(data) <- CRS( "+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs" )
data<-as.data.frame(data)

# 插件法计算核密度家域
h1 <- Hpi(data[,1:2], pilot = "samse", binned = T)
kernPI1 <- kde(data[,1:2], H = h1)
cont = contourLevels(kernPI1, cont = 95)
line = contourLines(x = kernPI1$eval.points[[1]], y =
                  kernPI1$eval.points[[2]], z = kernPI1$estimate,
                level = cont)
sldf = ContourLines2SLDF(line)
sldf = SpatialLines2PolySet(sldf)
sldf = PolySet2SpatialPolygons(sldf)

gArea(sldf) # 运行输出结果为2010204962

执行gArea(sldf)得到面积值2010204962，需要转换为km²单位；根据先验知识研究区域家域实际范围应为1000-10000 km²，初步推测结果量级偏差与UTM坐标单位设置有关，需要明确具体的转换与修正方法。

解答 #

你的计算结果不存在量级偏差，核心问题是单位换算逻辑不清晰，同时代码存在投影属性丢失的隐患：

• 单位转换方法：你设置的UTM投影单位为米（m），gArea()在正确识别投影时返回的面积单位就是平方米（m²）。1 km² = 1000m × 1000m = 10^6 m²，直接做数值除法即可完成单位转换：

  # 平方米转换为平方公里
  area_km2 <- gArea(sldf) / 10^6
  

  代入你的结果计算可得：2010204962 ÷ 1000000 = 2010.2 km²，完全落在1000-10000 km²的预期区间内，不存在坐标单位设置错误导致的量级偏差。
• 代码隐患修正：你在通过PolySet2SpatialPolygons生成最终家域面对象时，没有为输出对象绑定投影参考系，此时gArea()是直接基于坐标数值计算面积，本次结果数值恰好匹配米制单位计算结果，但后续开展投影转换、空间叠加等操作时会触发错误，需要补全投影赋值步骤：

  sldf = PolySet2SpatialPolygons(sldf)
  # 为生成的家域面绑定与原始点数据一致的投影信息
  proj4string(sldf) <- CRS("+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs")
  # 重新计算面积并转换单位
  gArea(sldf) / 10^6
  

注意：禁止直接使用未投影的WGS84经纬度坐标计算面积，该场景下gArea()返回值单位为平方度，和实际面积会存在多个数量级的偏差，必须先将数据转换为米制单位的平面投影再执行面积计算。

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