## multipolygon
library(sf)
## Linking to GEOS 3.5.1, GDAL 2.2.1, proj.4 4.9.3
p1 <- rbind(c(0,0), c(1,0), c(3,2), c(2,4), c(1,4), c(0,0))
p2 <- rbind(c(1,1), c(1,2), c(2,2), c(1,1))
pol <-st_polygon(list(p1,p2))
p3 <- rbind(c(3,0), c(4,0), c(4,1), c(3,1), c(3,0))
p4 <- rbind(c(3.3,0.3), c(3.8,0.3), c(3.8,0.8), c(3.3,0.8), c(3.3,0.3))[5:1,]
p5 <- rbind(c(3,3), c(4,2), c(4,3), c(3,3))
(mpol <- st_multipolygon(list(list(p1,p2), list(p3,p4), list(p5))))
## MULTIPOLYGON(((0 0, 1 0, 3 2, 2 4, 1 4, 0 0), (1 1, 1 2, 2 2, 1 1)), ((3 0, 4 0, 4 1, 3 1, 3 0), (3.3 0.3, 3.3 0.8, 3.8 0.8, 3.8 0.3, 3.3 0.3)), ((3 3, 4 2, 4 3, 3 3)))
mp <- st_sfc(mpol, crs = "+proj=merc +datum=WGS84")
library(leaflet)

## 1. standard sf example
leaflet() %>% leaflet::addPolygons(data = mp)
## Warning: sf layer is not long-lat data
library(leaflet)
library(sf)
## a couple of realer examples 850KB and 1.4Mb
uf <- c("https://github.com/r-gris/polyggon/raw/master/inst/extdata/water_lake_superior_basin.gpkg",
        "https://github.com/r-gris/polyggon/raw/master/inst/extdata/inlandwaters.gpkg")
td <- tempdir()
for (i in seq_along(uf)) download.file(uf[i], file.path(td, basename(uf[i])), mode = "wb")

## read the layers from the GeoPackage (it's a SQLite db, understood by GDAL)
lake_superior <- st_read(file.path(td, basename(uf[1])), layer = "lake_superior_basin")
## Reading layer `lake_superior_basin' from data source `/tmp/Rtmpqbop96/water_lake_superior_basin.gpkg' using driver `GPKG'
## Simple feature collection with 525 features and 5 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 565501.9 ymin: 5155179 xmax: 1165764 ymax: 5441497
## epsg (SRID):    26915
## proj4string:    +proj=utm +zone=15 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
lakey <- st_read(file.path(td, basename(uf[2])), layer = "inlandwaters")
## Reading layer `inlandwaters' from data source `/tmp/Rtmpqbop96/inlandwaters.gpkg' using driver `GPKG'
## Simple feature collection with 9 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -2239268 ymin: -2677570 xmax: 2477934 ymax: 2429213
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=-47 +lat_2=-17 +lat_0=-32 +lon_0=136 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## reproject to longlat :|
lake_superior_ll <- st_transform(lake_superior, "+proj=longlat +datum=WGS84")
lakey_ll <- st_transform(lakey, "+proj=longlat +datum=WGS84")

## 2. Lake Superior
##
## the holes have other features in them :)
leaflet() %>%  addPolygons(data = lake_superior_ll)
## 3. Lake Superior
## this is the main polygon only, which has holes in it
idx <- which.max(st_area(lake_superior))
leaflet() %>%  addPolygons(data = lake_superior_ll[idx, ])
## Inland waters, Australia
##
## these multipolygons with holes also work fine:
leaflet() %>%  addTiles() %>% addPolygons(data = lakey_ll)