在R语言中生成Shapefile指定区域随机点并判断归属

数据加载 #

先加载所需包并读取北卡罗来纳州地理数据：

library(sf)  
library(leaflet)
library(leafgl)
library(colourvalues)
library(leaflet.extras)

nc <- st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet = TRUE) %>% 
  st_transform(st_crs(4326)) %>% 
  st_cast('POLYGON')

数据结构如下：

> nc
Simple feature collection with 108 features and 14 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -84.32377 ymin: 33.88212 xmax: -75.45662 ymax: 36.58973
Geodetic CRS:  WGS 84
First 10 features:
     AREA PERIMETER CNTY_ CNTY_ID        NAME  FIPS FIPSNO CRESS_ID BIR74 SID74 NWBIR74 BIR79 SID79 NWBIR79                           geom
1   0.114     1.442  1825    1825        Ashe 37009  37009        5  1091     1      10  1364     0      19 POLYGON ((-81.47258 36.2344...
2   0.061     1.231  1827    1827   Alleghany 37005  37005        3   487     0      10   542     3      12 POLYGON ((-81.23971 36.3654...
3   0.143     1.630  1828    1828       Surry 37171  37171       86  3188     5     208  3616     6     260 POLYGON ((-80.45614 36.2426...
4   0.070     2.968  1831    1831   Currituck 37053  37053       27   508     1     123   830     2     145 POLYGON ((-76.00863 36.3196...
4.1 0.070     2.968  1831    1831   Currituck 37053  37053       27   508     1     123   830     2     145 POLYGON ((-76.02682 36.5567...
4.2 0.070     2.968  1831    1831   Currituck 37053  37053       27   508     1     123   830     2     145 POLYGON ((-75.90164 36.5562...
5   0.153     2.206  1832    1832 Northampton 37131  37131       66  1421     9    1066  1606     3    1197 POLYGON ((-77.21736 36.2410...
6   0.097     1.670  1833    1833    Hertford 37091  37091       46  1452     7     954  1838     5    1237 POLYGON ((-76.74474 36.2339...
7   0.062     1.547  1834    1834      Camden 37029  37029       15   286     0     115   350     2     139 POLYGON ((-76.00863 36.3196...
8   0.091     1.284  1835    1835       Gates 37073  37073       37   420     0     254   594     2     371 POLYGON ((-76.56218 36.3406...

目标1：生成Ashe区域内的随机点 #

思路1：先生成大范围随机点，再筛选Ashe区域内的点 #

先生成目标范围附近的随机点，转换为sf空间对象后筛选出落在Ashe区域内的点：

set.seed(123) # 设置随机种子保证结果可复现
id <- 1:100
longitude <- rnorm(100, -81, 0.15) # 注意：空间坐标顺序为 经度(x), 纬度(y)
latitude <- rnorm(100, 36.2, 0.15)
my_data <- data.frame(id, longitude, latitude)

# 转换为sf点对象
my_data_sf <- st_as_sf(my_data, coords = c("longitude", "latitude"), crs = st_crs(nc))

# 提取Ashe区域
ashe_area <- nc[nc$NAME == "Ashe", ]

# 筛选出落在Ashe内的点
ashe_points <- my_data_sf[ashe_area, ]

# 查看结果
head(ashe_points)

思路2：直接在Ashe区域内生成随机点 #

利用st_sample函数直接从Ashe多边形中采样生成随机点，无需后续筛选：

set.seed(123)
# 提取Ashe区域
ashe_area <- nc[nc$NAME == "Ashe", ]
# 生成100个随机点
ashe_points_direct <- st_sample(ashe_area, size = 100, type = "random")
# 转换为带id和经纬度的格式
ashe_points_direct <- st_sf(id = 1:100, geometry = ashe_points_direct)
ashe_points_direct$longitude <- st_coordinates(ashe_points_direct)[,1]
ashe_points_direct$latitude <- st_coordinates(ashe_points_direct)[,2]

# 查看结果
head(ashe_points_direct)

目标2：判断随机点所属区域 #

对于已生成的随机点（如思路1中的my_data_sf），使用st_join进行空间连接，获取每个点所属的区域名称：

# 空间连接：将点与nc区域匹配，仅保留点在多边形内的匹配结果
point_with_location <- st_join(my_data_sf, nc[, "NAME"], join = st_within)

# 转换为普通数据框查看（可选）
point_with_location_df <- as.data.frame(point_with_location) %>% 
  select(-geometry)

# 查看结果
head(point_with_location_df)

输出示例：

id  longitude latitude  NAME
1  1 -81.15816  36.42389  Ashe
2  2 -81.40090  36.23823  Ashe
3  3 -80.97732  35.97633 Surry
4  4 -80.80150  36.20300  <NA> # 不在任何区域内的点会显示NA
5  5 -81.26429  36.23899  Ashe
6  6 -81.13721  36.31100 Surry

内容的提问来源于stack exchange，提问作者stats_noob