Is it possible to colour outer space black when the projection view contains both Earth and outer space in ggplot with sf data?

I am making a few maps of the USA and some of them contains both Earth and outer space due to the projection, as shown below:

It is generated by the code below:

library(ggplot2)
library(sf)


# Import Data:
df <- read_sf("https://raw.githubusercontent.com/rfortherestofus/book/refs/heads/main/data/states.geojson")

# Create the map:
ggplot() +
  geom_sf(data = df, fill = "grey90") +
  
  coord_sf(crs = "+proj=ortho +lat_0=22.0 +lon_0=-162.5") +
  
  theme_minimal() +
  
  theme(
    panel.background = element_rect(fill = "lightblue")   # <- Set the colour for ocean
  )

I am just wondering if it is possible to colour outer space (highlighted) in black to reflect reality and also make the Earth, and therefore the USA, more distinguishable on the map?

One option would be to create a new sf object which is a grid covering the entire planet. You would need to transform that to your chosen projection and weed out the illegal panels. Thereafter you can simply plot it as a normal sf object:

library(ggplot2)
library(sf)

seas <- st_polygon(list(cbind(c(seq(-180, 180, len = 100),
                                rep(0, 100),
                                seq(180, -180, len = 100),
                                rep(-180, len = 100)), 
                              c(rep(-90, 100),
                                seq(-90, 90, len = 100),
                                rep(90, 100),
                                seq(90, -90, len = 100))))) |>
  st_sfc(crs = "WGS84") |>
  st_sf() |>
  st_make_grid(n = c(200, 200)) |>
  st_transform(crs = "+proj=ortho +lat_0=22.0 +lon_0=-162.5")

val_st <- st_is_valid(seas)
seas <- seas[!is.na(val_st) & val_st]

ggplot() +
  geom_sf(data = seas, fill = "lightblue", col = "lightblue",
          linewidth = 1) +
  geom_sf(data = df, fill = "grey90") +
  coord_sf(crs = "+proj=ortho +lat_0=22.0 +lon_0=-162.5",
           xlim = c(-513107, 6513107),
           ylim = c(-400000, 5000000)) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "black"),
        panel.grid = element_blank())

Note that this covers up the latitude / longitude lines, but you can explicitly add them back in (without the annoying artefacts) by specifying a multilinestring:

latlon <- st_multilinestring(x = list(cbind(rep(-120, 100), 
                                            seq(0, 90, len = 100)),
                                      cbind(rep(-60, 100), 
                                            seq(0, 90, len = 100)),
                                      cbind(rep(-180, 100), 
                                            seq(0, 90, len = 100)),
                                      cbind(seq(-180, 0, len = 100),
                                            rep(10, 100)),
                                      cbind(seq(-180, 0, len = 100),
                                            rep(20, 100)),
                                      cbind(seq(-180, 0, len = 100),
                                            rep(30, 100)),
                                      cbind(seq(-180, 0, len = 100),
                                            rep(40, 100)),
                                      cbind(seq(-180, 0, len = 100),
                                            rep(50, 100)),
                                      cbind(seq(-180, 0, len = 100),
                                            rep(60, 100)),
                                      cbind(seq(-180, 180, len = 100),
                                            rep(70, 100)))) |>
  st_sfc(crs = "WGS84") |>
  st_sf() |>
  st_transform(crs = "+proj=ortho +lat_0=22.0 +lon_0=-162.5")

which results in

ggplot() +
  geom_sf(data = seas, fill = "lightblue", col = "lightblue",
          linewidth = 1) +
  geom_sf(data = df, fill = "grey90") +
  geom_sf(data = latlon, linewidth = 0.2, col = "gray50") +
  coord_sf(crs = "+proj=ortho +lat_0=22.0 +lon_0=-162.5",
           xlim = c(-500000, 6000000),
           ylim = c(-400000, 5000000)) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "black"),
        panel.grid = element_blank(),
        panel.border = element_rect(fill = NA, colour = "black", linewidth = 1))

If you don't want to manually specify the co-ordinate limits you can get them from your sf object as follows:

crds <- df |> st_transform(crs = "+proj=ortho +lat_0=22.0 +lon_0=-162.5") |> 
  st_coordinates() |> 
  apply(2, range) 

ggplot() +
  geom_sf(data = seas, fill = "lightblue", col = "lightblue",
          linewidth = 1) +
  geom_sf(data = df, fill = "grey90") +
  geom_sf(data = latlon, linewidth = 0.2, col = "gray50") +
  coord_sf(crs = "+proj=ortho +lat_0=22.0 +lon_0=-162.5",
           xlim = crds[,1],
           ylim = crds[,2]) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "black"),
        panel.grid = element_blank(),
        panel.border = element_rect(fill = NA, colour = "black", linewidth = 1),
        axis.text.y = element_blank())

If you want more of a 3D appearance, the grid allows you to have an easy gradient fill, and I suppose you could add a gratuitous star field:

ggplot() +
  annotate("point", x = I(runif(300)), y = I(runif(300)), color = "white",
           alpha = runif(300), size = 0.3) +
  geom_sf(data = seas, aes(color = after_scale(fill),
            fill = sapply(seas[[1]], function(x) {
              sqrt(sum(colMeans(st_coordinates(x)[,1:2])^2))})),
          linewidth = 1) +
  geom_sf(data = latlon, linewidth = 0.2, col = "gray25") +
  geom_sf(data = df, fill = "#809060") +
  coord_sf(crs = "+proj=ortho +lat_0=22.0 +lon_0=-162.5",
           xlim = crds[,1],
           ylim = crds[,2]) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "black"),
        panel.grid = element_blank(),
        panel.border = element_rect(fill = NA, colour = "black", linewidth = 1),
        axis.text.y = element_blank(),
        axis.title = element_blank()) +
  scale_fill_gradient(low = "lightblue", high = "deepskyblue4", guide = "none")