安装 rgdal 和使用关联的 readOGR() 函数时出现问题

Question

使用 R 版本 4.3.2，我使用

getData()

包中的

raster

函数来读取数据集，以在

ggplot2

中制作“KUD 图”。

但是，我有一个带有几何形状的多边形形状文件（.shp），我想用

readOGR()

读取它，但我无法安装

rgdal

:

install.packages("rgdal")

install.packages 中的警告： 软件包“rgdal”不适用于此版本的 R 适用于您的 R 版本的此软件包的版本可能在其他地方可用， 请参阅以下位置的想法 https://cran.r-project.org/doc/manuals/r-patched/R-admin.html#Installing-packages

我尝试使用旧版本的 R 安装

rgdal

，但它不起作用。

我的问题是需要使用

geom_polygon()

添加海岸线。我有一个现有的 SpatialPolygonsDataFrame，但它不够详细。这是我尝试过的：

ggplot() +
  coord_sf(crs = 32632)+
  geom_sf(data = st_as_sf(ver95), color= NA, fill="#3B9AB2", alpha=0.5)+
  geom_sf(data = st_as_sf(ver90), color= NA, fill="#48A0B6", alpha=0.5)+
  geom_sf(data = st_as_sf(ver85), color= NA, fill="#56A6BA", alpha=0.5)+
  geom_sf(data = st_as_sf(ver80), color= NA, fill="#63ADBE", alpha=0.5)+
  geom_sf(data = st_as_sf(ver75), color= NA, fill="#71B3C2", alpha=0.5)+
  geom_sf(data = st_as_sf(ver70), color= NA, fill="#84B9B3", alpha=0.5)+
  geom_sf(data = st_as_sf(ver65), color= NA, fill="#9EBE91", alpha=0.5)+
  geom_sf(data = st_as_sf(ver60), color= NA, fill="#B7C26E", alpha=0.5)+
  geom_sf(data = st_as_sf(ver55), color= NA, fill="#D1C74C", alpha=0.5)+
  geom_sf(data = st_as_sf(ver50), color= NA, fill="#EBCC2A", alpha=0.5)+
  geom_sf(data = st_as_sf(ver45), color= NA, fill="#E8C520", alpha=0.5)+
  geom_sf(data = st_as_sf(ver40), color= NA, fill="#E6BF17", alpha=0.5)+
  geom_sf(data = st_as_sf(ver35), color= NA, fill="#E4B80E", alpha=0.5)+
  geom_sf(data = st_as_sf(ver30), color= NA, fill="#E2B204", alpha=0.5)+
  geom_sf(data = st_as_sf(ver25), color= NA, fill="#E29E00", alpha=0.5)+
  geom_sf(data = st_as_sf(ver20), color= NA, fill="#E67D00", alpha=0.5)+
  geom_sf(data = st_as_sf(ver15), color= NA, fill="#EA5C00", alpha=0.5)+
  geom_sf(data = st_as_sf(ver10), color= NA, fill="#EE3B00", alpha=0.5)+
  geom_sf(data = st_as_sf(ver05), color= NA, fill="#F21A00", alpha=0.5)+
  geom_sf(data = st_as_sf(mcp), color= "black", fill=NA, alpha=0.5)+
 geom_point(data = Stations, aes(x=station_longitude, y = station_latitude),
            shape=21, color="red", size = 2) +
  geom_point(data = COA, aes(x=Longitude.coa, y = Latitude.coa), shape=3, color="black", size = 1) + ### points COAs
  ylab("Latitude") +
  xlab("Longitude")+
  coord_sf(xlim = c(9.43, 9.55195), ylim = c(42.551657, 42.68), expand = FALSE)+   
  theme(legend.position = "none") +
  geom_polygon(data = Corsica , aes(x=long, y=lat))+

  ggtitle("KUD blue crab")+ theme_bw()

这就是我得到的：

这是没有海岸线的情节，以显示为什么我需要

rgdal

来阅读更详细的海岸线，例如一些数据变得模糊：

Answer 1

您遇到的主要问题是您需要为多边形设置

fill = NA

。这是因为默认的

ggplot2

行为是填充多边形。此外，正如您所观察到的，您的数据不够详细，因此一些 KUD 数据被模糊化。另外，我强烈建议在绘制空间数据时使用

geom_sf()

，因为使用其他几何选项 可能会导致问题。

正如评论中提到的，

rgdal

已于 2023 年 10 月弃用。此代表使用

sf

包，大多数相关

rgdal

函数已迁移到该包。使用

sf

更加面向未来，并且可以避免您的次要问题，例如安装

rgdal

。请注意，如果您之前安装过

rgdal

，它可能会 导致问题，因此您可能需要从库目录中删除它的所有痕迹。这包括

GDAL

、

GEOS

和

PROJ

绑定到的

rgdal

的版本。

我设法从此链接获取了您用于绘图的数据集，并提供了一个工作流程来实现您想要的效果。它涉及创建数据的 LINESTRING 和 POLYGON 版本，然后将两者裁剪到您感兴趣的区域。您需要首先将下载内容解压缩到您的工作目录中，并在必要时更新您的包。

请注意，如果您想在

coord_sf()

内使用 xlim/ylim 值，则需要将当前的 WGS84:EPSG4326 xlim/ylim 值转换为等效的 EPSG:32632 值。我通过在绘图之前将海岸线数据裁剪为原始 xlim/ylim 值来解决这个问题。但是，如果您的其他数据超出了这些 xlim/ylim 值，您将需要相应地编辑代码。

library(dplyr) # 1.1.4
library(sf) # 1.0.15
library(ggplot2) # 3.5.0

# Data downloaded from https://sdi.eea.europa.eu/datashare/s/esBJK983EZpAZJg/download
# and unzipped to getwd()
Corsica_line <- st_read("EEA_Coastline_20170228.shp") %>%
  st_cast("LINESTRING") %>%
  st_transform(4326) %>%
  st_crop(xmin = 9.43, ymin = 42.551657, xmax = 9.55195, ymax = 42.68) %>%
  st_transform(32632)

Corsica_poly <- st_read("EEA_Coastline_20170228.shp") %>%
  st_transform(4326) %>%
  st_crop(xmin = 9.43, ymin = 42.551657, xmax = 9.55195, ymax = 42.68) %>%
  st_transform(32632)

ggplot() +
  geom_sf(data = Corsica_poly,
          fill = "grey70",
          colour = NA) +
  geom_sf(data = Corsica_line,
          colour = "grey30") +
  coord_sf(expand = FALSE) +
  # coord_sf(xlim = c(9.43, 9.55195),
  #          ylim = c(42.551657, 42.68),
  #          expand = FALSE) +
  theme(axis.text = element_text(size = 5),
        axis.text.x = element_text(angle = 90, hjust = 1))

结果：