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

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)
``````
1. 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

2. 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
``````
• The proper way is to convert your coordinates to UTM as @Kartograaf suggests Commented Sep 7, 2022 at 17:25

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

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

pointfile <- 'points.csv'

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
``````

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
• Both answers helped greatly so thank you and @kartograaf. I ended up doing: - download map - get raster with `as.raster` - project the raster with `raster::projectRaster` to UTM zone 55 - project my centroids with the same UTM zone using `sp::spTransform` - get the distances with `gdistance::shortestPath`
– Sam
Commented Mar 13, 2020 at 10:09

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:
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)
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")