从 Lat Lon 数据创建距离海岸 (km) 变量？

Posted 2023-02-23

【中文标题】从 Lat Lon 数据创建距离海岸 (km) 变量？

【英文标题】：Create Distance to Shore (km) Variable from Lat Lon data?

【发布时间】：2018-07-31 15:42:00

【问题描述】：

我有一个包含 3k + 个数据点的数据框，分布在墨西哥湾北部（这里我只提供 6 个）。我正在尝试为到岸的距离（公里）创建一个新变量。我有一个想要使用的形状文件 (gulf.shape)，但我不清楚如何使用。

这是一些数据

require(maptools)
require(sp)
library(rgdal)
library(lubridate)


df <- data.frame(Lat = c(26.84853, 28.38329, 28.00364, 
                     29.53840, 29.32030, 26.81622, 25.28146), 
             Lon = c(-96.55716, -94.29307, -91.21581, 
                     -88.42556, -84.20031, -83.89737, -82.95665))

然后我加载 shapefile（提供 here）。

gulf.shape <- "Shape\\stateshigh.shp"
gulf.shape <- maptools::readShapePoly(gulf.shape)

还有一个快速的图表来可视化我所拥有的。

plot(df$Lon, df$Lat,
     xlim = c(-97.5, -80.7), ylim = c(25, 30.5),
     xlab ="Latitude", ylab = "Longitude",
     pch = 20, col="red", cex=1.5)
par(new=T)
sp::plot(gulf.shape, add= T,
         xlim = c(-97.5, -80.7), ylim = c(25, 30.5),
         xlab ="Latitude", ylab = "Longitude",
         col = "gray")

我发现了一个堆栈溢出帖子 (here)，它允许我使用下面的代码得到答案。他们使用的形状文件是可用的here。

require(rgdal)   # for readOGR(...); loads package sp as well
require(rgeos)   # for gDistance(...)


require(parallel) # for detect cores
require(foreach)   # for foreach(...)
require(snow)      # for makeCluster(...)
require(doSNOW)    # for resisterDoSNOW(...)


wgs.84    <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
mollweide <- "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
sp.points <- SpatialPoints(df[,c("Lon","Lat")], proj4string=CRS(wgs.84))

coast  <- rgdal::readOGR(dsn=".",layer="ne_10m_coastline",p4s=wgs.84);str(coast)

coast.moll <- spTransform(coast,CRS(mollweide))

point.moll <- spTransform(sp.points,CRS(mollweide))

no_cores <- detectCores()

cl <- makeCluster(no_cores,type="SOCK")  # create a 4-processor cluster
registerDoSNOW(cl)                # register the cluster

get.dist.parallel <- function(n) 
  foreach(i=1:n, .combine=c, .packages="rgeos", .inorder=TRUE, 
          .export=c("point.moll","coast.moll")) %dopar% gDistance(point.moll[i],coast.moll)


df$Dis.to.SHORE <- get.dist.parallel(length(sp.points))

df$Dis.to.SHORE <- df$Dis.to.SHORE/1000

df

plot(coast)

points(sp.points,pch=20,col="red")

但是，我不理解 wgs.84 和 mollweide 中使用的 CRS 代码，这让我对使用此代码生成的数据感到不安。我还想只使用 gulf.shape 文件而不是整个世界，因为在前面提到的堆栈溢出帖子中有一些建议，这会更好。

所以我的问题是：

我得到的离岸距离值是否相当准确（即在 5 - 10m 内）？

如何修改代码以利用gulf.shape 文件而不是整个世界？

谁能解释一下CRS 代码或指出一个好的参考？

请注意，我使用并行计算来加快速度，因为我实际上有超过 6 个数据点。

【参考方案1】：

我将尝试回答您的第三个问题。 CRS（坐标参考系统）在映射中非常重要，因为它定义了点的坐标系统。这是一个有用的概述。 https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/OverviewCoordinateReferenceSystems.pdf

对于您的特定情况，当您更改为不同的 shapefile 时，您需要 (1) 找出适用于您的 shapefile (gulf.shape) 的 CRS。通常，它位于 shapefile 附带的 .prj 文件或元数据中。 (2) 选择适合您目标的 CRS。您正在计算距离，因此等距投影可能对您最有帮助。 (3)在计算距离前将原始crs变换为目标crs。 您引用的代码也是这样做的。世界 shapefile 带有 wgs84 crs；选择的目标 crs 是 mollweide；它使用 spTransform() 函数将 wgs84 转换为 mollweide。

另一方面，与您的第一个问题有关，您计算的准确性与您使用的 crs 有关，但也与您的 shapefile 的比例和点的精度（纬度/经度）有关。