识别靠近多边形边界的点

Question

我试图识别位于英国海岸线附近(即多边形)的所有点(邮政编码)。我用R来处理这件事。

我从这里下载了英国的地理轮廓作为一个文件。英国的所有邮政编码列表都是从英国国家统计局( ONS 这里 )访问的。请注意，后一个文件非常大(211 is压缩)。

首先，我将两个文件加载到R中，然后将它们转换为相同的坐标参考系统(OSGB1936；27700)。对于联合王国的多边形，我将其转换为代表边界/海岸线的线(注意，虽然北爱尔兰与爱尔兰有共同的边界，但我将把任何误配为lat/long与海岸线附近的邮政编码进行子集)。然后我将这些点转换成空间点。

# Load libraries
library(sf)
library(data.table)

# Load data
uk_shp <- read_sf("./GBR_adm/GBR_adm0.shp") # Load UK shapefile (ignore the download file says GBR, it is UK)
uk_shp <- st_transform(uk_shp, crs = 27700) # Convert to co-ordinate reference system (CRS) that allow buffers in correct units later (note: 4326 is World CRS)
uk_coast <- st_cast(uk_shp,"MULTILINESTRING") # Convert polygon to a line (i.e., coastline)

# Load in postcodes
pcd <- fread("./ONSPD_FEB_2022_UK/Data/ONSPD_FEB_2022_UK.csv") # Load all postcodes for Great Britain - this is a very large file so I also create a single 
pcd <- pcd[, c(1:3, 43:44)] # Drop unnecessary information/columns to save memory

# Convert to spatial points data frame
pcd_sp <- pcd %>% # For object of postcodes
  st_as_sf(coords = c("long", "lat")) %>% # Define as spatial object and identify which columns tell us the position of points
  st_set_crs(27700) # Set CRS

我最初认为最有效的方法是定义沿海区域(这里定义为海岸线5公里以内)，创建一个缓冲区来表示海岸线周围的海岸线，然后使用点多边形函数来选择缓冲区中的所有点。然而，下面的代码并没有在一夜之间运行完，这可能表明这是不正确的方法，我不知道为什么要花这么长时间。

uk_coast <- st_buffer(uk_coast, 5000) # Create 5km buffer
pcd_coastal <- st_intersection(uk_buf, pcd_sp) # Point-in-polygon (i.e., keep only the postcodes that are located in the buffer region)

所以我改变了我的方法来计算每个点到最近的海岸线的直线距离。在运行下面的代码时，它提供了不正确的距离。例如，下面我选择了一个邮政编码(AB12 4XP)，它位于距海岸线2.6公里的地方，但是下面的代码给出了82公里的高度，这是非常错误的。我曾经尝试过st_nearest_feature()，但无法让它工作(它可能会起作用，但超出了我的尝试范围)。

test <- pcd_sp[pcd_sp$pcd == "AB124XP",] # Subset test postcode
dist <- st_distance(test, uk_coast, by_element = TRUE, which = "Euclidean") # Calculate distance

我不知道如何从这里开始-我不认为这是错误的CRS。可能是多行字符串转换导致了问题。有人有什么建议吗？

Answer 1

sf有一个st_is_within_distance函数，可以测试点是否在一条线的距离之内。我的测试数据是英国形状边界框中的10,000个随机点，以及OSGB网格坐标中的UK形状。

> system.time({indist = st_is_within_distance(uk_coast, pts, dist=5000)})
   user  system elapsed 
 30.907   0.003  30.928 

但这并不是在建立一个空间索引。docs说，如果坐标是“地理”的，并且设置了使用球形几何学的标志，它就会构建一个空间索引。我不明白为什么它不能为笛卡尔坐标建立一个，但让我们看看它有多快.

转换根本不需要时间：

> ukLL = st_transform(uk_coast, 4326)
> ptsLL = st_transform(pts, 4326)

那就测试..。

 system.time({indistLL = st_is_within_distance(ukLL, ptsLL, dist=5000)})
   user  system elapsed 
  1.405   0.000   1.404 

稍过一秒。两者之间有什么区别吗？让我们看看：

> setdiff(indistLL[[1]], indist[[1]])
[1] 3123
> setdiff(indist[[1]], indistLL[[1]])
integer(0)

因此，第3123点使用的是lat-长，而不是使用OSGB的集合。在OSGB中，没有什么不是长时间设置的。

显示所选点的快速绘图：

> plot(uk_coast$geometry)
> plot(pts$geometry[indistLL[[1]]], add=TRUE)

​

​