在 ggplot2 中的多层 shapefile 地图上显示多边形内的多边形

Show polygons within polygons on a multi-layered shapefile map in ggplot2

我想使用位于彼此之上的多层 shapefile。我有三个深度级别的测深图(即 200m、1000m 和 2000m)。多边形中有一些 "holes"(*不是白地多边形),如下面我在 QGIS 中生成的地图(即 200 米内有 2000 米深度的多边形):

这是我使用相同的 shapefile 在 ggplot2 中制作的:

ggplot 地图上未显示多边形内的多边形。我怎样才能解决这个问题?

提前致谢。

用于地图的 R 脚本是:

Thailand <- readShapePoly("Thailand.shp")
Thailand2 <- fortify(Thailand)
Bthy_200m <- readShapePoly("ne_10m_bathymetry_K_200.shp")
Bthy_1000m <- readShapePoly("ne_10m_bathymetry_all/ne_10m_bathymetry_J_1000.shp")
Bthy_2000m <- readShapePoly("ne_10m_bathymetry_all/ne_10m_bathymetry_I_2000.shp")

Bthy_200m_crop <- crop(Bthy_200m, extent(84.11236, 108.4594, -4.046979, 24.09534))
Bthy_1000m_crop <- crop(Bthy_1000m, extent(84.11236, 108.4594, -4.046979, 24.09534))
Bthy_2000m_crop <- crop(Bthy_2000m, extent(84.11236, 108.4594, -4.046979, 24.09534))

Bthy_200m_crop2<- fortify(Bthy_200m_crop)
Bthy_1000m_crop2<- fortify(Bthy_1000m_crop)
Bthy_2000m_crop2<- fortify(Bthy_2000m_crop)

ggplot()+geom_polygon(data = Bthy_200m_crop2, aes(long, lat, group = group), fill="#52958b", na.rm =TRUE)+geom_polygon(data = Bthy_1000m_crop2, aes(long, lat, group = group), fill="#128277", na.rm =TRUE)+geom_polygon(data = Bthy_2000m_crop2, aes(long, lat, group = group), fill="#004D47", na.rm =TRUE)+geom_polygon(data = est_contour, aes(long, lat, group = group), fill="#99CCCC")+  geom_path(data = est_contour, aes(long, lat, group = group), color = "black")+theme_bw()+theme(panel.background = element_rect(fill = "#5EA8A7"))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

形状文件位于此处:http://www.naturalearthdata.com/downloads/10m-physical-vectors/10m-bathymetry/

所以我下载了一些提到的文件。要查看多边形中的多边形,我建议使用 geom_sf。这里有一些代码以及输出图片,以确保它是您所需要的。这是使用 geom_sf 生成此代码的代码。请注意,我没有使用泰国,因为我找不到 shapefile,也没有使用河口 shapefile,但您应该能够将此代码应用于绘图

library(sf)
library(tidyverse)
library(maptools)
library(rgdal)
library(raster)

Bthy_200m <- read_sf("ne_10m_bathymetry_K_200.shp")
Bthy_1000m <- read_sf("ne_10m_bathymetry_J_1000.shp")
Bthy_2000m <- read_sf("ne_10m_bathymetry_I_2000.shp")


Bthy_200m_crop <- st_intersection(Bthy_200m,st_set_crs(st_as_sf(as(extent(84.11236, 108.4594, -4.046979, 24.09534),"SpatialPolygons")), st_crs(Bthy_200m)))
Bthy_1000m_crop <- st_intersection(Bthy_1000m, st_set_crs(st_as_sf(as(extent(84.11236, 108.4594, -4.046979, 24.09534),"SpatialPolygons")), st_crs(Bthy_1000m)))
Bthy_2000m_crop <- st_intersection(Bthy_2000m, st_set_crs(st_as_sf(as(extent(84.11236, 108.4594, -4.046979, 24.09534),"SpatialPolygons")), st_crs(Bthy_2000m)))


Bthy_200m_crop2<- fortify(Bthy_200m_crop)
Bthy_1000m_crop2<- fortify(Bthy_1000m_crop)
Bthy_2000m_crop2<- fortify(Bthy_2000m_crop)
ggplot()+
  geom_sf(data = Bthy_200m_crop2,fill = "#52958b",na.rm = TRUE)+
  geom_sf(data = Bthy_1000m_crop2, fill = "#128277",na.rm = TRUE)+
  geom_sf(data = Bthy_2000m_crop2, fill="#004D47", na.rm =TRUE)+
  coord_sf(crs = st_crs(Bthy_1000m_crop2), datum = NA)+
  theme_bw()+
  theme(panel.background = element_rect(fill = "#5EA8A7"))+ 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())