I would like to compute distances from all pixels in a global raster of a 5 arc minute resolution to the equator. The resulting distances should be in meters, take the curvature of Earth into account and should not be distorted by map projections.
Applications on small areas (e.g. individual countries) usually overcome the threat of map projection-induced distortions by selecting an EPSG code that fits the given region well. Since my application covers the entire globe, this trick does not work. However, I do not need a map that displays the planet adequately in every aspect. I only need a projection that does not distort the distance between each pixel and the nearest point at the equator, i.e. a projection that does not distort distances along meridians. As far as I am concerned the Plate Carrée projection ("+proj=longlat +datum=WGS84") fulfills this criterion. Given this information I tried the following code:
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
The distances computed in this example are unfortunately expressed in degrees rather than meters, which is caused by the longlat format of the projection. Furthermore, I could not find any information on whether rgeos::gDistance() takes the curvature of the planet into account. And according to a related post, gDistance() should not even be applied to longlat data.
Alternatively, I projected the grid into another projection (Mollweide) that produces output in meters:
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
In this case I am, however, not sure whether gDistance() corrects for the distance distortions implied by the Mollweide projection.
I am aware that distances from the equator can be calculated using a rule of thumb: latitude * 111 km. This approximation is, however, too imprecise for my applications that demand a function of higher precision.
It would be great if anyone could offer some advice. Feel free to rely on distance functions other than gDistance(). As long as distances are not distorted, are measured in meters and take the curvature of the earth into account I am fine with any function, be it from sf, gdistance, raster, geosphere, rgeos or any other package.
Here is an approach (that is not memory safe for very large rasters)
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)))
This would set the "rule of the thumb" to mean(dist/abs(lats[,2])) = 110.8032 km per degree latitude
Here is a memory-safe (but inefficient) method:
x <- init(r, "y")
f <- function(i) { pointDistance(cbind(0, i), cbind(0,0), lonlat=TRUE) / 1000}
z <- calc(x, fun=f)