创建由三个独立的形状文件组成的单个地图

Question

我有以下形状文件：

qmasp <- st_read("Rock_Lobster_(Spiny_Red)_QMAs.shp", quiet = TRUE)
statShp <- st_read("Rock_Lobster_Statistical_Areas.shp", quiet = TRUE)
NE_countries <- st_read("ne_110m_admin_0_countries", quiet = TRUE)

我希望它们作为一张地图一起出现。每个形状文件看起来像这样（参见图片）：

我正在使用以下代码：

将 QMA shapefile 转换为 sf 对象

qmasf=st_as_sf(qmasp)
qmasf4326=st_transform(qmasf, 4326)

然后 NZ UTM 进行绘图

qmasf2134=st_transform(qmasf, 2193)#2134
qmasf2134 %>% 
  ggplot()+
  geom_sf(aes(fill=FishstockC))

将 Stat Area Shapefile 转换为 sf 对象

statsf=st_as_sf(statShp)
statsf4326=st_transform(statsf, 4326)

然后 NZ UTM 进行绘图

statsf2134=st_transform(statsf, 2134)#check if NZ2000 is better
stats<-statsf2134 %>%
 ggplot()+
  geom_sf(fill= 'White')+
  geom_sf_text(
    aes(label = Statistica),
    size = 4 / .pt
  )
stats

情节地图

theme_set(theme_bw())

newzealand <- ne_countries(country="new zealand", type="countries", scale = 'large', returnclass = "sf")

ggplot(data = newzealand) +
  geom_sf() +
  coord_sf(crs = 2193)

然后，我可以使用以下代码，使用坐标数据将数据点绘制到地图上 QMA 形状文件中的特定多边形上：

准备在地图上绘制的数据

samples_sf <- sf::st_as_sf(d2, coords=c("Longitude","Latitude"), crs=4326)

QMA 形状文件中的 CRA3 多边形子集

CRA3_shapefile<-qmasf2134 %>%
  filter(FishstockC %in% c("CRA2","CRA3"))
CRA3_shapefile

看起来像这样（见下文）：

或者这取决于我使用的代码：

我的问题是，(i) 如何生成一张地图，该地图在一张图像上包含所有三个形状文件，并且基本上彼此层叠？ (ii) 为什么我的新西兰地图这么小，如何使其与其他 shapefile 大小相同，而又不丢失查塔姆岛（新西兰第三岛）？这个岛位于新西兰北岛和南岛的东部。

任何帮助将不胜感激！

Answer 1

Tēnā koe e hoa，鉴于我们无法访问您的数据，这是一个代表。您可以将

coord_sf()

视为缩放功能。要获取数据的范围，请使用

st_bbox()

并从那里继续前进。通过一些尝试和错误，您将获得所需的结果。示例图像使用如上所述的 xlim/ylim（长/纬度）。对于每一层，添加一个新的

geom_sf()

，如下所示。您每次都必须明确说明

data =

，但我更喜欢这种方法，因为我发现它更容易遵循。

geom_sf()

按顺序绘制，因此背景图层先行：

library(rnaturalearth)
library(sf)
library(dplyr)
library(ggplot2)

newzealand <- ne_countries(country = "new zealand", 
                           type = "countries", 
                           scale = 'large', 
                           returnclass = "sf") %>%
  st_transform(2193)

# Dummy CRA3_shapefile
CRA3_shapefile <- data.frame(cra = rep(c("CRA2", "CRA3"), each = 4),
                             lat = c(5700000, 5800000, 5800000, 5700000,
                                     5800000, 5900000, 5900000, 5800000),
                             long = c(2000000, 2000000, 2200000, 2200000)) %>%
  st_as_sf(coords = c("long","lat")) %>%
  group_by(cra) %>%
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("POLYGON") %>%
  st_set_crs(2193)

# Dummy point data. set.seed() for replicable result, only needed for dummy point data
set.seed(123)
sf_points <- data.frame(lat = sample(5700000:5900000, 20),
                        long = sample(2000000:2200000, 20)) %>%
  st_as_sf(coords = c("long","lat")) %>%
  st_cast("POINT") %>%
  st_set_crs(2193)

# Get extent of newzealand data as basis for setting coord_sf limits
st_bbox(newzealand)
xmin    ymin    xmax    ymax 
1092701 4165390 3357853 9024853 

ggplot() +
  geom_sf(data = CRA3_shapefile,
          aes(fill = cra)) +
  geom_sf(data = newzealand,
          colour = "black") +
  geom_sf(data = sf_points,
          colour = "grey40",
          size = 1.5) +
  coord_sf(xlim = c(1002701, 2557853),
           ylim = c(4165390, 6204853))

根据OP的评论进行更新

通常更需要过滤掉不必要的几何图形，而不是使用

coord_sf()

。在这种情况下，情况会更加棘手，因为您的 NZ 数据是多多边形，因此每个单独的多边形都将共享相同的非几何值（ArcGIS 术语中的属性）。换句话说，除了长/纬度值之外，无法知道哪些多边形代表 Avaiki Nui/库克群岛。

在您的情况下，最简单的方法是“借用”上一步中的最大 ylim 值，例如 6204853 并过滤掉 6204853 以北的所有内容。此方法涉及使用

st_cast()

:

将 MULTIPOLYGON 分解为单个多边形

# Convert MULTIPOLYGON TO POLYGON (you can safely ignore the resultant warning)
nz1 <- newzealand %>%
  st_cast("POLYGON") %>%
  mutate(id = 1:n()) # This column added as a key for filtering

# 1. Extract coordinates from nz1
# 2. Filter out polygon ids north of max ylim value
# 3. left_join() polygon ids to nz1
# 4. Filter out unwanted polygons
# NOTE: the Y and L2 columns referenced here are 
#       default names generated by data.frame(st_coordinates(nz1))
nz2 <- data.frame(st_coordinates(nz1)) %>%
  filter(Y <= 6204853) %>%
  select(L2) %>%
  distinct() %>%
  left_join(nz1, ., by = join_by(id == L2), keep = TRUE) %>%
  filter(!is.na(L2))