如何将栅格和多边形数据集添加到同一张地图，每个数据集都有单独的填充美学

Question

我正在尝试制作一张带有海洋测深和通过变量（每个区域不同颜色）着色（填充=填充）多边形的地图。我尝试了一些不同的软件包和方法，但我正在努力如何按区域进行着色/填充，同时在地图中进行测深。

我不确定如何共享我的数据以使其可重复，因为即使我只采样 1 个观察值，dput 返回也非常长。可以添加最有帮助的内容。简而言之，它是 NOAA 公开的数据。

library(sf)
library(tidyverse)
library(rnaturalearth) #for ne_states()
library(marmap) #for getNOAA.bathy()
library(ggOceanMaps) #for basemap()

那么，对于我的数据，我...

ecomap <- st_read("MYSHAPEFILE.shp") #how to share since dput even for 1 observation is very long
st_crs(ecomap) <- 4326 #set Coordinate Reference System CRS 
ecomap_sf <- st_as_sf(ecomap, wkt = "geom", crs = 4326)

custom_colors <- c("MAB" = "red", 
                   "SNE" = "blue", 
                   "GB" = "green", 
                   "GOM" = "orange")

o

> str(ecomap_sf)
Classes ‘sf’ and 'data.frame':  48 obs. of  7 variables:
 $ Name     : int  1 2 3 4 5 6 7 8 9 10 ...
 $ NumofPoly: int  1 1 1 1 1 1 1 1 1 1 ...
 $ NumofSta : int  1 2 2 1 5 2 2 4 1 3 ...
 $ Region   : chr  "MAB" "MAB" "MAB" "MAB" ...
 $ Area     : num  1540 4496 3190 2510 9620 ...
 $ Type     : chr  "SB" "MS" "IS" "SB" ...
 $ geometry :sfc_POLYGON of length 48; first list element: List of 1
  ..$ : num [1:387, 1:2] -74.8 -74.8 -74.7 -74.7 -74.7 ...
  ..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg"
 - attr(*, "sf_column")= chr "geometry"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA
  ..- attr(*, "names")= chr [1:6] "Name" "NumofPoly" "NumofSta" "Region" ...


> head(ecomap)
Simple feature collection with 6 features and 6 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -75.96791 ymin: 35.14219 xmax: -74.28861 ymax: 37.91444
Geodetic CRS:  WGS 84
  Name NumofPoly NumofSta Region     Area Type                       geometry
1    1         1        1    MAB 1540.432   SB POLYGON ((-74.76881 36.5052...
2    2         1        2    MAB 4496.485   MS POLYGON ((-74.81583 36.5052...
3    3         1        2    MAB 3190.414   IS POLYGON ((-75.72364 36.5060...
4    4         1        1    MAB 2510.323   SB POLYGON ((-74.3466 37.56465...
5    5         1        5    MAB 9620.144   MS POLYGON ((-74.63893 37.7118...
6    6         1        2    MAB 2726.496   IS POLYGON ((-75.54484 37.5418...

> unique(ecomap$Region)
[1] "MAB" "SNE" "GB"  "GOM"

我想向此地图添加测深功能（或类似地图）：

ggplot() +
  borders("world", color = "gray85", fill = "#575151") + 
  geom_sf(data = ecomap_sf1, aes(fill = Region), 
          color = "black", size = 0.5, alpha = 0.5) +    
  scale_fill_manual(values = custom_colors, name = "Region") + 
  coord_sf(xlim = c(-80, -62), ylim = c(34, 48))

我尝试过以下方法：

us <- ne_states(country = "united states of america", returnclass = "sf")

bathy <- getNOAA.bathy(lon1 = -81, lon2 = -62, 
                       lat1 = 24, lat2 = 45, 
                       resolution = 4)
bathy_df <- fortify.bathy(bathy)
bat_xyz <- as.xyz(bathy)


ggplot() +
  geom_sf(data = us) +
  geom_tile(data = bat_xyz, aes(x = V1, y = V2, fill = V3)) +
  geom_sf(data = us)+
  geom_sf(data = ecomap_sf, aes(fill = Region, color = "black"), #the issue is here
          color = "black", 
          size = 0.5, alpha = 0.5) +
  coord_sf(xlim = c(-80, -62), 
           ylim = c(34, 48))

但这会出现以下错误：

Error in `scale_fill_continuous()`:
! Discrete values supplied to continuous scale.
ℹ Example values: "MAB", "MAB", "MAB", "MAB", and "MAB"

如果我使用它，它会绘制：

aes(color = Region), 
          fill = "transparent"....

按区域（而不是整个区域/多边形）对边界线进行着色，这足以满足我的需要，但我仍然想知道为什么多边形中的填充/着色不能与测深同时工作

我尝试使用

fill=as.factor(Region)

和

fill=factor(Region)

代替，但仍然遇到相同的错误，即使在：

> class(ecomap$Region)
[1] "character"

或将其作为..

> class(ecomap$Region)
[1] "factor"

然后我使用 ggOceanMaps 包尝试了不同的方法

这绘制了地图（没有测深）

basemap(ecomap_sf) + 
  ggspatial::annotation_spatial(ecomap_sf, aes(fill = Region)) + 
  coord_sf(expand = FALSE)

但是如果我尝试添加

bathymetry=T

...

basemap(ecomap_sf, bathymetry = T) + 
  ggspatial::annotation_spatial(ecomap_sf, aes(fill = Region)) + 
  coord_sf(expand = FALSE)

给出以下错误：

Error in `palette()`:
! Insufficient values in manual scale. 12 needed but only 8 provided.

我尝试过：

basemap(limits = c(-82, -64, 30, 52),
        land.col = "gray21",
        land.border.col = NA,
        bathymetry = T,
        grid.col = NA) +
  geom_sf(data = ecomap_sf, aes(fill = Region), 
          color = "black", size = 0.5, alpha = 0.5) +    
  scale_fill_manual(values = custom_colors, name = "Region") + 
  coord_sf(xlim = c(-80, -62), ylim = c(34, 48))

绘制但以灰色“覆盖”测深。

我猜测测深颜色和多边形/几何颜色正在以一种不允许我控制其中一个的方式相互作用..？我不需要深度/等深线的图例。只是该地区的传奇。

我在网上研究了很多其他东西，但无法弄清楚。

来自在线人工智能聊天的解决方案有点有效，但我似乎无法控制太多大小，因为它相对于指定的限制将其绘制得非常小。

base_map1 <- basemap(data = ecomap_sf, 
                     bathymetry = TRUE, bathy.style = "rcb")


base_map1 +
  geom_sf(data = ecomap_sf[ecomap_sf$Region == "GOM", ], 
          fill = custom_colors["GOM"], 
          color = "black", size = 0.5, alpha = 0.2) +
  geom_sf(data = ecomap_sf[ecomap_sf$Region == "GB", ], 
          fill = custom_colors["GB"], 
          color = "black", size = 0.5, alpha = 0.2) +
  geom_sf(data = ecomap_sf[ecomap_sf$Region == "SNE", ], 
          fill = custom_colors["SNE"], 
          color = "black", size = 0.5, alpha = 0.2) +
  geom_sf(data = ecomap_sf[ecomap_sf$Region == "MAB", ], 
          fill = custom_colors["MAB"], 
          color = "black", size = 0.5, alpha = 0.2) +
  coord_sf(xlim = c(-80, -62), ylim = c(34, 48))

对我来说，这似乎是一种奇怪的方法，而且不是很简单

关于将测深添加到这样的地图并能够通过变量填充多边形的最佳方法是什么？有什么建议吗？

Answer 1

出现您的问题是因为您的 Region 值是离散字符串，但 V3 值是数字，是连续的。你不能把两者混为一谈。解决方案是使用

new_scale_fill()

包中的

ggnewscale

函数在每次调用时重置

aes(fill = )

。我组合并调整了各个图中的参数，调整了测深范围，并编辑了

coord_sf()

以避免 Bathy_xys 数据周围出现空白。修改此以满足您的需要：

library(sf)
library(dplyr)
library(rnaturalearth)
library(marmap)
library(ggOceanMaps)
library(ggnewscale)

# Create example sf
ecomap_sf <-  st_as_sfc(st_bbox(c(xmin = -75.96791,
                                  ymin = 35.14219,
                                  xmax = -65.28861,
                                  ymax = 44.91444))) |>
  st_make_grid(n = 4, crs = 4326) |>
  st_as_sf() |>
  mutate(Region = rep(c("MAB", "SNE", "GB",  "GOM"), each = 4))

custom_colors <- c("MAB" = "red", 
                   "SNE" = "blue", 
                   "GB" = "green", 
                   "GOM" = "orange")

us <- ne_states(country = "united states of america", returnclass = "sf")

bathy <- getNOAA.bathy(lon1 = -81, lon2 = -62, 
                       lat1 = 24, lat2 = 47, 
                       resolution = 4)
bathy_df <- fortify.bathy(bathy)
bat_xyz <- as.xyz(bathy)

ggplot() +
  geom_tile(data = bat_xyz, aes(x = V1, y = V2, fill = V3)) +
  geom_sf(data = us) +
  xlab("Longitude (decimal degrees)") + 
  ylab("Latitude (decimal degrees)") +
  guides(fill = guide_legend(title = "Depth(m)")) +
  new_scale_fill() +
  geom_sf(data = ecomap_sf, aes(fill = Region), 
          color = "black",
          size = 0.5,
          alpha = 0.5) +    
  scale_fill_manual(values = custom_colors) +
  coord_sf(xlim = c(-80, -62), 
           ylim = c(34, 46),
           expand = FALSE)