计算英国某点与海岸之间的最小距离

Posted 2023-02-23

【中文标题】计算英国某点与海岸之间的最小距离

【英文标题】：calculating minimum distance between a point and the coast in the UK

【发布时间】：2014-12-09 16:54:02

【问题描述】：

我一直按照here 所示的示例进行操作，但针对的是英国。出于这个原因，我使用的是 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"

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

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

也尝试过使用 datum=OSGB36 和 +ellps=airy。

完整代码如下：

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

当尝试完成投影线时，会显示错误。

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

我无法理解我在这里做错了什么。

【参考方案1】：

AFAICT 你没有做错任何事。错误消息告诉您，全球海岸线 shapefile 中的某些坐标映射到 Inf。我不确定为什么会发生这种情况，但通常期望特定区域中的平面投影在全球范围内起作用并不是一个好主意（尽管它确实在另一个问题中起作用......）。一种解决方法是将海岸线 shapefile 剪辑到特定投影的边界框，然后转换剪辑的 shapefile。 EPSG.27700的边界框可以在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 公里。我们还可以确定这是在海岸的哪个位置

# 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 

最后我们可以绘制结果。您应该始终这样做。

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

【讨论】：

只是对此的跟进，是否可以将 MAD、MAD.proj 和 dist.27700 放置在一个循环中，因为我有大约 400K 个位置来检查到海岸或任何其他的距离想法如何做到这一点？

我倾向于将所有点放入一个 spatialPoints 对象（为方便起见，我们再次称其为 MAD），然后将其转换为 MAD.proj，然后执行以下操作：@ 987654330@。这将产生一个从每个点到最近海岸的距离向量。

感谢您的建议。我正在为将所有点放入单个 spatialPoints 对象的概念而苦苦挣扎。我尝试过使用 MULTIPOINT，但这似乎只接受大约 200 分。我也尝试将数据放入数据帧并在此数据帧上使用坐标函数，但随后失败MAD.proj <- spTransform(points,CRS(epsg.27700)) 和No transformation possible from NA reference system 不确定我错过了什么？

这是因为spTransform(sp,...) 将从一个投影转换为另一个。它必须知道起始 CRS，这意味着 sp 必须具有已定义的 CRS。如果您的点位于名为 df 的矩阵或数据框中，请尝试以下操作：MAD <- SpatialPoints(df,proj4string=CRS(wgs.84))。然后spTransform(...) 应该可以工作。