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英里处.坐标应该在英国的某个地方，所以我知道我需要绘制英国的形状文件并在其上绘制点.

I would like to know whether the coordinates fall on land, in the sea, or are 3 miles inside a coastline. The coordinates should fall somewhere within the UK, so I know that I need to draw a shape file of the UK and plot the points onto it.

我只是不知道如何测量这些点是落在海洋，陆地还是在距海岸2英里的地方.显然，从查看地图可以看出它们落在哪里，但是我想在数据集中添加另一列，如下所示:

I just don't know how to measure whether the points fall in the sea, land or 2 miles from the coast. Obviously I can tell from looking at a map where they fall, but what I would like is to add another column to my data set so it looks like this:

 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"中的数据是说明性的，这些经/纬度点可能全部在陆地上，或者其他地方
• note that the data in d2 column 'where' is illustrative, these lat/long points may all be on land, or something else

推荐答案

可以通过下载openstreetmap海岸线数据来计算到海岸线的距离.然后，您可以使用geosphere::dist2Line获取从您的点到海岸线的距离.

Distance to the coastline can be calculated by downloading openstreetmap coastline data. You can then use geosphere::dist2Line to get the distance from your points to the coastline.

我注意到您的示例点之一是在法国，因此您可能需要将海岸线数据扩展到仅英国以外(可以通过使用边界框的范围来完成).

I noticed that one of your example points was in France so you may need to expand the coastline data beyond just the UK (can be done by playing with the extents of the bounding box).

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-多边形，然后查看您的每个点是否都位于海中.

That's just for the distance to the coastline. You also need to know how whether it is inland or at sea. For this, you would need a separate shapefile for the sea: http://openstreetmapdata.com/data/water-polygons and see if each point of your points sits in the sea or not.

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

由于您要求一种专门使用shapefile的方法，因此我在这里下载了该方法:openstreetmapdata.com/data/coastlines，我将使用它来执行与上述相同的方法.

Since you asked for an approach specifically using a shapefile I have downloaded this one here: openstreetmapdata.com/data/coastlines which I will use to carry out the same approach as above.

clines <- read_sf('lines.shp') #path to shapefile

接下来，我创建了一个自定义边界框，以便我们可以减少shapefile的大小，使其仅包括合理接近点的海岸线.

Next I created a custom bounding box so that we can cut down the size of the shapefile to only include coastlines reasonably close to the points.

# 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函数，因为它很准确，并且可以为您提供海岸线上要测量的点，这对于检查错误非常有用.不利的一面是，对于我们相当大的海岸线文件来说，它非常慢.

Now I used the dist2Line function here because it is accurate and it gives you the points on the coastline it is measuring to which makes it good for checking errors. The downside is, it's very slow for our rather large coastline file.

运行此过程花了我8分钟:

Running this took me 8 minutes:

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％)，并且对知道要测量的海岸线确切位置不感兴趣，则可以测量到多边形:

If you don't mind about the slight loss of accuracy (results differ by around 0.3% on your data), and are not fussed about knowing where exactly on the coastline it is measuring to, you can measure the distance to the polygon:

# 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秒钟即可运行！

By contrast this took just 8 seconds to run!

其余与上一个答案相同.

The rest is as in the previous answer.

这篇关于R测量距海岸线的距离的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！