Convert geometry column and retain other columns

I have imported a shapefile (.shp) into R to define a set of polygons. The shapefile takes the form of a grouped dataframe with ageometry column containing a MULTIPOLYGON list of XY coordinates. I have added some extra columns of data to this dataframe:

ID   Value   geometry
A    5       MULTIPOLYGON (((519270 1639...
B    10      MULTIPOLYGON (((519553 1642...
C    25      MULTIPOLYGON (((519000 1666...

This works fine for generating static choropleth maps using ggplot and geom_sf. However, to generate dynamic maps using leaflet I need to convert the XY coordinates that define the polygons into lat/long coordinates.

I have tried using st_transform to do this and, while it does generate a set of polygons, I cannot colour those polygons as there are no values associated with them:

inter_map2 <- leaflet(df) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(data = st_transform(df$geometry, 4326), stroke = TRUE, weight = 0.5, opacity = 1, fillOpacity = 0.05, highlightOptions = highlightOptions(color = "red", weight = 2, bringToFront = TRUE))

I have also tried using extract but leaves the resulting lat and long columns with NA values:

df_latlong <- df %>%
  extract(geometry, c('lat', 'lon'), '\\((.*), (.*)\\)', convert = TRUE)

Is there a way to do this while also retaining the ID and Value columns associated with the polygons?

Solution

sf object is bit more than just a regular tibble / dataframe with coordinates in a column named geometry, attributes and geometries are meant to be used together and in most cases there's no reason to split those by force.

extract() will not work as geometry is not really a string, printing it just outputs geometries as a WKT (well-known-text). And by calling st_transform(df$geometry, 4326) you are indeed left with just the transformed geometry list, attribute columns get dropped; it has its uses but its more common to transform the sf objects and not just the geometry column, i.e. do st_transform(df, 4326).

data parameter in leaflet() sets the default data object, in your case it's df. In following addPolygons() calls you don't need to override it unless you want to add more layers, just use formula notation (~) and reference df columns. Something like this:

library(leaflet)
library(sf)

# sf exampl shapefile, tranform to projected CRS
nc <- st_read(system.file("shape/nc.shp", package="sf")) %>% st_transform(6542)
nc[,1:2]
#> Simple feature collection with 100 features and 2 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 123829.8 ymin: 14740.06 xmax: 930518.6 ymax: 318255.5
#> Projected CRS: NAD83(2011) / North Carolina
#> First 10 features:
#>     AREA PERIMETER                       geometry
#> 1  0.114     1.442 MULTIPOLYGON (((387344.7 27...
#> 2  0.061     1.231 MULTIPOLYGON (((408601.4 29...
#> 3  0.143     1.630 MULTIPOLYGON (((478715.7 27...
#> 4  0.070     2.968 MULTIPOLYGON (((878193.4 28...
#> 5  0.153     2.206 MULTIPOLYGON (((769834.9 27...
#> 6  0.097     1.670 MULTIPOLYGON (((812327.7 27...
#> 7  0.062     1.547 MULTIPOLYGON (((878193.4 28...
#> 8  0.091     1.284 MULTIPOLYGON (((828444.5 29...
#> 9  0.118     1.421 MULTIPOLYGON (((671746.3 27...
#> 10 0.124     1.428 MULTIPOLYGON (((517435.1 27...

# pallete for leaflet from nc$AREA values
pal <- colorNumeric(
  palette = "Greens",
  domain = nc$AREA
)

# pass transformed sf object directly to leaflet and in 
# addPolygons() reference its columns though formula (~) interface:
nc %>% 
  st_transform("WGS84") %>% 
  leaflet() %>%
  addTiles() %>% 
  addPolygons(color = ~pal(AREA), stroke = FALSE, fillOpacity = 1)

^{Created on 2023-06-14 with reprex v2.0.2}