how to get the 'mean' point of the polygon's coordinates?

I have created a data frame (projects) by merging a shape file (US Zip codes) and a dataset. The below line of code returns the X and Y coordinates of my data frame. I would like to get the 'mean' point of the polygon's coordinates (lat and long) and add them to the projects data frame. I tried st_centroid() but got an error:

Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented, : Loop 0 is not valid: Edge 1616 has duplicate vertex with edge 2124 –

#> head(st_coordinates(projects))
             X        Y L1 L2 L3
[1,] -71.36374 41.85854  1  1  1
[2,] -71.36449 41.85847  1  1  1
[3,] -71.36473 41.85844  1  1  1
[4,] -71.36580 41.85834  1  1  1
[5,] -71.36573 41.85828  1  1  1
[6,] -71.36562 41.85818  1  1  1

Is this about a single centroid for your whole dataset or finding a centroid for each (multi)polygon? Provided data seems to hint former as it lists just coordinate pairs, though question states polygons.

Using dataset from sf package. First example starts with sf object of multipolygons, takes an union and from that finds a single centroid:

library(sf)
nc = st_read(system.file("shape/nc.shp", package="sf"))

### centroid of multiploygons
head(nc[c("NAME", "geometry")])
#> Simple feature collection with 6 features and 1 field
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -81.74107 ymin: 36.07282 xmax: -75.77316 ymax: 36.58965
#> Geodetic CRS:  NAD27
#>          NAME                       geometry
#> 1        Ashe MULTIPOLYGON (((-81.47276 3...
#> 2   Alleghany MULTIPOLYGON (((-81.23989 3...
#> 3       Surry MULTIPOLYGON (((-80.45634 3...
#> 4   Currituck MULTIPOLYGON (((-76.00897 3...
#> 5 Northampton MULTIPOLYGON (((-77.21767 3...
#> 6    Hertford MULTIPOLYGON (((-76.74506 3...
c_multiploy <- st_union(nc) |> st_centroid()
c_multiploy
#> Geometry set for 1 feature 
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -79.40065 ymin: 35.55937 xmax: -79.40065 ymax: 35.55937
#> Geodetic CRS:  NAD27
#> POINT (-79.40065 35.55937)

# multiploygons (grey):
plot(nc$geometry, lwd = .15, border = "grey70", main = paste0("c_multiploy = ", c_multiploy))
# union of multiploygons (blue):
plot(st_union(nc), lwd = 2, border = "blue", add =TRUE)
# centroid of union (red):
plot(st_union(nc) |> st_centroid(), pch = 4, cex = 3, lwd = 5, col = "red", add =TRUE)

2nd example first generates centroids of mulitpolygons, so if this is what you are after, you can just stop there. Then finds a union (multipoint) to get a single centroid.

### centroid of points:
nc_points <- st_centroid(nc)
head(nc_points[c("NAME", "geometry")])
#> Simple feature collection with 6 features and 1 field
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -81.49823 ymin: 36.36142 xmax: -76.02719 ymax: 36.49111
#> Geodetic CRS:  NAD27
#>          NAME                   geometry
#> 1        Ashe  POINT (-81.49823 36.4314)
#> 2   Alleghany POINT (-81.12513 36.49111)
#> 3       Surry POINT (-80.68573 36.41252)
#> 4   Currituck POINT (-76.02719 36.40714)
#> 5 Northampton POINT (-77.41046 36.42236)
#> 6    Hertford POINT (-76.99472 36.36142)
c_points <- st_union(nc_points) |> st_centroid()
c_points
#> Geometry set for 1 feature 
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -79.5111 ymin: 35.64247 xmax: -79.5111 ymax: 35.64247
#> Geodetic CRS:  NAD27
#> POINT (-79.5111 35.64247)

# points (grey):
plot(nc_points$geometry, pch = 1, cex = 3, col = "grey70", main = paste0("c_points = ", c_points))
# union of points (blue):
plot(st_union(nc_points), pch = 3, cex = 1, col = "blue", add = TRUE)
# centroid of union (red):
plot(st_union(nc_points) |> st_centroid(), pch = 4, cex = 3, lwd = 5, add =TRUE, col = "red")

^{Created on 2023-01-21 with reprex v2.0.2}