提取多边形质心的坐标并用多边形编号标记它们

Question

我导入了经度和纬度坐标的 csv，并将它们转换为多边形 shapefile。我在多边形上放置一个网格并找到每个网格正方形的质心。然后我提取质心的坐标并将其放在数据框中，但我需要能够说出特定质心所在的多边形。

#Create shapefile of polygons
polygon <- lapply(split(df, df$shape), function(x) { coords <- 
as.matrix(cbind(x$longitude, x$latitude)); list(rbind(coords, coords[1,]))}) 
Coord_Ref <- st_crs(4326)
polygon <-  st_sfc(st_multipolygon(x=polygon), crs = Coord_Ref)
polygon <-  st_cast(polygon, "POLYGON")

#Create grid and centroids
PolygonBits <- st_make_grid(polygon, cellsize=0.0002)
PolygonBitCentroids <- st_centroid(st_make_grid(polygon, cellsize=0.0002))

#Extract coordinates and place them in dataframe
PolygonBitCentroids <- st_coordinates(PolygonBitCentroids)
PolygonBitCentroids <- as.data.frame(PolygonBitCentroids)

PolygonBitCentroids 数据框的前三行如下所示：

         X      Y
1   -0.0014 0.1990
2   -0.0012 0.1990
3   -0.0010 0.1990

但我需要这样的东西：

          X      Y  Shape
1   -0.0014 0.1990  Polygon 1
2   -0.0012 0.1990  Polygon 1
3   -0.0010 0.1990  Polygon 1

可重现的数据：

structure(list(shape = c("polygon 1", "polygon 1", "polygon 1", 
"polygon 1", "polygon 2", "polygon 2", "polygon 2", "polygon 2", 
"polygon 3", "polygon 3", "polygon 3", "polygon 3", "polygon 4", 
"polygon 4", "polygon 4", "polygon 4"), longitude = c(0, 1, 1, 
0, 1.5, 2, 2, 1.5, -2, -2, -1, -1, 0, 1, 1, 0), latitude = c(1, 
1, 0, 0, 1, 1, 0, 0, -0.5, -2, -2, -0.5, 1.5, 1.5, 2, 2)), class = 
"data.frame", row.names = c(NA, 
-16L), spec = structure(list(cols = list(shape = structure(list(), 
class = c("collector_character", 
"collector")), longitude = structure(list(), class = 
c("collector_double", 
"collector")), latitude = structure(list(), class = 
c("collector_double", 
"collector"))), default = structure(list(), class = 
c("collector_guess", 
"collector")), skip = 1), class = "col_spec"))

Answer 1

这个问题的解决方案是 point-in-polygon 通过 st_join。

这对于 tidyverse 来说非常简单，我相信您可以使用以下内容转换为基础 R。

~~（我冒昧地稍微更改了您的可重现数据，因为 polygon 4 不是有效的多边形，因为它只有 3 个点）：~~

首先我们从xy数据帧

创建一个sf

library(sf)
library(tidyverse)

polygons <- polygons %>%
  st_as_sf(coords = c('longitude', 'latitude')) %>%
  st_set_crs(4326)

绘制时，它看起来像这样

polygons <- polygons %>%
  group_by(shape) %>%
  summarise(do_union=FALSE) %>%
  st_cast("POLYGON")

这正确地从点创建了多边形。

对 plot(polygons) 的调用会生成以下图：

（do_union=FALSE 参数很重要，否则不会保留顺序）。接下来我们为网格创建一个单独的 sf 对象：

grids <- polygons %>%
  st_make_grid(cellsize = 0.2) %>%
  st_centroid() %>%
  st_sf()

最后，我们加入两个sf objects using st_within`

grids %>% st_join(polygons, join = st_within)

您得到的与您要求的完全一样。

Simple feature collection with 92 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: -1.9 ymin: -1.9 xmax: 1.9 ymax: 1.9
CRS:            EPSG:4326
First 10 features:
       shape          geometry
1       <NA> POINT (-1.9 -1.9)
2       <NA> POINT (-1.1 -1.9)
3       <NA> POINT (-0.9 -1.9)
4  polygon 3 POINT (-1.9 -1.7)
5       <NA> POINT (-1.7 -1.7)
6       <NA> POINT (-1.3 -1.7)
7  polygon 3 POINT (-1.1 -1.7)
8       <NA> POINT (-0.9 -1.7)
9  polygon 3 POINT (-1.9 -1.5)
10 polygon 3 POINT (-1.7 -1.5)

如果你plot输出，你会得到

提取多边形质心的坐标并用多边形编号标记它们

Extract coordinates of polygon centroids and label them by polygon number

r

shapefile

shapes

sf