计算英国点与海岸之间的最小距离 [英] calculating minimum distance between a point and the coast in the UK

问题描述

我一直在此处，但适用于英国.因此，我在EPSG:27700的英国使用CRS，其投影字符串如下:

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"

但是，我不确定要遵循的wgs.84代码.目前我正在使用:

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"

还尝试过使用datum = OSGB36和+ ellps = airy.

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

完整代码如下:

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.

推荐答案

AFAICT您没有做错任何事情.该错误消息告诉您全局海岸线shapefile中的某些坐标映射到Inf.我不确定为什么会发生这种情况，但是通常期望在特定区域进行平面投影才能在全球范围内工作不是一个好主意(尽管它确实在其他问题中起作用了...).一种解决方法是将海岸线shapefile剪切到特定投影的边界框，然后转换剪切的shapefile.可以在此处找到EPSG.27700的边框.

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

所以最近的海岸线在32.2公里外.我们还可以确定这是在海岸上的什么地方

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

这篇关于计算英国点与海岸之间的最小距离的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！