将 ggplot 转换为地理参考 pdf 会产生镜像

Question

我正在尝试创建在 ggplot 中创建的地图的地理参考 pdf 版本，以便它们可以在平板电脑上的现场使用。 基于这篇文章sf：使用ggplot2输出对象创建GeoPDF似乎应该是可能的，但最终结果是原始ggplot地图的镜像。下面的代码复制了这个问题。我最初认为这个问题是由于我使用底图而引起的，这需要使用 CRS 3857 进行操作，但即使在 4326 (WGS84) 投影中复制，我仍然得到镜像。 我尝试更改提供的顺序范围，更改导出到 geotiff 函数的数据类型。 直到导出到 geotiff 后才会出现差异。它似乎非常接近工作，我想知道我的数据以及北美的 nc 数据是否具有负经度（以及 x 在 3857 中）值导致了问题。 任何见解将不胜感激。

library(sf)
library(raster)
library(gdalUtilities)
library(tidyverse)

nc <- st_read(system.file("shape/nc.shp", package="sf"),
              quiet=TRUE) %>%
  st_transform(4326)

nc_points <- nc %>%
  st_centroid()


GGPLOT_MAP <- 
  # basemaps::basemap_ggplot(
  #   ext = nc %>%
  #     st_buffer(5000),
  #   map_service = "esri" ,
  #   map_type = "world_topo_map",
  #   map_res = 1,
  #   dpi = 300,
  #   interpolate = T
  # )+
  ggplot()+
  geom_sf(
    data = nc,
    inherit.aes = F,
    aes(colour = NAME),
    fill = NA, lwd=1
  ) +
  geom_sf(
    data = nc_points,
    inherit.aes = F,
    aes(colour = as.character(NAME))
  ) +
  geom_sf_text(
    data = nc_points,
    inherit.aes = F,
    aes(label =NAME
    ),size = 1)+
  coord_sf(expand = F )+
  theme(legend.position = "none")

GGPLOT_MAP


ggsave(plot=GGPLOT_MAP, "gg.tiff", device = "tiff", dpi = 600)

# Create a StackedRaster object from the saved plot
stackedRaster <- stack("gg.tiff")

# Get the GeoSpatial Components
lat_long <- ggplot_build(GGPLOT_MAP)$layout$panel_params[[1]][c("x_range","y_range")] 

# Supply GeoSpatial  data to the StackedRaster 
extent(stackedRaster) <- c(lat_long$x_range,lat_long$y_range)
projection(stackedRaster) <- CRS("+proj=longlat +datum=WGS84")

# Create the GeoTiff
writeRaster(stackedRaster, "myGeoTiffgg.tif", options="PHOTOMETRIC=RGB", datatype="INT1U",overwrite=TRUE)

#using gdalutilities vs gdalUtil

#Save raster as GeoPDF
gdalUtilities::gdal_translate("myGeoTiffgg.tif","myGeoTiffgg.pdf",of="PDF", ot="Byte",
                          co="TILED=YES",  a_srs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

Answer 1

如果无法访问您的数据，我不能绝对肯定地说这就是问题所在，但它可能是因为初始 gg.tiff 中没有定义绘图尺寸。因此，代码的输出如下所示：

注意数据周围的空白以及轴的包含。换句话说，这些 lat_long 范围值已分配给完整绘图范围，而不是数据范围。

要纠正此问题，您可以使用

tmaptools::get_asp_ratio

定义输出图的宽度/高度。由于

raster

软件包已被弃用，以下使用其替代品 -

terra

。

library(terra)
library(sf)
library(gdalUtilities)
library(tmaptools)
library(ggplot2)

# Example data
nc <- st_read(system.file("shape/nc.shp", package="sf"), quiet = TRUE) %>%
  st_transform(4326)

nc_points <- nc %>%
  st_centroid() %>%
  select(NAME)

# Define extent for plot ratio
pr <- st_as_sfc(st_bbox(nc))  %>%
  st_as_sf(crs = 4326) %>%
  get_asp_ratio()

# Generate initial plot (modified for better output quality)
ggplot() +
  geom_sf(
    data = nc,
    aes(colour = NAME),
    fill = NA,
    lwd = 1
  ) +
  geom_sf(
    data = nc_points,
    aes(colour = NAME),
    size = 4
  ) +
  geom_sf_text(
    data = nc_points,
    aes(label = NAME),
    size = 5)+
  coord_sf(expand = F)+
  theme_void() +
  theme(legend.position = "none")

# Save plot using pr value
ggsave("gg1.tiff",
       bg = "white",
       width = 10 * pr, 
       height = 10,
       dpi = 600)

# Load gg1.tiff, assign extent and CRS values
stackedRaster <- rast("gg1.tiff")
ext(stackedRaster) <- ext(nc)
crs(stackedRaster) <- "EPSG:4326"

# Save modified stackedRaster
writeRaster(stackedRaster, "gg2.tiff", overwrite = TRUE)

# Save as georeferenced PDF
gdal_translate("gg2.tiff",
               "gg2.pdf",
               of = "PDF",
               ot = "Byte",
               co = "TILED=YES",
               a_srs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

无法发布 PDF 结果，但这是 gg2.tif 的一个版本（由于原始太大而无法共享，因此缩小了）：

我对照原始 nc sf 文件检查了 gg2.tif，它们按预期对齐。但是，我无法验证这是否是正确对齐的 GeoPDF。有问题请在下方评论。