R测量距海岸线的距离
R measuring distance from a coastline
我有一组坐标:
d1 <- data_frame(
title = c("base1", "base2", "base3", "base4"),
lat = c(57.3, 58.8, 47.2, 57.8, 65.4, 56.7, 53.3),
long = c(0.4, 3.4, 3.5, 1.2, 1.5, 2.6, 2.7))
我想知道坐标是落在陆地上、海洋中还是海岸线内 3 英里处。坐标应该落在英国境内的某个地方,所以我知道我需要绘制英国的形状文件并将点绘制到它上面。
我只是不知道如何测量这些点是落在海里、陆地上还是离海岸 2 英里的地方。显然我可以通过查看地图来判断它们落在哪里,但我想要的是在我的数据集中添加另一列,因此它看起来像这样:
d2 <- data_frame(
title = c("base1", "base2", "base3", "base4", "base5", "base6", "base7"),
lat = c(57.3, 58.8, 47.2, 57.8, 65.4, 56.7, 53.3),
long = c(0.4, 3.4, 3.5, 1.2, 1.5, 2.6, 2.7),
where = c("land", "land", "sea", "coast", "land", "sea", "coast"))
- 注意d2列'where'中的数据是说明性的,这些lat/long点可能都在陆地上,或者其他地方
首先您需要一个包含英国海岸线的文件。
您可以使用此question中提供的方法来了解某个点是落在英国海岸线内还是外部。
然后,对于落入英国的点,您可以计算它们与海岸线点之间的 Haversine distance,以确定它们是否在距离海岸 3 英里的范围内。
到海岸线的距离可以通过下载openstreetmap
海岸线数据来计算。然后,您可以使用 geosphere::dist2Line
获取点到海岸线的距离。
我注意到您的示例点之一在法国,因此您可能需要将海岸线数据扩展到英国以外(可以通过调整边界框的范围来完成)。
library(tidyverse)
library(sf)
library(geosphere)
library(osmdata)
#get initial data frame
d1 <- data_frame(
title = c("base1", "base2", "base3", "base4",
"base5", "base6", "base7"),
lat = c(57.3, 58.8, 47.2, 57.8, 65.4, 56.7, 53.3),
long = c(0.4, 3.4, 3.5, 1.2, 1.5, 2.6, 2.7))
# convert to sf object
d1_sf <- d1 %>% st_as_sf(coords = c('long','lat')) %>%
st_set_crs(4326)
# get bouding box for osm data download (England) and
# download coastline data for this area
osm_box <- getbb (place_name = "England") %>%
opq () %>%
add_osm_feature("natural", "coastline") %>%
osmdata_sf()
# use dist2Line from geosphere - only works for WGS84
#data
dist <- geosphere::dist2Line(p = st_coordinates(d1_sf),
line =
st_coordinates(osm_box$osm_lines)[,1:2])
#combine initial data with distance to coastline
df <- cbind( d1 %>% rename(y=lat,x=long),dist) %>%
mutate(miles=distance/1609)
# title y x distance lon lat miles
#1 base1 57.3 0.4 219066.40 -2.137847 55.91706 136.15065
#2 base2 58.8 3.4 462510.28 -2.137847 55.91706 287.45201
#3 base3 47.2 3.5 351622.34 1.193198 49.96737 218.53470
#4 base4 57.8 1.2 292210.46 -2.137847 55.91706 181.60998
#5 base5 65.4 1.5 1074644.00 -2.143168 55.91830 667.89559
#6 base6 56.7 2.6 287951.93 -1.621963 55.63143 178.96329
#7 base7 53.3 2.7 92480.24 1.651836 52.76027 57.47684
#plot
p <- ggplot() +
geom_sf(data=osm_box$osm_lines) +
geom_sf(data=d1_sf) +
geom_segment(data=df,aes(x=x,y=y,xend=lon,yend=lat))
那只是到海岸线的距离。您还需要知道它是内陆还是海上。为此,您需要一个单独的海洋 shapefile:http://openstreetmapdata.com/data/water-polygons 并查看您的点的每个点是否位于海中。
#read in osm water polygon data
sea <- read_sf('water_polygons.shp')
#get get water polygons that intersect our points
in_sea <- st_intersects(d1_sf,sea) %>% as.data.frame()
#join back onto original dataset
df %>% mutate(row = row_number()) %>%
#join on in_sea data
left_join(in_sea,by=c('row'='row.id')) %>%
mutate(in_sea = if_else(is.na(col.id),F,T)) %>%
#categorise into 'sea', 'coast' or 'land'
mutate(where = case_when(in_sea == T ~ 'Sea',
in_sea == F & miles <=3 ~ 'Coast',
in_sea == F ~ 'Land'))
# title y x distance lon lat miles row col.id in_sea where
#1 base1 57.3 0.4 219066.40 -2.137847 55.91706 136.15065 1 24193 TRUE Sea
#2 base2 58.8 3.4 462510.28 -2.137847 55.91706 287.45201 2 24194 TRUE Sea
#3 base3 47.2 3.5 351622.34 1.193198 49.96737 218.53470 3 NA FALSE Land
#4 base4 57.8 1.2 292210.46 -2.137847 55.91706 181.60998 4 24193 TRUE Sea
#5 base5 65.4 1.5 1074644.00 -2.143168 55.91830 667.89559 5 25417 TRUE Sea
#6 base6 56.7 2.6 287951.93 -1.621963 55.63143 178.96329 6 24193 TRUE Sea
#7 base7 53.3 2.7 92480.24 1.651836 52.76027 57.47684 7 24143 TRUE Sea
ggplot() +
geom_sf(data=osm_box$osm_lines) +
geom_sf(data=d1_sf) +
geom_segment(data=df,aes(x=x,y=y,xend=lon,yend=lat)) +
ggrepel::geom_text_repel(data=df,
aes(x=x,y=y,label=paste0(where,'\n',round(miles,0),'miles')),size=2)
2018 年 8 月 16 日更新
由于您询问了一种专门使用 shapefile 的方法,因此我已在此处下载了该方法:openstreetmapdata。com/data/coastlines 我将使用它来执行与上述相同的方法。
clines <- read_sf('lines.shp') #path to shapefile
接下来,我创建了一个自定义边界框,这样我们就可以减小 shapefile 的大小,使其仅包含合理靠近这些点的海岸线。
# create bounding box surrounding points
bbox <- st_bbox(d1_sf)
# write a function that takes the bbox around our points
# and expands it by a given amount of metres.
expand_bbox <- function(bbox,metres_x,metres_y){
box_centre <- bbox %>% st_as_sfc() %>%
st_transform(crs = 32630) %>%
st_centroid() %>%
st_transform(crs = 4326) %>%
st_coordinates()
bbox['xmin'] <- bbox['xmin'] - (metres_x / 6370000) * (180 / pi) / cos(bbox['xmin'] * pi/180)
bbox['xmax'] <- bbox['xmax'] + (metres_x / 6370000) * (180 / pi) / cos(bbox['xmax'] * pi/180)
bbox['ymin'] <- bbox['ymin'] - (metres_y / 6370000) * (180 / pi)
bbox['ymax'] <- bbox['ymax'] + (metres_y / 6370000) * (180 / pi)
bbox['xmin'] <- ifelse(bbox['xmin'] < -180, bbox['xmin'] + 360, bbox['xmin'])
bbox['xmax'] <- ifelse(bbox['xmax'] > 180, bbox['xmax'] - 360, bbox['xmax'])
bbox['ymin'] <- ifelse(bbox['ymin'] < -90, (bbox['ymin'] + 180)*-1, bbox['ymin'])
bbox['ymax'] <- ifelse(bbox['ymax'] > 90, (bbox['ymax'] + 180)*-1, bbox['ymax'])
return(bbox)
}
# expand the bounding box around our points by 300 miles in x and 100 #miles in y direction to make nice shaped box.
bbox <- expand_bbox(bbox,metres_x=1609*200, metres_y=1609*200) %>% st_as_sfc
# get only the parts of the coastline that are within our bounding box
clines2 <- st_intersection(clines,bbox)
现在我在这里使用了 dist2Line 函数,因为它很准确,并且可以为您提供它正在测量的海岸线上的点,这有助于检查错误。缺点是,对于我们相当大的海岸线文件来说,它非常慢。
运行 这花了我 8 分钟:
dist <- geosphere::dist2Line(p = st_coordinates(d1_sf),
line = as(clines2,'Spatial'))
#combine initial data with distance to coastline
df <- cbind( d1 %>% rename(y=lat,x=long),dist) %>%
mutate(miles=distance/1609)
df
# title y x distance lon lat ID miles
#1 base1 57.3 0.4 131936.70 -1.7711149 57.46995 4585 81.99919
#2 base2 58.8 3.4 98886.42 4.8461433 59.28235 179 61.45831
#3 base3 47.2 3.5 340563.02 0.3641618 49.43811 4199 211.66129
#4 base4 57.8 1.2 180110.10 -1.7670712 57.50691 4584 111.93915
#5 base5 65.4 1.5 369550.43 6.2494627 62.81381 9424 229.67709
#6 base6 56.7 2.6 274230.37 5.8635346 58.42913 24152 170.43528
#7 base7 53.3 2.7 92480.24 1.6518358 52.76027 4639 57.47684
剧情:
ggplot() +
geom_sf(data=clines2) +
geom_sf(data=bbox,fill=NA)+
geom_sf(data=d1_sf) +
geom_segment(data=df,aes(x=x,y=y,xend=lon,yend=lat))
如果您不介意精度的轻微损失(结果与您的数据相差约 0.3%),并且不关心它测量到的海岸线上的确切位置,您可以测量距离到多边形:
# make data into polygons
clines3 <- st_intersection(clines,bbox) %>%
st_cast('POLYGON')
#use rgeos::gDistance to calculate distance to nearest polygon
#need to change projection (I used UTM30N) to use gDistance
dist2 <- apply(rgeos::gDistance(as(st_transform(d1_sf,32630), 'Spatial'),
as(st_transform(clines3,32630),'Spatial'),
byid=TRUE),2,min)
df2 <- cbind( d1 %>% rename(y=lat,x=long),dist2) %>%
mutate(miles=dist2/1609)
df2
# title y x dist2 miles
#1 base1 57.3 0.4 131917.62 81.98733
#2 base2 58.8 3.4 99049.22 61.55949
#3 base3 47.2 3.5 341015.26 211.94236
#4 base4 57.8 1.2 180101.47 111.93379
#5 base5 65.4 1.5 369950.32 229.92562
#6 base6 56.7 2.6 274750.17 170.75834
#7 base7 53.3 2.7 92580.16 57.53894
相比之下,这只用了 8 秒 运行!
其余同上一个回答
为了更快地实现 geosphere:::dist2Line
,使用 purrr
进行高效循环,使用 progress
作为进度条,从而保持 的准确性,请参见下文:
library(geosphere)
library(purr)
library(progress)
spDistPoint2Line <- function (p, line, distfun)
{
## rewrite of internal function from geosphere
test <- !sp::is.projected(line)
if (!isTRUE(test)) {
if (is.na(test)) {
warning("Coordinate reference system of SpatialPolygons object is not set. Assuming it is degrees (longitude/latitude)!")
}
else {
stop("Points are projected. They should be in degrees (longitude/latitude)")
}
}
x <- line@lines
n <- length(x)
res <- matrix(nrow = nrow(p), ncol = 3)
colnames(res) <- c("distance", "lon", "lat")
line_coords <- map(x, ~(map(.@Lines, ~(.@coords)))) #basically an unlist
pb <- progress_bar$new(
total = length(line_coords),
format = "(:spin) :current of :total, :percent, eta: :eta"
)
res[] <- Inf
result <- reduce(
.x = line_coords,
.init = res,
.f = function(res, crd){
pb$tick()
crd <- crd[[1]]
r <- dist2Line(p, crd, distfun) # have to live without ID
k <- r[, 1] < res[, 1]
res[k, ] <- r[k, ]
return(res)
}
)
return(result)
}
dist2Line <- function (p, line, distfun = distGeo)
{
p <- geosphere:::.pointsToMatrix(p)
if (inherits(line, "SpatialPolygons")) {
line <- methods::as(line, "SpatialLines")
}
if (inherits(line, "SpatialLines")) {
return(spDistPoint2Line(p, line, distfun))
}
line <- geosphere:::.pointsToMatrix(line)
line1 <- line[-nrow(line), , drop = FALSE]
line2 <- line[-1, , drop = FALSE]
seglength <- distfun(line1, line2)
res <-
p %>%
array_branch(1) %>%
map(
function(xy){
crossdist <- abs(dist2gc(line1, line2, xy))
trackdist1 <- alongTrackDistance(line1, line2, xy)
trackdist2 <- alongTrackDistance(line2, line1, xy)
mintrackdist <- pmin(trackdist1, trackdist2)
maxtrackdist <- pmax(trackdist1, trackdist2)
crossdist[maxtrackdist >= seglength] <- NA
nodedist <- distfun(xy, line)
warnopt = getOption("warn")
options(warn = -1)
distmin1 <- min(nodedist, na.rm = TRUE)
distmin2 <- min(crossdist, na.rm = TRUE)
options(warn = warnopt)
if (distmin1 <= distmin2) {
j <- which.min(nodedist)
return(c(distmin1, line[j, ]))
}
else {
j <- which.min(crossdist)
if (trackdist1[j] < trackdist2[j]) {
bear <- bearing(line1[j, ], line2[j, ])
pt <- destPoint(line1[j, ], bear, mintrackdist[j])
return(c(crossdist[j], pt))
}
else {
bear <- bearing(line2[j, ], line1[j, ])
pt <- destPoint(line2[j, ], bear, mintrackdist[j])
return(c(crossdist[j], pt))
}
}
}
) %>%
simplify %>%
matrix(ncol = 3, byrow = TRUE)
colnames(res) <- c("distance", "lon", "lat")
return(res)
}
有一种简单的方法可以使用 ggOceanMaps package
library(ggOceanMaps)
library(dplyr)
d1 <- data.frame(
title = c("base1", "base2", "base3", "base4", "base5", "base6", "base7"),
lat = c(57.3, 58.8, 47.2, 57.8, 65.4, 56.7, 53.3),
long = c(0.4, 3.4, 3.5, 1.2, 1.5, 2.6, 2.7))
ggOceanMaps::dist2land(d1) %>%
mutate(where = ifelse(ldist == 0, "land",
ifelse(ldist < 100*1.852, "coast", "sea")))
#> Used long and lat as input coordinate column names in data
#> Using ArcticStereographic as land shapes.
#> Calculating distances with parallel processing...
#> title lat long ldist where
#> 1 base1 57.3 0.4 138.62996 coast
#> 2 base2 58.8 3.4 105.43322 coast
#> 3 base3 47.2 3.5 0.00000 land
#> 4 base4 57.8 1.2 189.78665 sea
#> 5 base5 65.4 1.5 382.64410 sea
#> 6 base6 56.7 2.6 289.76883 sea
#> 7 base7 53.3 2.7 99.18336 coast
The ggOceanMaps::dist2land
function returns distances in kilometers。我转换为海里并提高了限制以获得所有不同的类别。您的示例坐标会产生一个案例“陆地”,其余案例为“海洋”。
由 reprex package (v2.0.1)
创建于 2021-12-01
我有一组坐标:
d1 <- data_frame(
title = c("base1", "base2", "base3", "base4"),
lat = c(57.3, 58.8, 47.2, 57.8, 65.4, 56.7, 53.3),
long = c(0.4, 3.4, 3.5, 1.2, 1.5, 2.6, 2.7))
我想知道坐标是落在陆地上、海洋中还是海岸线内 3 英里处。坐标应该落在英国境内的某个地方,所以我知道我需要绘制英国的形状文件并将点绘制到它上面。
我只是不知道如何测量这些点是落在海里、陆地上还是离海岸 2 英里的地方。显然我可以通过查看地图来判断它们落在哪里,但我想要的是在我的数据集中添加另一列,因此它看起来像这样:
d2 <- data_frame(
title = c("base1", "base2", "base3", "base4", "base5", "base6", "base7"),
lat = c(57.3, 58.8, 47.2, 57.8, 65.4, 56.7, 53.3),
long = c(0.4, 3.4, 3.5, 1.2, 1.5, 2.6, 2.7),
where = c("land", "land", "sea", "coast", "land", "sea", "coast"))
- 注意d2列'where'中的数据是说明性的,这些lat/long点可能都在陆地上,或者其他地方
首先您需要一个包含英国海岸线的文件。
您可以使用此question中提供的方法来了解某个点是落在英国海岸线内还是外部。
然后,对于落入英国的点,您可以计算它们与海岸线点之间的 Haversine distance,以确定它们是否在距离海岸 3 英里的范围内。
到海岸线的距离可以通过下载openstreetmap
海岸线数据来计算。然后,您可以使用 geosphere::dist2Line
获取点到海岸线的距离。
我注意到您的示例点之一在法国,因此您可能需要将海岸线数据扩展到英国以外(可以通过调整边界框的范围来完成)。
library(tidyverse)
library(sf)
library(geosphere)
library(osmdata)
#get initial data frame
d1 <- data_frame(
title = c("base1", "base2", "base3", "base4",
"base5", "base6", "base7"),
lat = c(57.3, 58.8, 47.2, 57.8, 65.4, 56.7, 53.3),
long = c(0.4, 3.4, 3.5, 1.2, 1.5, 2.6, 2.7))
# convert to sf object
d1_sf <- d1 %>% st_as_sf(coords = c('long','lat')) %>%
st_set_crs(4326)
# get bouding box for osm data download (England) and
# download coastline data for this area
osm_box <- getbb (place_name = "England") %>%
opq () %>%
add_osm_feature("natural", "coastline") %>%
osmdata_sf()
# use dist2Line from geosphere - only works for WGS84
#data
dist <- geosphere::dist2Line(p = st_coordinates(d1_sf),
line =
st_coordinates(osm_box$osm_lines)[,1:2])
#combine initial data with distance to coastline
df <- cbind( d1 %>% rename(y=lat,x=long),dist) %>%
mutate(miles=distance/1609)
# title y x distance lon lat miles
#1 base1 57.3 0.4 219066.40 -2.137847 55.91706 136.15065
#2 base2 58.8 3.4 462510.28 -2.137847 55.91706 287.45201
#3 base3 47.2 3.5 351622.34 1.193198 49.96737 218.53470
#4 base4 57.8 1.2 292210.46 -2.137847 55.91706 181.60998
#5 base5 65.4 1.5 1074644.00 -2.143168 55.91830 667.89559
#6 base6 56.7 2.6 287951.93 -1.621963 55.63143 178.96329
#7 base7 53.3 2.7 92480.24 1.651836 52.76027 57.47684
#plot
p <- ggplot() +
geom_sf(data=osm_box$osm_lines) +
geom_sf(data=d1_sf) +
geom_segment(data=df,aes(x=x,y=y,xend=lon,yend=lat))
那只是到海岸线的距离。您还需要知道它是内陆还是海上。为此,您需要一个单独的海洋 shapefile:http://openstreetmapdata.com/data/water-polygons 并查看您的点的每个点是否位于海中。
#read in osm water polygon data
sea <- read_sf('water_polygons.shp')
#get get water polygons that intersect our points
in_sea <- st_intersects(d1_sf,sea) %>% as.data.frame()
#join back onto original dataset
df %>% mutate(row = row_number()) %>%
#join on in_sea data
left_join(in_sea,by=c('row'='row.id')) %>%
mutate(in_sea = if_else(is.na(col.id),F,T)) %>%
#categorise into 'sea', 'coast' or 'land'
mutate(where = case_when(in_sea == T ~ 'Sea',
in_sea == F & miles <=3 ~ 'Coast',
in_sea == F ~ 'Land'))
# title y x distance lon lat miles row col.id in_sea where
#1 base1 57.3 0.4 219066.40 -2.137847 55.91706 136.15065 1 24193 TRUE Sea
#2 base2 58.8 3.4 462510.28 -2.137847 55.91706 287.45201 2 24194 TRUE Sea
#3 base3 47.2 3.5 351622.34 1.193198 49.96737 218.53470 3 NA FALSE Land
#4 base4 57.8 1.2 292210.46 -2.137847 55.91706 181.60998 4 24193 TRUE Sea
#5 base5 65.4 1.5 1074644.00 -2.143168 55.91830 667.89559 5 25417 TRUE Sea
#6 base6 56.7 2.6 287951.93 -1.621963 55.63143 178.96329 6 24193 TRUE Sea
#7 base7 53.3 2.7 92480.24 1.651836 52.76027 57.47684 7 24143 TRUE Sea
ggplot() +
geom_sf(data=osm_box$osm_lines) +
geom_sf(data=d1_sf) +
geom_segment(data=df,aes(x=x,y=y,xend=lon,yend=lat)) +
ggrepel::geom_text_repel(data=df,
aes(x=x,y=y,label=paste0(where,'\n',round(miles,0),'miles')),size=2)
2018 年 8 月 16 日更新
由于您询问了一种专门使用 shapefile 的方法,因此我已在此处下载了该方法:openstreetmapdata。com/data/coastlines 我将使用它来执行与上述相同的方法。
clines <- read_sf('lines.shp') #path to shapefile
接下来,我创建了一个自定义边界框,这样我们就可以减小 shapefile 的大小,使其仅包含合理靠近这些点的海岸线。
# create bounding box surrounding points
bbox <- st_bbox(d1_sf)
# write a function that takes the bbox around our points
# and expands it by a given amount of metres.
expand_bbox <- function(bbox,metres_x,metres_y){
box_centre <- bbox %>% st_as_sfc() %>%
st_transform(crs = 32630) %>%
st_centroid() %>%
st_transform(crs = 4326) %>%
st_coordinates()
bbox['xmin'] <- bbox['xmin'] - (metres_x / 6370000) * (180 / pi) / cos(bbox['xmin'] * pi/180)
bbox['xmax'] <- bbox['xmax'] + (metres_x / 6370000) * (180 / pi) / cos(bbox['xmax'] * pi/180)
bbox['ymin'] <- bbox['ymin'] - (metres_y / 6370000) * (180 / pi)
bbox['ymax'] <- bbox['ymax'] + (metres_y / 6370000) * (180 / pi)
bbox['xmin'] <- ifelse(bbox['xmin'] < -180, bbox['xmin'] + 360, bbox['xmin'])
bbox['xmax'] <- ifelse(bbox['xmax'] > 180, bbox['xmax'] - 360, bbox['xmax'])
bbox['ymin'] <- ifelse(bbox['ymin'] < -90, (bbox['ymin'] + 180)*-1, bbox['ymin'])
bbox['ymax'] <- ifelse(bbox['ymax'] > 90, (bbox['ymax'] + 180)*-1, bbox['ymax'])
return(bbox)
}
# expand the bounding box around our points by 300 miles in x and 100 #miles in y direction to make nice shaped box.
bbox <- expand_bbox(bbox,metres_x=1609*200, metres_y=1609*200) %>% st_as_sfc
# get only the parts of the coastline that are within our bounding box
clines2 <- st_intersection(clines,bbox)
现在我在这里使用了 dist2Line 函数,因为它很准确,并且可以为您提供它正在测量的海岸线上的点,这有助于检查错误。缺点是,对于我们相当大的海岸线文件来说,它非常慢。
运行 这花了我 8 分钟:
dist <- geosphere::dist2Line(p = st_coordinates(d1_sf),
line = as(clines2,'Spatial'))
#combine initial data with distance to coastline
df <- cbind( d1 %>% rename(y=lat,x=long),dist) %>%
mutate(miles=distance/1609)
df
# title y x distance lon lat ID miles
#1 base1 57.3 0.4 131936.70 -1.7711149 57.46995 4585 81.99919
#2 base2 58.8 3.4 98886.42 4.8461433 59.28235 179 61.45831
#3 base3 47.2 3.5 340563.02 0.3641618 49.43811 4199 211.66129
#4 base4 57.8 1.2 180110.10 -1.7670712 57.50691 4584 111.93915
#5 base5 65.4 1.5 369550.43 6.2494627 62.81381 9424 229.67709
#6 base6 56.7 2.6 274230.37 5.8635346 58.42913 24152 170.43528
#7 base7 53.3 2.7 92480.24 1.6518358 52.76027 4639 57.47684
剧情:
ggplot() +
geom_sf(data=clines2) +
geom_sf(data=bbox,fill=NA)+
geom_sf(data=d1_sf) +
geom_segment(data=df,aes(x=x,y=y,xend=lon,yend=lat))
如果您不介意精度的轻微损失(结果与您的数据相差约 0.3%),并且不关心它测量到的海岸线上的确切位置,您可以测量距离到多边形:
# make data into polygons
clines3 <- st_intersection(clines,bbox) %>%
st_cast('POLYGON')
#use rgeos::gDistance to calculate distance to nearest polygon
#need to change projection (I used UTM30N) to use gDistance
dist2 <- apply(rgeos::gDistance(as(st_transform(d1_sf,32630), 'Spatial'),
as(st_transform(clines3,32630),'Spatial'),
byid=TRUE),2,min)
df2 <- cbind( d1 %>% rename(y=lat,x=long),dist2) %>%
mutate(miles=dist2/1609)
df2
# title y x dist2 miles
#1 base1 57.3 0.4 131917.62 81.98733
#2 base2 58.8 3.4 99049.22 61.55949
#3 base3 47.2 3.5 341015.26 211.94236
#4 base4 57.8 1.2 180101.47 111.93379
#5 base5 65.4 1.5 369950.32 229.92562
#6 base6 56.7 2.6 274750.17 170.75834
#7 base7 53.3 2.7 92580.16 57.53894
相比之下,这只用了 8 秒 运行!
其余同上一个回答
为了更快地实现 geosphere:::dist2Line
,使用 purrr
进行高效循环,使用 progress
作为进度条,从而保持
library(geosphere)
library(purr)
library(progress)
spDistPoint2Line <- function (p, line, distfun)
{
## rewrite of internal function from geosphere
test <- !sp::is.projected(line)
if (!isTRUE(test)) {
if (is.na(test)) {
warning("Coordinate reference system of SpatialPolygons object is not set. Assuming it is degrees (longitude/latitude)!")
}
else {
stop("Points are projected. They should be in degrees (longitude/latitude)")
}
}
x <- line@lines
n <- length(x)
res <- matrix(nrow = nrow(p), ncol = 3)
colnames(res) <- c("distance", "lon", "lat")
line_coords <- map(x, ~(map(.@Lines, ~(.@coords)))) #basically an unlist
pb <- progress_bar$new(
total = length(line_coords),
format = "(:spin) :current of :total, :percent, eta: :eta"
)
res[] <- Inf
result <- reduce(
.x = line_coords,
.init = res,
.f = function(res, crd){
pb$tick()
crd <- crd[[1]]
r <- dist2Line(p, crd, distfun) # have to live without ID
k <- r[, 1] < res[, 1]
res[k, ] <- r[k, ]
return(res)
}
)
return(result)
}
dist2Line <- function (p, line, distfun = distGeo)
{
p <- geosphere:::.pointsToMatrix(p)
if (inherits(line, "SpatialPolygons")) {
line <- methods::as(line, "SpatialLines")
}
if (inherits(line, "SpatialLines")) {
return(spDistPoint2Line(p, line, distfun))
}
line <- geosphere:::.pointsToMatrix(line)
line1 <- line[-nrow(line), , drop = FALSE]
line2 <- line[-1, , drop = FALSE]
seglength <- distfun(line1, line2)
res <-
p %>%
array_branch(1) %>%
map(
function(xy){
crossdist <- abs(dist2gc(line1, line2, xy))
trackdist1 <- alongTrackDistance(line1, line2, xy)
trackdist2 <- alongTrackDistance(line2, line1, xy)
mintrackdist <- pmin(trackdist1, trackdist2)
maxtrackdist <- pmax(trackdist1, trackdist2)
crossdist[maxtrackdist >= seglength] <- NA
nodedist <- distfun(xy, line)
warnopt = getOption("warn")
options(warn = -1)
distmin1 <- min(nodedist, na.rm = TRUE)
distmin2 <- min(crossdist, na.rm = TRUE)
options(warn = warnopt)
if (distmin1 <= distmin2) {
j <- which.min(nodedist)
return(c(distmin1, line[j, ]))
}
else {
j <- which.min(crossdist)
if (trackdist1[j] < trackdist2[j]) {
bear <- bearing(line1[j, ], line2[j, ])
pt <- destPoint(line1[j, ], bear, mintrackdist[j])
return(c(crossdist[j], pt))
}
else {
bear <- bearing(line2[j, ], line1[j, ])
pt <- destPoint(line2[j, ], bear, mintrackdist[j])
return(c(crossdist[j], pt))
}
}
}
) %>%
simplify %>%
matrix(ncol = 3, byrow = TRUE)
colnames(res) <- c("distance", "lon", "lat")
return(res)
}
有一种简单的方法可以使用 ggOceanMaps package
library(ggOceanMaps)
library(dplyr)
d1 <- data.frame(
title = c("base1", "base2", "base3", "base4", "base5", "base6", "base7"),
lat = c(57.3, 58.8, 47.2, 57.8, 65.4, 56.7, 53.3),
long = c(0.4, 3.4, 3.5, 1.2, 1.5, 2.6, 2.7))
ggOceanMaps::dist2land(d1) %>%
mutate(where = ifelse(ldist == 0, "land",
ifelse(ldist < 100*1.852, "coast", "sea")))
#> Used long and lat as input coordinate column names in data
#> Using ArcticStereographic as land shapes.
#> Calculating distances with parallel processing...
#> title lat long ldist where
#> 1 base1 57.3 0.4 138.62996 coast
#> 2 base2 58.8 3.4 105.43322 coast
#> 3 base3 47.2 3.5 0.00000 land
#> 4 base4 57.8 1.2 189.78665 sea
#> 5 base5 65.4 1.5 382.64410 sea
#> 6 base6 56.7 2.6 289.76883 sea
#> 7 base7 53.3 2.7 99.18336 coast
The ggOceanMaps::dist2land
function returns distances in kilometers。我转换为海里并提高了限制以获得所有不同的类别。您的示例坐标会产生一个案例“陆地”,其余案例为“海洋”。
由 reprex package (v2.0.1)
创建于 2021-12-01