创建网格单元和网格单元ID

Question

我有带有计数的纬度和经度的数据框

我想在瑞典数据中创建 10 x10 km 的网格单元，并生成网格单元 ID 以及该单元内的采样点数量。

下面的代码行似乎没有解决问题

瑞典_MM4

library(sf)
# Convert Sweden_MM4  to an sf object
my_data_sf <- st_as_sf(Sweden_MM4, coords = c("longitude", "latitude"), crs = 4326)
my_data_sf <- st_transform(my_data_sf, crs = "+proj=utm +zone=33 +datum=WGS84")  # Set the CRS to an appropriate projection that uses meters 
# Create a grid of 10x10 km using st_make_grid
grid_size_km <- 10
grid_size_m <- grid_size_km * 1000
grid_polygons <- st_make_grid(my_data_sf, cellsize = grid_size_m, what = "polygons")
# Set the CRS of the grid to match your data
grid_polygons <- st_set_crs(grid_polygons, st_crs(my_data_sf))
# Perform a spatial overlay to associate each data point with the corresponding grid cell ID
grid_overlay <- st_intersection(grid_polygons, my_data_sf)
str(my_data_sf)

预期输出

Answer 1

您还需要为网格单元生成ID并构建一个实际的

sf

对象，

st_make_grid()

的返回值只是一个几何列表列。为了将网格 ID 分配给您的点，您可能需要使用

st_join()

，默认情况下它使用

st_intersects

谓词。

另请注意

st_intersects

（二元谓词）和

st_intersection

（几何运算）之间的区别。

准备代表

library(sf)
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(dplyr)
library(ggplot2)

swe <- giscoR::gisco_get_countries(country = "Sweden")

# sample points
set.seed(123)
my_data_sf <- st_geometry(swe) %>% 
  st_sample(20) %>% 
  st_as_sf() %>% 
  st_set_geometry("geometry") %>% 
  mutate(random_point_id = row_number(), .before = 1) %>% 
  st_transform("+proj=utm +zone=33 +datum=WGS84")
print(my_data_sf, n = 3)
#> Simple feature collection with 20 features and 1 field
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 291587.2 ymin: 6421365 xmax: 850403.5 ymax: 7480341
#> Projected CRS: +proj=utm +zone=33 +datum=WGS84
#> First 3 features:
#>   random_point_id                 geometry
#> 1               1 POINT (493032.3 6437133)
#> 2               2   POINT (585428 6482855)
#> 3               3 POINT (836136.5 7424094)

从网格生成

sf

对象，添加

grid_id

列

# grid, 200x200km for better visualization
grid_size_km <- 200
grid_size_m <- grid_size_km * 1000
grid_polygons_sfc <- st_make_grid(my_data_sf, cellsize = grid_size_m, what = "polygons")

# st_make_grid() does not return sf object but a geometry column:
print(grid_polygons_sfc, n = 3)
#> Geometry set for 18 features 
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 291587.2 ymin: 6421365 xmax: 891587.2 ymax: 7621365
#> Projected CRS: +proj=utm +zone=33 +datum=WGS84
#> First 3 geometries:
#> POLYGON ((291587.2 6421365, 491587.2 6421365, 4...
#> POLYGON ((491587.2 6421365, 691587.2 6421365, 6...
#> POLYGON ((691587.2 6421365, 891587.2 6421365, 8...

# create sf object, add grid_id attribute column:
grid_polygons_sf <- st_sf(grid_id = seq_along(grid_polygons_sfc), 
                          geometry = grid_polygons_sfc)
print(grid_polygons_sf, n = 3)
#> Simple feature collection with 18 features and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 291587.2 ymin: 6421365 xmax: 891587.2 ymax: 7621365
#> Projected CRS: +proj=utm +zone=33 +datum=WGS84
#> First 3 features:
#>   grid_id                       geometry
#> 1       1 POLYGON ((291587.2 6421365,...
#> 2       2 POLYGON ((491587.2 6421365,...
#> 3       3 POLYGON ((691587.2 6421365,...

加入数据集

# spatial join, join grid_polygons_sf$grid_id values to my_data_sf
my_data_sf <- st_join(my_data_sf, grid_polygons_sf)
my_data_sf
#> Simple feature collection with 20 features and 2 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 291587.2 ymin: 6421365 xmax: 850403.5 ymax: 7480341
#> Projected CRS: +proj=utm +zone=33 +datum=WGS84
#> First 10 features:
#>    random_point_id grid_id                 geometry
#> 1                1       2 POINT (493032.3 6437133)
#> 2                2       2   POINT (585428 6482855)
#> 3                3      18 POINT (836136.5 7424094)
#> 4                4      14 POINT (637972.7 7301930)
#> 5                5      14 POINT (595969.3 7221551)
#> 6                6      18 POINT (713168.4 7480341)
#> 7                7       6 POINT (701661.1 6638757)
#> 8                8      10 POINT (375548.9 7076110)
#> 9                9       4 POINT (463189.1 6647020)
#> 10              10       1 POINT (307018.8 6498219)

ggplot() +
  geom_sf(data = swe) +
  geom_sf(data = grid_polygons_sf, alpha = .2) +
  geom_sf_text(data = grid_polygons_sf, aes(label = grid_id)) +
  geom_sf(data = my_data_sf, alpha = .8, color = "gold") +
  geom_sf_label(data = my_data_sf, aes(label = grid_id), size = 3, alpha = .5) +
  coord_sf(crs = "+proj=utm +zone=33")

^{创建于 2023-07-28，使用 reprex v2.0.2}