计算一个点到英国海岸的最小距离

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:

除了英国，我一直在遵循这里显示的示例。由于这个原因，我使用的是英国EPSG:27700的CRS，它有以下投影字符串:

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

然而，我对wgs没有把握。84代码。目前我用:

"+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.

也尝试使用datum=OSGB36和+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 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.

我很难理解我在这里做错了什么。

1 个解决方案

#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.

你肯定没有做错什么。错误消息告诉你,一些全球海岸线shapefile坐标映射到正无穷。我不确定为什么这发生,但通常期望全球平面投影在一个特定的地区工作并不是一个好主意(尽管它工作在另一个问题…)。一个解决方案是将海岸线的shapefile剪切到特定投影的边界框中，然后转换剪切好的shapefile。可以在这里找到epss .27700的边界框。

# 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

所以最近的海岸线距离是32.2公里。我们也可以确定这是在海岸的什么地方

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