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

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 EPSG. 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")

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

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