Trimming (difference) polygon by ggplot2 states polygons in R-CodePudding

I am trying to trim (take difference of) a polygon by another polygon.

I've created SpatialPolygons with the package sp, and I can use rgeos::gDifference (from

We want to keep the red that is off the west coast (not overlapped by the states).

library("rgeos")
Real <- gDifference(p1,p2)

Here I get a huge output, with errors:

p2 is invalid
Input geom 1 is INVALID: Self-intersection at or near point -122.33795166 48.281641190312918 (-122.33795166000000165 48.2816411903129179)
<A>
...
Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false,  : 
  TopologyException: Input geom 1 is invalid: Self-intersection at -122.33795166 48.281641190312918
In addition: Warning messages:
1: In RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid") :
  Self-intersection at or near point -122.33795166 48.281641190312918
2: In gDifference(p1, p2) :
  Invalid objects found; consider using set_RGEOS_CheckValidity(2L)

I am assuming this is because the state polygon isn't one cohesive polygon (state borders are intact somehow and lines intersect one another, even though I only fed Lat/Lon to Polygon.

If I create a test polygon just by shifting the initial polygon it works:

 testPoly<-Polygon(as.matrix(cbind(WCG$Lon 1,WCG$Lat)))
 p3 <- SpatialPolygons(list(Polygons(list(testPoly), "p3")),proj4string=crdref)

plot(p1,col="red")
plot(p3,add=T)

Test <- gDifference(p1,p3)
plot(Test,col="red")

I've tried combining the "states" polygon with various forms of "union" functions (

library(rgeos)
trim<-gDifference(p1,coast_lines)

plot(trim)

CodePudding user response：

In case you are willing to switch from {sp} to {sf}, you might be looking for st_difference():

library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE
library(dplyr)

# load data
states <- ggplot2::map_data("state") 

# trim to continental west coast, create sf from df
states <- states[states$region=="washington"|states$region=="oregon"|states$region=="california",] |> 
  st_as_sf(coords = c("long", "lat"), crs = "epsg:4326") 

# point to polygon
states <- states |> 
  group_by(region) |> 
  summarise(geometry = st_combine(geometry)) |> 
  st_cast("POLYGON") |> 
  st_make_valid()

# load data for primary polygon
WCG <- structure(list(X = c(
  665004L, 665232L, 661983L, 663266L, 660980L,
  666562L, 660979L, 659316L, 661115L, 665803L, 663685L, 666535L,
  667728L, 660758L, 661000L, 665903L, 664469L, 659077L, 665725L
), Survey = c(
  "WCGBTS", "WCGBTS", "WCGBTS", "WCGBTS", "WCGBTS",
  "WCGBTS", "WCGBTS", "WCGBTS", "WCGBTS", "WCGBTS", "WCGBTS", "WCGBTS",
  "WCGBTS", "WCGBTS", "WCGBTS", "WCGBTS", "WCGBTS", "WCGBTS", "WCGBTS"
), Lat = c(
  33.07667, 32.91278, 32.70306, 32.57472, 32.0075, 31.99861,
  32.01028, 32.28583, 32.38222, 32.81528, 40.13528, 40.25611, 48.07222,
  48.175, 48.42278, 48.44444, 48.45083, 48.41556, 48.37028
), Lon = c(
  -117.3383,
  -117.2867, -117.2897, -117.3006, -118.3397, -118.6144, -118.8803,
  -119.6567, -119.885, -120.2967, -125.07, -125.1383, -125.9614,
  -125.9075, -125.1892, -124.9861, -124.8464, -124.7778, -124.755
), Date = c(
  "7/15/2011", "7/17/2012", "7/17/2012", "7/17/2015",
  "7/14/2010", "10/12/2016", "7/14/2010", "7/15/2007", "10/13/2010",
  "10/9/2017", "6/22/2015", "9/18/2016", "5/29/2019", "5/26/2010",
  "8/24/2010", "5/23/2017", "5/29/2009", "8/22/2007", "8/25/2017"
)), row.names = c(
  478258L, 478486L, 475237L, 476520L, 474234L,
  479816L, 474233L, 472570L, 474369L, 479057L, 476939L, 479789L,
  480982L, 474012L, 474254L, 479157L, 477723L, 472331L, 478979L
), class = "data.frame")

# df to sf again
WCG <- st_as_sf(WCG, coords = c("Lon", "Lat"), crs = "epsg:4326") 

# point to polygon
WCG <- WCG |> 
  group_by(Survey) |> 
  summarise(geometry = st_combine(geometry)) |> 
  st_cast("POLYGON") 

# inspect 
plot(st_geometry(states))
plot(st_geometry(WCG), border = "red", add = TRUE)

Note the topological artifact in Washington. I was hoping to remove this via st_make_valid(), but was not able to do so unfortunately. However, I hope this does not cause significant problems in demonstrating the workflow.

I assume this is what causes st_intersection() and st_difference() to fail with the following message, but there seems to be a workaround using sf_use_s2(FALSE) (c.f.

# difference
y <- st_combine(states) |> st_union() |> st_difference(WCG, y = _) 
plot(st_geometry(y))