Extract coordinates of polygon centroids and label them by polygon number

Question

I import a csv of longitude and latitude coordinates and convert them to polygon shapefiles. I place a grid over the polygons and find the centroid of each grid square. I then extract the coordinates of the centroids and place it in a dataframe, but I need to be able to say which polygon a particular centroid is in.

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

The first three rows of the PolygonBitCentroids dataframe looks as follows:

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

But I need something like this:

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

Reproducible data:

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

yeedle yeedle · Accepted Answer · 2020-08-24T20:22:12

The solution to this problem is to do point-in-polygon via st_join.

This is pretty straightforward with the tidyverse, and I'm sure you can use the following to translate to base R.

~~(I took the liberty to change your reproducible data slightly, since polygon 4 is not a valid polygon given that it only has 3 points):~~

First we create an sf from the xy dataframe

library(sf)
library(tidyverse)

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

When plotted, this looks like this

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

This correctly creates the polygons from the points.

a call to plot(polygons) produces the following plot:

(the do_union=FALSE argument is important because otherwise order is not preserved). Next we create a separate sf object for the grids:

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

Finally, we join the two sf objects using st_within`

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

What you get looks exactly as you asked for.

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)

If you plot the output, you'll get

Extract coordinates of polygon centroids and label them by polygon number

1 Answers