Creating a mask for natural earth R giving unexpected results

Question

I am wanting to create a mask to color all parts of the image that are NOT within the specified land areas. I can get this to work successfully for two countries (e.g Burkina Faso and Nigeria). However, the same code does not work when I try a longer list of 9 countries in the same area. The error I am receiving is : Error in [[<-.data.frame(tmp, attr(x, "sf_column"), value = list( : replacement has 2 rows, data has 1

Being relatively new to GIS, I am not quite sure what the issue is or how to resolve it. I would really appreciate some help if people could offer some insight. The method I used was based on the accepted answer from this question

Relevant code is below. Thankyou in advance for anyone who is able to help.

This code works perfectly fine

library(sf)
library(dplyr)
library(tidyverse)
library(rnaturalearth)
library(rnaturalearthdata)

# https://stackoverflow.com/questions/49266736/clip-spatial-polygon-by-world-map-in-r
x_coord = c(-20,-20,27,27)
y_coord = c(-1,28,28,-1)

# Create rectalgular polygon for a mask (slightly bigger than the image)
polygon <-  cbind(x_coord, y_coord) %>%
  st_linestring() %>%
  st_cast("POLYGON") %>%
  st_sfc(crs = 4326, check_ring_dir = TRUE) %>%
  st_sf() %>%
  st_wrap_dateline()

land = ne_countries(scale = "medium", 
             returnclass = "sf") %>% 
     filter(admin %in% c("Nigeria", "Burkina Faso"))
    #filter(admin %in% c("Niger", "Burkina Faso", "Nigeria", "Mali", "Chad", "Mauritania",
    #              "Gambia", "Senegal", "Guinea", "Cameroon"))


# Get the difference between the bounding box and the area I want to plot
# This works fine when using only two countries 
land = st_union(land)
polygon_land_diff = st_difference(polygon, land)

# Plot the data 
plot(st_geometry(land))
plot(st_geometry(polygon), add=TRUE)
plot(st_geometry(polygon_land_diff),add=TRUE, col="green")

This code produces an error and I am not sure how to resolve it. The only difference is the number of countries selected

x_coord = c(-20,-20,27,27)
y_coord = c(-1,28,28,-1)

# Create rectalgular polygon for a mask (slightly bigger than the image)
polygon <-  cbind(x_coord, y_coord) %>%
  st_linestring() %>%
  st_cast("POLYGON") %>%
  st_sfc(crs = 4326, check_ring_dir = TRUE) %>%
  st_sf() %>%
  st_wrap_dateline()

land = ne_countries(scale = "medium", 
             returnclass = "sf") %>% 
    filter(admin %in% c("Niger", "Burkina Faso", "Nigeria", "Mali", "Chad", "Mauritania",
                  "Gambia", "Senegal", "Guinea", "Cameroon"))


# get the difference between the bounding box and the area I want to plot
# not quite working. The rror occurs when running  `polygon_land_diff = st_difference(polygon, land)`

land = st_union(land)
polygon_land_diff = st_difference(polygon, land)

# Plot the data 
plot(st_geometry(land))
plot(st_geometry(polygon), add=TRUE)
plot(st_geometry(polygon_land_diff),add=TRUE, col="green")

Answer 1

st_difference() is generic function and I guess it faild to choice appropriate method at your latter case. (I didn't read source code, this is an estimate.)
Maybe it choice not sfc but sf method. (because Polygon is class sf)
you can get desired output by making polygon sfc class or using sf:::st_difference.sfc directly. Below is an example;

x_coord = c(-20,-20,27,27)
y_coord = c(-1,28,28,-1)


polygon <-  cbind(x_coord, y_coord) %>%
  st_linestring() %>%
  st_cast("POLYGON") %>%
  st_sfc(crs = 4326, check_ring_dir = TRUE)# %>%  # keep sfc class
#  st_sf() %>%
#  st_wrap_dateline()

land = ne_countries(scale = "medium", 
                    returnclass = "sf") %>% 
  filter(admin %in% c("Niger", "Burkina Faso", "Nigeria", "Mali", "Chad", "Mauritania",
                      "Gambia", "Senegal", "Guinea", "Cameroon"))

land = st_union(land)
polygon_land_diff = st_difference(polygon, land)

# Plot the data 
plot(st_geometry(land))
plot(st_geometry(polygon), add=TRUE)
plot(st_geometry(polygon_land_diff),add=TRUE, col="green")