Sf: `st_intersect`: Errors in st_intersect

Created on 18 Oct 2018  路  2Comments  路  Source: r-spatial/sf

Hi there!

First, thanks for developing sf.
I wanted to mention an issue I am experiencing while using st_intersect().
I am trying to overlap a grid of hundreds of circles with countries' polygons, available from the map_data() function from ggplot2 v3. My intention is to see how many countries overlap each circle. The way I came around this was using st_intersect(), as its output for each circle consists of a data.framewhere the number of rows tells the number of polygons overlapped i.e number of countries.
However, in some cases, this method does not work, and I get the following errors

Error in CPL_geos_op2(op, st_geometry(x), st_geometry(y)) : Evaluation error: TopologyException: Input geom 0 is invalid: Nested shells at or near point...
or
Error in CPL_geos_op2(op, st_geometry(x), st_geometry(y)) : Evaluation error: TopologyException: Input geom 0 is invalid: Self-intersection at or near point ...
In the error messages above, I replaced the actual point with ....

Let me attach an R notebook (.Rmd file) showing a detailed MWE for this problem. I'll put below a summarised version of the notebook. All I have been able to figure out is that the nested shells error disappears if the circle size is changed. But then different circles can get the same error, so there's no unique size that works with all circles. I have no clue about the second error. Thanks for looking into this.

library(tidyverse)
library(ggplot2)
library(dplyr)
library(sp)
library(raster)
library(rgdal)
library(rgeos)
library(maptools)
library(maps)
library(rvest)
library(sf)
library(units)
library(mapview)
library(ggmap)

# Function to define circles based on a long-lat dataframe
# giving the circle ceners and a radius in kilometers
create_circle <- function(center_df, radius=100) {
  dat_sf <- st_as_sf(center_df, coords = c("lon", "lat"), crs = 4326) %>%
    # transform to an equal area projection so that the circle stays the same regardless
    # of the coordinates
    st_transform(4088)
  r <- set_units(radius, km)
  cent <- st_transform(dat_sf, 4088)
  cent_buffer <-  st_buffer(cent, r)
  cent_buffer <- st_transform(cent_buffer, 4326)
  return(cent_buffer)
}
world_map_sf <- sf::st_as_sf(maps::map("world", plot = FALSE, fill = TRUE))

center_df <- data.frame(lon = c(45), lat=c(10))
circles <- create_circle(center_df,radius=1)
st_intersection(world_map_sf, circles)

This results in the circle below, which the error indicates cannot be intersected with the country polygon (happens to be Region Somalia, subregion Somaliland, but it is irrelevant).
error

Most helpful comment

Your world_map_sf contains invalid geometries. Try

library(lwgeom)
st_intersection(st_make_valid(world_map_sf), circles)

All 2 comments

Your world_map_sf contains invalid geometries. Try

library(lwgeom)
st_intersection(st_make_valid(world_map_sf), circles)

Yes, it fixes the issue!
Using this trick I get a table where I get a row for every country and circle overlapped, i.e there's no grouping variable that tells me which countries are overlapped by the same circle.

However, that kind of output can be retrieved with st_intersects(), which besides does not need to be passed valid geometries, at least in my case. This function returns a list where for the ith country, I get in the ith element a numeric vector with the index of the circles that overlap it, so one can easily do the match.

Thank you
Antonio

Was this page helpful?
0 / 5 - 0 ratings