United States missing using geostationary satellite view projection only

[TeresaPegan]

Hello, I am using rnaturalearth for its polygons of the western hemisphere. I have typically used lambert equal area projections in the past, but I am trying to make an image using the Geostationary satellite view projection (+proj="geos"). When I do this, the United States is missing.

I am attaching reproducible code and images illustrating the issue. The "geos" projection shows the missing US polygon, while it is present in the version of the image made with the "laea" projection.

Thanks! -Teresa

library("tidyverse") # version 2.0.0
library("rnaturalearth") # version 1.0.1
library("rnaturalearthdata") # version 0.1.0
library("sf") # version 1.0-14

world <- ne_countries(continent = c("north america","south america"), scale = "medium", returnclass = "sf")
  worldGeos <- st_transform(world, crs="+proj=geos +h=35785831.0 +lat_0=45 +lon_0=-90.476194 +sweep=y")
  worldLaea <- st_transform(world, crs="+proj=laea +lat_0=-30 +lon_0=-90.476194 +x_0=4321000 +y_0=3210000 +ellps=WGS84 +units=m +no_defs")

  ggplot(data = worldGeos) +
    geom_sf(color = "lightgray", fill = "lightgray") +
    ggtitle("Geostationary Satellite View")

  ggplot(data = worldLaea) +
    geom_sf(color = "lightgray", fill = "lightgray") +
    ggtitle("Lambert Azimuthal Equal Area")

[PMassicotte]

I am also getting the same results. I am asking around for some pointer and I will get back to you.

[PMassicotte]

Does this work for you?

Basically, I am creating points regularly spaced over your area. Then, I project these points into your coord system and remove points that are out of view (NAs). Finally, I create a convex hull and use it to crop world before the projection.

library(tidyverse)
library(rnaturalearth)
library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.3, PROJ 9.1.1; sf_use_s2() is TRUE
library(s2)

proj <- "+proj=geos +h=35785831.0 +lat_0=45 +lon_0=-90.476194 +sweep=y"

world <- ne_countries(
  continent = c("north america", "south america"),
  scale = "medium",
  returnclass = "sf"
)

pts <- world |>
  # st_bbox() |>
  st_sample(size = 2e3, method = "regular") |>
  st_transform(proj)

ch <- pts[!st_is_empty(pts), drop = FALSE] |>
  s2_convex_hull() |>
  st_as_sf() |>
  st_transform(st_crs(world))

world2 <- world |>
  st_crop(ch)
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries

worldGeos <- st_transform(world2, crs = proj)

ggplot(data = worldGeos) +
  geom_sf(color = "lightgray", fill = "lightgray") +
  ggtitle("Geostationary Satellite View")

^{Created on 2024-02-12 with reprex v2.1.0}

[TeresaPegan]

That does indeed work for me, thank you! I think it makes sense to leave the issue open for now, since the underlying bug remains?

[PMassicotte]

I am not sure if it is a bug of sf or the underlying implementation. Tagging @edzer just in case this is related to sf.

[PMassicotte]

Much simpler solution (credit to @mdsumner)

library(ggplot2)
library(rnaturalearth)
library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.0, PROJ 9.2.0; sf_use_s2() is TRUE

world <- ne_countries(
  continent = c("north america", "south america"),
  scale = "medium",
  returnclass = "sf"
)

world <- st_cast(world, "POLYGON")
#> Warning in st_cast.sf(world, "POLYGON"): repeating attributes for all
#> sub-geometries for which they may not be constant

worldGeos <- st_transform(
  world,
  crs = "+proj=geos +h=35785831.0 +lat_0=45 +lon_0=-90.476194 +sweep=y"
)
ggplot(data = worldGeos) +
  geom_sf(color = "lightgray", fill = "lightgray") +
  ggtitle("Geostationary Satellite View")

^{Created on 2024-02-12 with reprex v2.1.0}