且构网

分享程序员开发的那些事...
且构网 - 分享程序员编程开发的那些事

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

更新时间:2023-02-10 11:16:46

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