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

Question

我想计算从5弧分的全局光栅中的所有像素到赤道的距离。得到的距离应以米为单位，考虑地球的曲率，不应被地图投影扭曲。

小区域（例如各个国家）的应用通常通过选择适合给定区域的EPSG代码来克服地图投影引起的扭曲的威胁。由于我的应用程序覆盖了整个地球，这个技巧不起作用。但是，我不需要在各个方面充分展示地球的地图。我只需要一个不会扭曲每个像素和赤道最近点之间距离的投影，即不会扭曲沿着经线的距离的投影。就我而言，PlateCarré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()是否考虑到地球曲率的任何信息。根据相关帖子，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)))

这将把“经验法则”设置为每度纬度为mean(dist/abs(lats[,2])) = 110.8032 km

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

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