R 测量距海岸线的距离

Question

我有一组坐标：

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' 中的数据是说明性的，这些纬度/经度点可能都在陆地上或其他地方

Answer 1

到海岸线的距离可以通过下载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))

这只是到海岸线的距离。您还需要知道是内陆还是海上。为此，您需要一个单独的海形状文件：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 秒即可运行！

其余的同上一个答案。

Answer 2

有一种简单的方法可以使用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 创建于 2021-12-01 (v2.0.1)}

Answer 3

首先，您需要一个包含英国海岸线的文件。

您可以使用question 中提供的方法来了解某个点是在英国海岸线内还是在外。

然后，对于落入英国的点，您可以计算它们与海岸线点之间的Haversine distance，以确定它们是否在距海岸 3 英里的范围内。

Answer 4

为了更快地实现geosphere:::dist2Line，它使用purrr 进行高效循环，使用progress 作为进度条，从而保持Chris' first answer 的准确性，见下文：

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