Create buffer from coastline 25km in R and calculate nearshore area

Question

I am trying to create a map containing polygons that represent 25km buffer into the ocean from every country's coastline, so that I can calculate the area within this buffer for each country. I am having trouble doing this both with an uploaded coastline file, and with naturalearth data. Here's what I've tried so far:

library(rgdal)
library(sf)
library(rnaturalearthdata)
library(rnaturalearth)
library(rgeos)
library(maps)
library(countrycode)

# World polygons from the maps package

world_shp <- sf::st_as_sf(maps::map("world", plot = T, fill = TRUE))

world_shp_df <- world_shp %>%
  mutate(Alpha.3.code = countrycode(ID, "country.name", "iso3c", warn = FALSE))

# Get coastlines using rnaturalearth
coastline <- ne_coastline(scale = "medium", returnclass = "sf")

# Buffer the coastlines by 25 km
buffer_distance <- 25000  # 25 km in meters
coastline_buffers <- st_buffer(coastline, dist = buffer_distance) %>%
  st_make_valid() 

ggplot() +
  geom_sf(data = coastline_buffers , fill = "lightblue")

This results in a map with horizontal lines running across it, see the image:

I have tried simplifying the geometry, using a different crs, and I just can't seem to figure it out. Any help would be really appreciated!

Answer 1

When the buffer crosses 180° latitude, st_buffer limits latitude to -180°, causing an inversion of what is inside and out of the buffer polygon, and thus horizontal stripes.

The polygons with stripes can be sorted out by their bounding box, because their xmin = -180 :

stripes <- sapply(coastline_buffers$geometry,\(x) {unname(st_bbox(x)$xmin == -180)})
which(stripes)
#[1]    1   61  138  297  298  299  811 1352 1353 1387 1388 1404

ggplot() +
  geom_sf(data = coastline_buffers[stripes,] , fill = "lightblue")

As there aren't too many faulty points, a quick workaround is to remove these faulty points from the polygons :

# Create a band from latitudes -175° to 180° 
cleanup <- data.frame(lon=c(-175,-175,180,180,-175),lat=c(90,-90,-90,90,90)) %>% 
           st_as_sf(coords = c("lon", "lat"), 
           crs = st_crs(p[1,])) %>% 
           st_bbox() %>% 
           st_as_sfc()

# Remove points in this band
coastline_buffers_cleaned <- st_difference(coastline_buffers,cleanup)


ggplot() +
  geom_sf(data = coastline_buffers_cleaned , fill = "lightblue")

This cleanup is not 100% exact, as some coastal area at both ends of the map disappeared, but stripes are away and removed area is negligible in comparison to total buffer area.