用R中的ggplot2态多边形裁剪(差分)多边形

Question

我正试图用另一个多边形修剪(取不同的)一个多边形。

我用包SpatialPolygons创建了sp，我可以使用rgeos::gDifference (来自https://gis.stackexchange.com/a/54146/171788)来复制(但移动)多边形，但是对于来自ggplot2的状态的多边形不起作用(参见下面)。

## Load data
library("ggplot2") 
states <- map_data("state") ## load state data
states<-states[states$region=="washington"|states$region=="oregon"|states$region=="california",]  ## trim to continental west coast

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

可能有一种更简单的方法可以互相减除多边形，但是sp::Polygon不适用于rgeos::gDifference，所以我转换为SpatialPolygons

library(sp)
WCGPoly<-Polygon(as.matrix(cbind(WCG$Lon,WCG$Lat)))
statesPoly<-Polygon(as.matrix(cbind(states$long,states$lat)))

crdref <- CRS('+proj=longlat +datum=WGS84')

p1 <- SpatialPolygons(list(Polygons(list(WCGPoly), "p1")),proj4string=crdref)
p2 <- SpatialPolygons(list(Polygons(list(statesPoly), "p2")),proj4string=crdref)

为了形象化我们试图保持的差异：

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

​

​

我们希望保留位于西海岸的红色(而不是与各州重叠)。

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

在这里，我得到了一个巨大的输出，带有错误：

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)

我假设这是因为状态多边形不是一个内聚多边形(状态边界在某种程度上是完整的，直线彼此相交，即使我只将Lat/Lon提供给Polygon。

如果我只是通过移动初始多边形来创建一个测试多边形，它可以工作：

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

​

​

我试着将“状态”多边形与各种形式的“联合”函数(https://gis.stackexchange.com/a/63696/171788，https://gis.stackexchange.com/a/385360/171788)组合在一起，但这对我来说没有任何意义，因为大多数示例使用的是两个不同的多边形，而不是像这些状态那样具有多个“分离”区域的多边形。

尝试这个不同的答案(https://gis.stackexchange.com/a/169597/171788)：

x<- p1-p2
x<- erase(p1,p2)

我得到以下错误：

x is invalid
Attempting to make x valid by zero-width buffering
Warning message:
In RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid") :
  Self-intersection at or near point -122.33795166 48.281641190312918

它再次告诉我，我的状态多边形有一个自交。有没有更好的方法来修剪我原来在北美大陆的多边形呢？

我在这里发布了另一个(相关的)问题(在多边形上)，以了解如何摆脱这个自交，创建一个实体多边形。

任何指示都将不胜感激。

Answer 1

如果您愿意从{sp}切换到{sf}，您可能正在寻找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)

​

​

注意华盛顿的拓扑伪像。我本来希望通过st_make_valid()删除它，但不幸的是我无法做到这一点。不过，我希望这不会对工作流程的演示造成重大问题。

我认为这是导致st_intersection()和st_difference()失败的原因，但是使用sf_use_s2(FALSE) (c.f )似乎有一个解决办法。R-空间/sf#1762)。

Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented,  : 
  Loop 0 is not valid: Edge 344 crosses edge 346 

因此，当st_intersection()返回两个多边形数据集之间的重叠区域时，st_difference()返回不重叠的区域，据我所知，这正是您所要寻找的。

sf_use_s2(FALSE)

# intersect
x <- st_intersection(states, WCG)
plot(st_geometry(x))

​

​

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

​

​

更新：

我不得不进一步挖掘，我想说，map_data()所获得的几何图形似乎从一开始就被打破了。

# cast MULTIPOLYGON to POLYGON
wash <- states[3, ] |> st_cast("POLYGON")

wash
#> Simple feature collection with 5 features and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -124.6813 ymin: 45.53296 xmax: -116.9235 ymax: 49.00508
#> Geodetic CRS:  WGS 84
#> # A tibble: 5 × 2
#>   region                                                                geometry
#>   <chr>                                                            <POLYGON [°]>
#> 1 washington ((-123.0198 48.56963, -123.0542 48.61547, -123.0943 48.60974, -123…
#> 2 washington ((-123.0198 48.56963, -122.825 48.40347, -122.7677 48.60974, -122.…
#> 3 washington ((-122.825 48.40347, -122.7963 48.40921, -122.7677 48.40921, -122.…
#> 4 washington ((-122.7677 48.60974, -122.7505 48.63839, -122.7161 48.65557, -122…
#> 5 washington ((-116.9521 46.06581, -116.9235 46.16894, -116.9579 46.20905, -116…

# get coordinates of doubtful polygon
c <- wash[2, ] |> st_coordinates() 

head(c)
#>              X        Y L1 L2
#> [1,] -123.0198 48.56963  1  1
#> [2,] -122.8250 48.40347  1  1
#> [3,] -122.7677 48.60974  1  1
#> [4,] -122.6703 48.17429  1  1
#> [5,] -116.9292 45.99705  1  1
#> [6,] -123.0198 48.56963  1  1

plot(c)
lines(c, type = "l", col = "red")

您可以清楚地看到，至少有四个移位点，因此产生了工件，特别是因为第五个点(-116.9292，45.99705)。

​

​

然而，我不明白为什么这些(你仍然可以看到小三角形.)当你绘制sf对象而不是它的几何图形时消失。我以为你只会以这种方式丢弃属性，而不会改变(只是)几何信息。

plot(states[3, ])

​

​

Answer 2

我无法用来自ggplot2的原始状态数据集来解决这个问题，这就是导致问题的原因。

相反，我遵循了这篇文章(https://pjbartlein.github.io/REarthSysSci/RMaps2.html)的开头，从naturalearthdata.com (代码如下)中获得了一个世界陆地的空间多边形。

形状文件：land.zip

## download shape file: https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_land.zip

library(rgdal)
layer <- ogrListLayers("../../../Downloads/ne_10m_land/ne_10m_land.shp")
ogrInfo("../../../Downloads/ne_10m_land/ne_10m_land.shp", layer=layer)
coast_lines <- readOGR("../../../Downloads/ne_10m_land/ne_10m_land.shp", layer=layer)

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


library(sp)
WCGPoly<-Polygon(as.matrix(cbind(WCG$Lon,WCG$Lat)))

crdref <- CRS('+proj=longlat +datum=WGS84')
p1 <- SpatialPolygons(list(Polygons(list(WCGPoly), "p1")),proj4string=crdref)

plot(p1)
plot(coast_lines,add=T)

​

​

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

plot(trim)

​

​