AFAICT 你没有做错任何事。错误消息告诉您全球海岸线形状文件中的某些坐标映射到Inf
。我不确定为什么会发生这种情况,但通常期望特定区域的平面投影在全球范围内工作并不是一个好主意(尽管它确实在另一个问题中起作用......)。一种解决方法是将海岸线形状文件剪切到特定投影的边界框,然后转换剪切的形状文件。可以找到 EPSG.27700 的边界框here http://spatialreference.org/ref/epsg/osgb-1936-british-national-grid/.
# 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
最后我们可以绘制结果。你应该always做这个。
# 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")