当前位置：首页 > article >正文

【R 自定义函数总结】投影转换等

article 2025/1/7 20:02:49

函数1：投影坐标系转换

changeCoord 是一个用于 将地理坐标（经纬度） 从一个投影坐标系转换为另一个投影坐标系的 R 函数。它基于 R 的 sp 包，通过投影转换实现地理坐标的变换。

函数功能概述：

输入：一组经纬度坐标（v_lon 和 v_lat），以及输入和输出的投影坐标系（proj_in 和 proj_out）。
检查：首先验证输入经纬度的范围是否合理（经度在 -180 到 180，纬度在 -90 到 90）。
转换：将输入的经纬度坐标按照指定的输入投影（proj_in）进行定义，并利用 spTransform 将其转换为目标投影（proj_out）。
输出：返回经过投影转换后的坐标。

# 函数说明：将地理坐标（经纬度）从一个投影坐标系转换为另一个投影坐标系
changeCoord <- function(v_lon, v_lat, proj_in, proj_out) {
  # 检查输入数据是否合理
  if (any(v_lon < -180 | v_lon > 180) || any(v_lat < -90 | v_lat > 90)) {
    stop("经纬度坐标超出有效范围")
  }
  
  dat0 = data.frame(lon=v_lon, lat=v_lat)
  coordinates(dat0) = c("lon", "lat")
  proj4string(dat0) = proj_in
  dat1 = spTransform(dat0, proj_out)
  return(dat1)
}

函数的代码逐行解析

changeCoord <- function(v_lon, v_lat, proj_in, proj_out) {

定义函数：函数名为 changeCoord，接受以下四个参数：

v_lon：一组经度数据（向量）。
v_lat：一组纬度数据（向量）。
proj_in：输入坐标系的投影字符串（通常为 PROJ.4 格式）。
proj_out：目标坐标系的投影字符串（也是 PROJ.4 格式）。

if (any(v_lon < -180 | v_lon > 180) || any(v_lat < -90 | v_lat > 90)) {
    stop("经纬度坐标超出有效范围")
  }

输入数据检查：

if：检查是否有经度超出 [-180, 180] 或纬度超出 [-90, 90] 的数据。
any：检测是否有任何一个元素不满足条件。
如果超出范围，使用 stop 函数终止程序，并报错 “经纬度坐标超出有效范围”。

dat0 = data.frame(lon=v_lon, lat=v_lat)

创建数据框：
data.frame：创建一个包含经纬度坐标的数据框，列名为 lon（经度）和 lat（纬度）。
例如，输入 v_lon = c(120, 121) 和 v_lat = c(30, 31) 时，dat0 的结果为：

lon lat
1 120  30
2 121  31

coordinates(dat0) = c("lon", "lat")

定义为空间点数据：

coordinates：将数据框 dat0 转换为空间点对象（SpatialPointsDataFrame）。
这里指定 lon 和 lat 列为地理坐标。
转换后，dat0 成为空间点对象，具备空间坐标属性。

proj4string(dat0) = proj_in

指定输入投影坐标系：
proj4string：为 dat0 指定输入的投影信息。
proj_in 是输入坐标系的 PROJ.4 字符串，例如：

WGS84 坐标系：“+proj=longlat +datum=WGS84 +no_defs”。
墨卡托投影：“+proj=merc +datum=WGS84 +units=m +no_defs”。

这一步定义了 dat0 的坐标系，但尚未进行坐标转换。

dat1 = spTransform(dat0, proj_out)

投影变换：
spTransform：将 dat0 的坐标从输入投影（proj_in）转换为目标投影（proj_out）。
proj_out 是目标坐标系的 PROJ.4 字符串，例如：

WGS84 坐标系：“+proj=longlat +datum=WGS84 +no_defs”。
UTM 坐标系（如 50N 带）：“+proj=utm +zone=50 +datum=WGS84 +units=m +no_defs”。

结果 dat1 是一个新的空间点对象，投影已经被转换。

return(dat1)
}

返回结果：
函数返回经过投影转换后的 dat1，它是一个空间点对象，包含转换后的坐标。

示例

示例 1：将 WGS84 地理坐标转换为 UTM 坐标

library(sp)

# 经纬度坐标
lon <- c(120.0, 121.5)
lat <- c(30.5, 31.0)

# 输入和输出坐标系
proj_in <- "+proj=longlat +datum=WGS84 +no_defs"  # WGS84 地理坐标
proj_out <- "+proj=utm +zone=50 +datum=WGS84 +units=m +no_defs"  # UTM 50N 坐标系

# 调用函数
result <- changeCoord(lon, lat, proj_in, proj_out)

# 查看结果
print(result)