Projection across multiple UTM zones in R

Question

I have a set of GPS points in decimal degrees for 120 samples from 13 sampling sites (data below are the centroids)

data.frame(x=c(122.4978, 122.1819, 113.7757, 123.0335, 113.0614, 116.8117, 113.6997, 118.8820, 116.8047, 114.0167, 115.0560, 118.6114, 113.3973), y=c(-16.90613, -17.89392, -23.10492, -16.57618, -25.63859, -20.56915, -25.76536, -19.02967, -19.79827, -21.91585, -21.62384, -20.28056, -26.03371))

I now need to calculate the shortest distance between each pair of points (in km). This part is clear and I've successfully reprojected points in a different project where all points were in the same UTM zone and then using gdistance::shortestPath with a raster and transition matrix. These points here, however, span 3 UTM zones (49 - 51).

I'm unfortunately not very experienced in these matters so I'd appreciate if someone could point me towards a resource that, ideally, explains how to do a projection over multiple zones on a step by step basis. Or, if this is not the right approach, give me a hint on how to proceed.

How do you define "shortest distance"? As the proverbial crow files? Or by minimizing cost given a cost-surface? — Robert Hijmans, Jun 18 '20 at 15:40

Robert Hijmans · Answer 1 · 2020-06-18T22:27:25.390

From your question I am not sure what is clear and what is not. But I will just show you how to get "the shortest distance between each pair of points (in km)"

Your example data

df <- data.frame(
  x=c(122.4978, 122.1819, 113.7757, 123.0335, 113.0614, 116.8117, 113.6997, 118.8820, 116.8047, 114.0167, 115.0560, 118.6114, 113.3973),
  y=c(-16.90613, -17.89392, -23.10492, -16.57618, -25.63859, -20.56915, -25.76536, -19.02967, -19.79827, -21.91585, -21.62384, -20.28056, -26.03371)
)

Solution

library(raster)
p <- pointDistance(df, lonlat=T)
p[1:3, 1:3] / 1000
#          [,1]     [,2] [,3]
#[1,]    0.0000       NA   NA
#[2,]  114.3593    0.000   NA
#[3,] 1141.3329 1049.202    0

Or the equivalent, but with both triangles of the matrix filled.

library(geosphere)
d <- distm(df)
d[1:3, 1:3] / 1000
#          [,1]      [,2]     [,3]
#[1,]    0.0000  114.3593 1141.333
#[2,]  114.3593    0.0000 1049.202
#[3,] 1141.3329 1049.2017    0.000

Another one, with terra, returning a distance matrix (dist object).

library(terra)
v <- vect(df, crs="+proj=longlat +datum=WGS84")
s <- distance(v)
class(s)
#[1] "dist"
s[1:3]/1000
#[1]  114.35935 1141.33295   67.79552

Projection across multiple UTM zones in R

1 Answers1