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

1

6

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:

enter image description here

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!

Milly answered 16/11, 2023 at 20:11 Comment(5)
Have you tried this? rdrr.io/github/UW-Upwelling-Project/imageryML/src/R/… You would need to find alternatives for rgeos and spatialEco functions.Godden
Yes, I have tried that but I think the the spatialEco package is out of date? i also get this error as I move through it: Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented, : Loop 12 edge 1841 has duplicate near loop 16 edge 3"Milly
I have also tried many iterations of that function, but cannot get past the buffer stage. the error message is: Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented, : Loop 0 is not valid: Edge 0 crosses edge 78Milly
Yes, of course. You need to find alternatives. The given routine may help in approaching your goal.Godden
Ok, I understand but I am just not able to find the alternatives, I have been trying to find a workaround for days now but no luck. let me know if anything comes to mind! thank you.Milly
J
3

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")

enter image description here

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")

enter image description here

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.

Jurywoman answered 23/11, 2023 at 22:6 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.