Блог

Я импортирую CSV-файл с координатами долготы и широты и преобразую их в шейп-файлы полигонов. Я размещаю сетку над полигонами и нахожу центроид каждого квадрата сетки. Затем я извлекаю координаты центроидов и помещаю их во фрейм данных, но мне нужно иметь возможность сказать, в каком полигоне находится конкретный центроид.

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

Решение этой проблемы состоит в том, чтобы выполнить функцию point-in-polygon через st_join .

С tidyverse это довольно просто, и я уверен, что вы можете использовать следующее для перевода в базу R.

(Я взял на себя смелость немного изменить ваши воспроизводимые данные, поскольку polygon 4 это недопустимый полигон, учитывая, что он имеет только 3 точки):

Сначала мы создаем sf из фрейма данных xy

 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 введете выходные данные, вы получите