How to find centroids of multi-polygons using sf library

Question

I have a shapefile of multipolygons and I want to find the centroid of each multipolygon, but instead I am getting the centroid of the whole shapefile.

I converted my csv file into a shapefile (see this question Converting CSV file to shapefile - but want polygon not points) using the following lines of code:

df <- as.data.frame(read_csv("/Users/alan/Desktop/shapes.csv"))
df1 <- lapply(split(df, df$shape), function(x) { coords <- as.matrix(cbind(x$longitude, 
x$latitude)); rbind(coords, coords[1,])}) 
Coord_Ref <- st_crs(3035)
plot_locations_df <-  st_multipolygon(x=list(df1))
plot_locations_df <-  st_sfc(plot_locations_df, crs = Coord_Ref)
plot(st_geometry(plot_locations_df))
plot(st_geometry(plot_locations_df, col = sf.colors(12, categorical=TRUE), border='grey', 
axes=TRUE))
plot(st_geometry(st_centroid(plot_locations_df)), pch = 3, col = 'red', add=TRUE)

My dataframe looks like this:

structure(list(shape = c(1.1, 1.1, 1.1, 1.1, 2.1, 2.1, 2.1, 2.1, 
3.1, 3.1, 3.1, 3.1, 4.1, 4.1, 4.1, 4.1, 5.1, 5.1, 5.1, 5.1, 6.1, 
6.1, 6.1, 6.1, 7.1, 7.1, 7.1, 7.1, 8.1, 8.1, 8.1, 8.1, 9.1, 9.1, 
9.1, 9.1), longitude = c(43, 43, 40, 40, 23, 23, 20, 20, 25, 
25, 38, 38, 25, 25, 38, 38, 45, 50, 50, 45, 65, 60, 60, 65, 60, 
60, 80, 80, 60, 60, 80, 80, 20, 20, 80, 80), latitude = c(10, 
13, 13, 10, 10, 13, 13, 10, 10, 10.3, 10.3, 10, 12.7, 13, 13, 
12.7, 13, 13, 10, 10, 13, 13, 10, 10, 9.8, 9.5, 9.5, 9.8, 6, 
5.7, 5.7, 6, 5, 4.5, 4.5, 5)), row.names = c(NA, 36L), class = "data.frame")

Answer 1

You can create your polygons directly from your .csv using sfheaders

library(sfheaders)
library(sf)

sf <- sfheaders::sf_polygon(
  obj = df
  , polygon_id = "shape"
  , x = "longitude"
  , y = "latitude"
)

## Then you can get the centroids of each polygon
sf_centers <- sf::st_centroid( sf )

Answer 2

you have two different problems, here. First of all, the method you are using to create the sf object leads to an invalid geometry (st_is_valid(plot_locations_df) returns FALSE). To get a valid MULTIPOLYGON you can use:

df1 <- lapply(split(df, df$shape), function(x) { coords <- as.matrix(cbind(x$longitude, 
x$latitude)); list(rbind(coords, coords[1,]))}) 
names(df1)<- NULL
Coord_Ref <- st_crs(3035)
plot_locations_df <-  st_sfc(st_multipolygon(x=df1), crs = Coord_Ref)

st_is_valid(plot_locations_df)
[1] TRUE

However, this still does not help, because your geometry is still a MULTIPOLYGON (i.e., a single feature made of multiple POLYGONS), and the centroid of a MULTIPOLYGON is a single point, computed taking into account all of its polygons.

plot_locations_df
> Geometry set for 1 feature 
> geometry type:  MULTIPOLYGON
> dimension:      XY
> bbox:           xmin: 20 ymin: 4.5 xmax: 80 ymax: 13
> projected CRS:  ETRS89-extended / LAEA Europe
> MULTIPOLYGON (((43 10, 43 13, 40 13, 40 10, 43 ...

st_centroid(plot_locations_df)
> Geometry set for 1 feature 
> geometry type:  POINT
> dimension:      XY
> bbox:           xmin: 49.10736 ymin: 8.969325 xmax: 49.10736 ymax: 8.969325
> projected CRS:  ETRS89-extended / LAEA Europe
> POINT (49.10736 8.969325)

To get what you want, you will have to split the polygons, by recasting to POLYGON:

plot_locations_df <- st_cast(plot_locations_df, "POLYGON")
plot_locations_df

> Geometry set for 9 features 
> geometry type:  POLYGON
> dimension:      XY
> bbox:           xmin: 20 ymin: 4.5 xmax: 80 ymax: 13
> projected CRS:  ETRS89-extended / LAEA Europe
> First 5 geometries:
> POLYGON ((43 10, 43 13, 40 13, 40 10, 43 10))
> POLYGON ((23 10, 23 13, 20 13, 20 10, 23 10))
> POLYGON ((25 10, 25 10.3, 38 10.3, 38 10, 25 10))
> POLYGON ((25 12.7, 25 13, 38 13, 38 12.7, 25 12...
> POLYGON ((45 13, 50 13, 50 10, 45 10, 45 13))

st_centroid(plot_locations_df)
> Geometry set for 9 features 
> geometry type:  POINT
> dimension:      XY
> bbox:           xmin: 21.5 ymin: 4.75 xmax: 70 ymax: 12.85
> projected CRS:  ETRS89-extended / LAEA Europe
> First 5 geometries:
> POINT (41.5 11.5)
> POINT (21.5 11.5)
> POINT (31.5 10.15)
> POINT (31.5 12.85)
> POINT (47.5 11.5)

plot(st_geometry(plot_locations_df))
plot(st_geometry(plot_locations_df, col = sf.colors(12, categorical=TRUE), border='grey', 
axes=TRUE))
plot(st_geometry(st_centroid(plot_locations_df)), pch = 3, col = 'red', add=TRUE)

HTH!