Calculation of the minimum distance between a point and the coast in the UK

Question

I followed the example below, but for the UK. For this reason, I use CRS for UK EPSG: 27700, which has the following forecast line:

"+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'm not sure about the wgs.84 code. I am currently using:

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

Also tried to use datum = OSGB36 and + ellps = airy.

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 you try 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 did wrong here.

Answer 1

Afatik, you are not doing anything wrong. The error message states that some of the coordinates in the card file of the shapefile of the global coast correspond to Inf . I’m not sure why this is happening for sure, but, as a rule, waiting for a planar projection in a specific region to work on a global scale is not a good idea (although it worked on a different issue ...). One way is to crop the shapefile of the shoreline file into a bounding box for your specific projection, and then convert the cropped shapefile. The bounding box for the 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 

Thus, the nearest coastline is 32.2 km. We can also determine where on the coast it 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: # xy # 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 build 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")