calculating minimum distance between a point and the coast

Question

I am trying to get the minimum distance between a given point and the coast. My example is the distance of Madrid to the coast:

library(rgeos)
library(maptools)
coast <- readShapeLines("Natural_Earth_quick_start/10m_physical/ne_10m_coastline.shp")
MAD = readWKT("POINT(-3.716667 40.383333)")
gDistance(MAD,coast)
[1] 3.021808

I am having trouble understanding what gDistance() returns. The docs say it's in the units of the projection. Does that mean it is in latlong degrees? How can I convert that to kilometers?

Aruba

Can you provide a link to your map? As a start, type proj4string(coast) to see what the projection is. If there is a `+units=' string, then those are the units. – Leola

The map I'm using is the coastline from naturalearthdata.com/downloads/10m-physical-vectors. proj4string(coast) returns NA – Aruba

Have you looked at the definition of Hausdorff distance linked from the help page?. I get the idea that this answer is in "degrees" but I'm not at all convinced that is meaningful since the number of km in a "degree" varies from roughly 100 at the equator to zero at the poles – Forage

Answer 1

gDistance(...) returns the minimum Cartesian (Euclidean) distance between the point and feature set provided as arguments. Since your map is in long/lat coordinates, you get distance in "degrees", e.g.

d = sqrt { (long₁ - long₂)² + (lat₁ - lat₂)² }

where long and lat are in decimal degrees. As was pointed out, this doesn't mean much because converting to planar distance (say, km) depends on where you are. So we need to transform your data into a CRS which is approximately planar in the region of interest. It turns out that the appropriate CRS for Spain is EPSG-2062. The projection string for EPSG-2062 is:

+proj=lcc +lat_1=40 +lat_0=40 +lon_0=0 +k_0=0.9988085293 +x_0=600000 +y_0=600000 +a=6378298.3 +b=6356657.142669561 +pm=madrid +units=m +no_defs 

which has +units=m (meters). So we need to reproject both the point (MAD) and the borders to EPSG-2062.

library(rgeos)
library(maptools)

epsg.2062 <- "+proj=lcc +lat_1=40 +lat_0=40 +lon_0=0 +k_0=0.9988085293 +x_0=600000 +y_0=600000 +a=6378298.3 +b=6356657.142669561 +pm=madrid +units=m +no_defs"
wgs.84    <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

coast <- readShapeLines("ne_10m_coastline",CRS(wgs.84))
MAD   <- readWKT("POINT(-3.716667 40.383333)",p4s=CRS(wgs.84))
gDistance(MAD,coast)            # WGS-84 (long/lat) projection, units in "degrees"
# [1] 3.021808

coast.proj <- spTransform(coast,CRS(epsg.2062))
MAD.proj   <- spTransform(MAD,CRS(epsg.2062))
gDistance(MAD.proj,coast.proj)  #EPSG-2062 projection, units are in meters.
# [1] 305171.2

So the minimum distance is ~305.2km.

Finally, note that your coastline file has all the coastlines of the world, so this is the minimum distance to some coastline, not necessarily the Spanish coast (although in this case it does turn out to be on the northern coast of Spain). If your reference point was very near the border with Portugal, the the nearest coastal point would be to the western coast of Portugal.

Leola

Thanks!, that is exactly what I wanted. How do you know that it's the northern coast? – Aruba

The short answer is that I looked at a map, but looking again, Madrid might be closer to the Mediterranean coast near Castellon. There is a way to identify the coordinates of the closest point but it's complicated. If you need to know that you should post as a separate question. – Leola

This looks great but what projection should I use when calculating the distance of points (in km) to nearest shore at a global scale? I have ~250 points! – Strutting

This answer doesn't generalize well. – Yorgo