ggplot2修剪(差异)多边形表示R中的多边形

l7wslrjt 于 2023-05-26 发布在其他

关注(0)|答案(2)|浏览(194)

我正试图用一个多边形修剪另一个多边形。
我已经用sp包创建了SpatialPolygons，我可以使用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")

我尝试过将“states”多边形与各种形式的“union”函数（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

这再次告诉我，我的状态多边形有一个自交。是否有更好的方法按北美大陆修剪我的原始多边形？
我在这里发布了另一个（相关）问题（How to use polyclip::polysimplify on a polygon），以了解如何摆脱这种自相交并创建一个实心多边形。
任何指针将不胜感激。

来源：https://stackoverflow.com/questions/74249703/trimming-difference-polygon-by-ggplot2-states-polygons-in-r

2条答案

按热度按时间

vm0i2vca1#

如果你想从{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)的解决方案（参见r-spatial/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, ])

赞(0）回复(0）举报 2023-05-26

zzzyeukh2#

我无法用ggplot2的原始状态数据集解决这个问题，这就是导致问题的原因。
相反，我按照这篇文章的开头（https://pjbartlein.github.io/REarthSysSci/RMaps2.html）从www.example.com获得了世界土地的空间多边形naturalearthdata.com（代码如下）
形状文件：https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_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)

赞(0）回复(0）举报 2023-05-26

我来回答

ggplot2修剪(差异)多边形表示R中的多边形

2条答案

相关问题

热门标签

最新问答