2

Im trying to crop a pacific centred world map using the naturalearth maps data as an sf object. The map was created following this answer and the following code:


worldMap <- ne_countries(scale = "medium", returnclass = "sf") %>%
  st_make_valid()

target_crs <- st_crs("+proj=eqc +x_0=0 +y_0=0 +lat_0=0 +lon_0=133")

# define a long & slim polygon that overlaps the meridian line & set its CRS to match
# that of world
# Centered in lon 133

offset <- 180 - 133


polygon <- st_polygon(x = list(rbind(
  c(-0.0001 - offset, 90),
  c(0 - offset, 90),
  c(0 - offset, -90),
  c(-0.0001 - offset, -90),
  c(-0.0001 - offset, 90)
))) %>%
  st_sfc() %>%
  st_set_crs(4326)


# modify world dataset to remove overlapping portions with world's polygons
world2 <- worldMap %>% st_difference(polygon)
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries

# Transform
world3 <- world2 %>% st_transform(crs = target_crs)

ggplot(data = world3, aes(group = admin)) +
  geom_sf(fill = "grey")

pacific centered map

I would like to crop to an area of interest covering the Indian and Pacific Oceans (like this screenshot taken from the "maps" world2 projection):

desired crop extent

To crop, I tried the following:


world4 <- st_crop(world3, c(xmin= 30, ymin = -40, xmax = 230, ymax = 40))

However, the code above produces the incorrect crop with crazy small lat/long increments:

ggplot(data = world3, aes(group = admin)) +
  geom_sf(fill = "grey")

too small?!

If anyone can help me to figure out how to crop using lat and long to define the area for this pacific centred map it would be greatly apriciated. I know that there are other ways to obtain a pacific centred map but the future planned analysis for this project requires it to be an rnaturalearth map.

Gemma Galbraith
  • 55
  • 3

1 Answers1

2

The problem with your code is that you need to define the bbox for the crop filter using the same CRS as the world3 object. For example:

world4 <- st_crop(
  x = world3, 
  y = st_as_sfc(
    st_bbox(c(xmin= 30, ymin = -40, xmax = 230, ymax = 40), crs = 4326)
  ) %>% st_transform(target_crs)
)
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries
ggplot(data = world4, aes(group = admin)) +
  geom_sf(fill = "grey")

agila
  • 3,289
  • 2
  • 9
  • 20