在R中通过世界地图剪辑空间多边形答案

【问题标题】：Clip spatial polygon by world map in R在R中通过世界地图剪辑空间多边形
【发布时间】：2018-03-13 22:17:53
【问题描述】：

这是我第一次在 R 中进行任何类型的空间数据可视化，我被困在一个特定的问题上。我想根据世界地图裁剪一个空间多边形（由一系列纬度/经度坐标指定），以便删除与地图多边形重叠的多边形的任何部分。以下面代码中的内容为例，我想剪裁矩形空间多边形，以便只保留多边形的海洋部分。

我找到了如何保留两个空间多边形之间的交集的示例，但我想做相反的事情。也许有一种方法可以定义交集，然后从我希望剪辑的多边形中“减去”？

这可能是一个非常基本的问题，但我们将不胜感激任何提示！

指定纬度/经度数据：

x_coord <- c(25, 25, 275, 275)
y_coord <- c(20, -50, -50, 20)
xy.mat <- cbind(x_coord, y_coord)
xy.mat

转换为空间多边形对象：

library(sp)
poly = Polygon(xy.mat)
polys = Polygons(list(poly),1)
spatial.polys = SpatialPolygons(list(polys))
proj4string(spatial.polys) = CRS("+proj=longlat +datum=WGS84 +no_defs 
    +ellps=WGS84 +towgs84=0,0,0")

转换为空间多边形数据框并导出为 shapefile：

df = data.frame(f=99.9)
spatial.polys.df = SpatialPolygonsDataFrame(spatial.polys, df)
spatial.polys.df

library(GISTools)
library(rgdal)
writeOGR(obj=spatial.polys.df, dsn="tmp", layer="polygon", 
    driver="ESRI Shapefile")

绘制世界地图并添加 .shp 文件：

map("world", wrap=c(0,360), resolution=0, ylim=c(-60,60))
map.axes()
shp <- readOGR("polygon.shp")
plot(shp, add=TRUE, col="blue", border=FALSE)

【问题讨论】：

标签： r gis geospatial spatial

【解决方案1】：

这是一个始终保留在sf 中的解决方案（我不知道sp），并说明了从头开始构造sf 对象。 st_difference 精确创建您想要的几何图形，然后可以使用基本绘图方法或具有geom_sf 的ggplot 的开发版本进行绘图。为此，我使用了来自maps 和rnaturalearth 的地图数据，您可以根据自己的具体情况进行调整。不幸的是，环绕日期线有点挑剔。

library(tidyverse)
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.0, proj.4 4.9.3
library(rnaturalearth)
library(maps)
#> 
#> Attaching package: 'maps'
#> The following object is masked from 'package:purrr':
#> 
#>     map

x_coord <- c(25, 25, 275, 275)
y_coord <- c(20, -50, -50, 20)
polygon <-  cbind(x_coord, y_coord) %>%
  st_linestring() %>%
  st_cast("POLYGON") %>%
  st_sfc(crs = 4326, check_ring_dir = TRUE) %>%
  st_sf() %>%
  st_wrap_dateline(options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180"))

land <- rnaturalearth::ne_countries(returnclass = "sf") %>%
  st_union()
ocean <- st_difference(polygon, land)
#> although coordinates are longitude/latitude, st_difference assumes that they are planar

plot(st_geometry(land))
plot(st_geometry(polygon), add = TRUE)
plot(st_geometry(ocean), add = TRUE, col = "blue")

ggplot() +
  theme_bw() +
  borders("world") +
  geom_sf(data = ocean)

由reprex package (v0.2.0) 于 2018 年 3 月 13 日创建。

【讨论】：

【解决方案2】：

如果我正确理解您想要什么，您可以使用 sf 包使用 st_difference() 和 st_union()` 来实现。

根据您的代码，您可以执行以下操作。

# world data
data("wrld_simpl", package = 'maptools')

# load sf package
library('sf')

# coerce sp object to sf
world <- st_as_sf(wrld_simpl)
rectangle <- st_as_sf(spatial.polys)

# difference between world polygons and the rectangle
difference <- st_difference(rectangle, st_union(world))

# coerce back to sp
difference <- as(difference, 'Spatial')

# plot the result
plot(difference)

【讨论】：