I have two data frames, DF1 and DF2, containing spatial data with longitude/latitude at the scale of North America, and their projection is EPSG:4326. I would like to calculate the difference between the mean value (M1) of the y coordinates of DF1 and the mean value (M2) of the y coordinates of DF2, and I want the difference to be expressed in km. I am not sure which projection system to use to convert longitude and latitude to km in my case.
I discovered the function terra::expanse() to compute the area covered by raster cells that are not NA by using the longitude/latitude CRS. It is stated that it is better to use the longitude/latitude CRS to calculate areas.
I could not find the source code for this function to see if there might be a way to convert the values M1 and M2 to km. I am a novice in GIS, so thank you for your help. Here’s a summary of my data frames and the values M1 and M2:
> summary(DF1 [,c(“x”, “y”)])
decimalLongitude decimalLatitude
Min. :-171.62 Min. :14.88
1st Qu.:-118.62 1st Qu.:41.88
Median :-103.38 Median :53.38
Mean :-104.81 Mean :51.95
3rd Qu.: -87.88 3rd Qu.:63.62
Max. : -52.88 Max. :76.38
> summary(DF2 [,c(“x”, “y”)])
decimalLongitude decimalLatitude
Min. :-171.62 Min. :15.62
1st Qu.:-118.62 1st Qu.:42.12
Median :-103.12 Median :53.38
Mean :-104.70 Mean :52.04
3rd Qu.: -87.88 3rd Qu.:63.62
Max. : -52.88 Max. :76.38
> M1 <- mean(DF1[,c(“y”)], na.rm = TRUE)
> M1
[1] 51.95062
> M2 <- mean(DF2[,c(“y”)], na.rm = TRUE)
> M2
[1] 52.03755
> M2 - M1
[1] 0.08692784 ## how to convert to km ?
You can compute the distance between longitude/latitude points with terra::distance or geosphere::distm
Example data
d1 <- data.frame(lon=0:1, lat=0:1)
d2 <- data.frame(lon=30:33, lat=20:23)
To get the actual distances between the points, you could do
x1 <- terra::distance(d1, d2, lonlat=TRUE, unit="km")
mean(x1)
#[1] 4088.824
x2 <- geosphere::distm(d1, d2, distGeo)
mean(x2) / 1000
#[1] 4088.824
Since you only care about the mean latitudinal distance, we can set longitude to zero and use the mean latitudes.
md1 <- data.frame(lon=0, lat=mean(d1$lat))
md2 <- data.frame(lon=0, lat=mean(d2$lat))
terra::distance(md1, md2, lonlat=TRUE, unit="km")
#[1] 2323.15
geosphere::distGeo(md1, md2) / 1000
#[1] 2323.15
You can also compute that like this
d1$lon <- 0
d2$lon <- 0
dm <- terra::distance(d1, d2, lonlat=TRUE, unit="km")
mean(dm)
#[1] 2323.158
Your data are longitude/latitude coordinates. You should not project them to compute distance, as that could lead to inaccurate results due to distortion.