What is the proper way to calculate distance in km from degrees

Question

For a genetic analysis of Australian humpback dolphins I am currently trying to assess whether or not genetic differentiation correlates with geographic distance in water. For this I have done the following:

1. Get bathymetric data using the R package marmap

getNOAA.bathy(lon1 = 145, lon2=154, lat1 = -28, lat2 = -17, resolution = 1)

2. Convert it to a raster and subsequently to a transition matrix with low cost for moving in water and prohibitively high cost for traversing land

3. Get the shortest distance in water between each pair of points using gdistance::shortestPath, and converting this into a distance matrix

I am now faced with the following: The lengths I am getting from shortestPath are in the same dimension as the coordinate system I am using, which is decimal degrees. So my question is: What is the proper way to convert this distance matrix into km?

I understand these degrees must be projected somehow but I am currently failing at figuring out how exactly I can do this.

my centroids are:

data.frame(lon = c(146.1333, 146.8698, 148.7083, 150.8709, 151.3147, 152.9791, 153.0279, 153.1816), lat = c(-18.33419, -19.22425, -20.55407, -23.43784, -23.84044, -25.26728, -25.81851, -27.35664))

and the corresponding distance matrix is:

1.3629804 
3.9769987  2.7628686
7.7416532  6.5480830  4.0348631
8.1950482  6.9809181  4.4882582  0.8048330
10.4037986  9.1896685  6.6970086  3.0371065  2.3537371
10.9785875  9.7644574  7.2650158  3.6118972  2.9285279  0.6491764
12.7441285 11.5573399  9.0441220  5.3774364  4.7008506  2.4420589  1.7487997

Answer 1

Convert your points to UTM first using Rgdal, for example:

dir <- 'C:/your/path'
setwd(dir)#set working directory

pointfile <- 'points.csv'
data <- read.csv(pointfile)#read in data

library(rgdal)

#load in points as lat long
coords <- SpatialPoints(cbind(data$lon, data$lat), proj4string = CRS("+proj=longlat"))

#convert points to WGS84 UTM zone 56S (western AUS)
coord.UTM <- spTransform(coords, CRS("+init=epsg:32756")) 

Or, you can use the package 'raster' to compute the distances in meters directly from the latlong points, though the estimates will vary slightly due to lack of precision in the lat and long coordinates.

library(raster)

dists <- pointDistance(coord.UTM, lonlat = FALSE)

dists_lonlat <- pointDistance(coords, lonlat = TRUE) #specify for lonlat

Answer 2

I'm not familiar with the marmap package and what it returns to the user. But when calculating distances, I think it's safer to have your source data projected to a CRS that is already in meters. So, what I would try is:

1. Get the data using the getNOAA.bathy()
2. Convert the data to raster using the as.raster() from the same package
3. Project it to a CRS suitable for your area
4. Calculate the distances

Answer 3

This is a needed addition to this post and a good example of how "things are always changing". With the pending release of sf 1.0 the interface with the s2 geometry library and changes in GDAL and PROJ handle projections, the flat earth model will no longer be supported for geographic coordinate operations. A spherical geometry will now be used for operations on geographic projections.

The most important thing to take note of is that how projections will be handled in sp, raster and sf are changing along with the changes in GDAL/PROJ moving to a full geodetic library and WKT2 coordinate reference systems. You will see this now with errors associated with certain definitions in proj4 and dwindling support for PROJ4. Here is a quick example of what WKT2 projection descriptions look like based on the familiar proj4 string.

sf::st_crs("+proj=longlat +datum=NAD83 +no_defs")

Coordinate Reference System:
  User input: +proj=longlat +datum=NAD83 +no_defs 
  wkt:
GEOGCRS["unknown",
    DATUM["North American Datum 1983",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]],
        ID["EPSG",6269]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8901]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]]

In sf versions < 1.0, ellipsoidal geometries were supported through the lwgeom package. The function lwgeom::st_geod_distance will return spheroid distances (meters) for data in geographic projections. There is also support for areas lwgeom::st_geod_area and lengths lwgeom::st_geod_length. These functions are automatically called by the associated sf functions if sf::st_is_longlat(nc) == TRUE so, one can just use st_dist.

library(sf)
library(lwgeom)
nc = read_sf(system.file("gpkg/nc.gpkg", package="sf")) 
lwgeom::st_geod_distance(nc[1:10,],nc[1:10,])
sf::st_distance(nc[1:10,],nc[1:10,])

There is great information on these changes here and here. If you want to play around with the new spherical geometry implementations in sf you can install the development versions of s2 and sf using (although the current dev version of sf is failing build but should resolved soon):

remotes::install_github("r-spatial/s2")
remotes::install_github("r-spatial/sf", ref = "s2")