R：以公制格式计算沿子午线的长距离

Question

我想计算分辨率为 5 弧分的全局栅格中所有像素到赤道的距离。结果距离应以米为单位，考虑到地球的曲率，并且不应因地图投影而变形。

小区域（例如个别国家）的应用程序通常通过选择适合给定区域的 EPSG 代码来克服地图投影引起的扭曲的威胁。因为我的应用程序覆盖了整个地球，所以这个技巧不起作用。然而，我不需要一张能充分展示地球各个方面的地图。我只需要一个不扭曲每个像素与赤道最近点之间距离的投影，即不扭曲沿子午线距离的投影。就我而言，Plate Carrée 投影 ("+proj=longlat +datum=WGS84") 满足此标准。鉴于此信息，我尝试了以下代码：

r_res <- 1/12
r <- raster(resolution = c(r_res, r_res), crs = "+proj=longlat +datum=WGS84")
p <- as(r, "SpatialPoints")
equator <- st_sfc(st_point(c(-180,0)), st_point(c(0,0)), st_point(c(180,0))) %>% st_combine() %>% st_cast(., "LINESTRING") %>% st_sf(., crs = 4326) %>% st_transform(crs = "+proj=longlat +datum=WGS84") %>% as(., "Spatial")
d <- gDistance(p, equator, byid = T)
dmin <- apply(d, 2, min)
r[] <- dmin

不幸的是，此示例中计算的距离以度而不是米表示，这是由投影的 longlat 格式引起的。此外，我找不到任何关于 rgeos::gDistance() 是否考虑了行星曲率的信息。根据相关 post，gDistance() 甚至不应应用于 longlat 数据。

或者，我将网格投影到另一个以米为单位产生输出的投影 (Mollweide)：

r_res <- 1/12
r <- raster(resolution = c(r_res, r_res), crs = "+proj=longlat +datum=WGS84")
r <- projectRaster(r, crs="+proj=moll +ellps=WGS84")
p <- as(r, "SpatialPoints")
equator <- st_sfc(st_point(c(-180,0)), st_point(c(0,0)), st_point(c(180,0))) %>% st_combine() %>% st_cast(., "LINESTRING") %>% st_sf(., crs = 4326) %>% st_transform(crs="+proj=moll +ellps=WGS84") %>% as(., "Spatial")
d <- gDistance(p, equator, byid = T)
dmin <- apply(d, 2, min)
r[] <- dmin

然而，在这种情况下，我不确定 gDistance() 是否纠正了 Mollweide 投影隐含的距离失真。

我知道可以使用经验法则计算距赤道的距离：latitude * 111 km。但是，对于需要更高精度函数的应用程序来说，这种近似值太不精确了。

如果有人能提供一些建议就太好了。随意依赖 gDistance() 以外的距离函数。只要距离不失真，以米为单位测量并考虑地球的曲率，我对任何功能都很好，无论是来自 sf、gdistance、raster、geosphere、rgeos 或任何其他包。

Answer 1

这是一种方法（对于非常大的栅格来说内存不安全）

library(raster)
r_res <- 1/12
r <- raster(resolution = c(r_res, r_res), crs = "+proj=longlat +datum=WGS84")
# latitudes for one column
lats <- cbind(0, yFromRow(r, 1:nrow(r)))
#distances (in km) for one column
dist <- pointDistance(lats, cbind(0,0), lonlat=TRUE) / 1000
# assign to all cells
d <- setValues(r, rep(dist, each=ncol(r)))

这会将 "rule of the thumb" 设置为 mean(dist/abs(lats[,2])) = 110.8032 公里每度纬度

这是一个内存安全（但效率低下）的方法：

x <- init(r, "y")
f <- function(i) { pointDistance(cbind(0, i), cbind(0,0), lonlat=TRUE) / 1000} 
z <- calc(x, fun=f)

R：以公制格式计算沿子午线的长距离

R: Compute long distances along meridians in metric format

r

distance

map-projections

r-raster