calculating minimum distance between a point and the coast in the UK

Question

I have been following the example shown here, but for the UK. For this reason I am using the CRS for the UK of EPSG:27700, which has the following projection string:

"+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs"

However, I am unsure on the wgs.84 code to follow. Currently I am using:

"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

Have also tried to use the datum=OSGB36 and +ellps=airy as well.

The full code is as follows:

library(rgeos)
library(maptools)
library(rgdal)

epsg.27700 <- '+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000    +ellps=airy +datum=OSGB36 +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)) #have tried other shapefiles with the same issue
MAD   <- readWKT("POINT(-0.1830372 51.1197467)",p4s=CRS(wgs.84)) #Crawley, West Sussex
gDistance(MAD,coast)

[1] 0.28958
Warning messages:
1: In RGEOSDistanceFunc(spgeom1, spgeom2, byid, "rgeos_distance") :
  Spatial object 1 is not projected; GEOS expects planar coordinates
2: In RGEOSDistanceFunc(spgeom1, spgeom2, byid, "rgeos_distance") :
  Spatial object 2 is not projected; GEOS expects planar coordinates
3: In RGEOSDistanceFunc(spgeom1, spgeom2, byid, "rgeos_distance") :
  spgeom1 and spgeom2 have different proj4 strings

When trying to complete the projection line an error is displayed.

coast.proj <- spTransform(coast,CRS(Epsg.27700))

non finite transformation detected:
 [1] 111.01051  19.68378       Inf       Inf
Error in .spTransform_Line(input[[i]], to_args = to_args, from_args = from_args,  : 
 failure in Lines 22 Line 1 points 1
In addition: Warning message:
In .spTransform_Line(input[[i]], to_args = to_args, from_args = from_args,  :
 671 projected point(s) not finite

I am having trouble understanding what I have done wrong here.

Forerun

Answer 1

AFAICT you are not doing anything wrong. The error message is telling you that some of the coordinates in the global coastline shapefile map to Inf. I'm not sure why this is happening exactly, but generally expecting a planar projection in a specific region to work globally is not a good idea (although it did work in the other question...). One workaround is to clip the coastline shapefile to the bounding box for your specific projection, then transform the clipped shapefile. The bounding box for EPSG.27700 can be found here.

# use bounding box for epsg.27700
# found here: http://spatialreference.org/ref/epsg/osgb-1936-british-national-grid/
bbx <- readWKT("POLYGON((-7.5600 49.9600, 1.7800 49.9600, 1.7800 60.8400, -7.5600 60.8400, -7.5600 49.9600))",
               p4s=CRS(wgs.84)) 
coast <- gIntersection(coast,bbx)  # just coastlines within bbx
# now transformation works
coast.proj   <- spTransform(coast,CRS(epsg.27700))
MAD.proj     <- spTransform(MAD,CRS(epsg.27700))
dist.27700   <- gDistance(MAD.proj,coast.proj)   # distance in sq.meters
dist.27700
# [1] 32153.23

So the closest coastline is 32.2 km away. We can also identify where on the coast this is

# identify the closest point
th     <- seq(0,2*pi,len=1000)
circle <- cbind(1.00001*dist.wgs84*cos(th)+MAD$x,1.00001*dist.wgs84*sin(th)+MAD$y)
sp.circle <- SpatialLines(list(Lines(list(Line(circle)),ID="1")),proj4string=CRS(wgs.84))
sp.pts <- gIntersection(sp.circle,coast)
sp.pts
# SpatialPoints:
#            x        y
# 1 -0.2019854 50.83079
# 1 -0.1997009 50.83064
# Coordinate Reference System (CRS) arguments: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 

And finally we can plot the results. You should always do this.

# plot the results: plot is in CRS(wgs.84)
plot(coast, asp=1)
plot(bbx,add=T)
plot(MAD,add=T, col="red", pch=20)
plot(sp.circle,add=T, col="blue", lty=2)
lines(c(MAD$x,sp.pts$x),c(MAD$y,sp.pts$y),col="red")

Ostend

Just a follow up on this, is it possible to place the MAD, MAD.proj and dist.27700 within a loop as I have circa 400K locations to check the distance to the coast, or any other ideas how to do this? – Forerun

I'd be inclined to put all the points into a single spatialPoints object (let's call it MAD again for convenience), then transform that to MAD.proj, then do this: sapply(1:length(MAD.proj),function(i)gDistance(MAD.proj[i],coast.proj)). This will produce a vector of distances from each point to the nearest coast. – Ostend

Thank you for that advice. I am struggling with the concept of placing all the points into a single spatialPoints object. I have tried using MULTIPOINT, but this seems to only accept circa 200 points. I have also tried bringing the data into a dataframe and using the coordinates function on this dataframe, but the following then fails MAD.proj <- spTransform(points,CRS(epsg.27700)) with No transformation possible from NA reference system Not sure what I am missing? – Forerun

This is because spTransform(sp,...) transforms from one projection to another. It has to know the starting CRS, which means that sp has to have a defined CRS. If your points are in a matrix or data frame called df, try this: MAD <- SpatialPoints(df,proj4string=CRS(wgs.84)). Then spTransform(...) should work. – Ostend