带有 Lambert 投影的 geom_contour答案

【问题标题】：geom_contour with Lambert projection带有 Lambert 投影的 geom_contour
【发布时间】：2022-08-03 21:33:11
【问题描述】：

我正在尝试绘制东北大西洋的地图，突出显示特定的等深线。鉴于该区域，我更喜欢使用 LAEA 投影。但是，这似乎会导致 geom_contour 失败。

下面是后台代码：

library(tidyverse)
library(marmap)
library(rnaturalearth)

#Set boundaries
bbox <- tibble(x = c(-20, 10), y = c(45, 60))
#Add coordinated converted to LAEA 
bbox <- bbox %>%
      bind_cols(bbox %>% 
                  st_as_sf(coords = c(\"x\", \"y\")) %>%
                  st_set_crs(4326) %>% #current CRS is WSG84
                  st_transform(3035) %>% #transform CRS to 3035 (Lambert)
                  mutate(x_laea = unlist(map(geometry, 1)),
                         y_laea = unlist(map(geometry, 2))) %>%
                  st_set_geometry(NULL))

#Extract bathymetry for area of interest
nea <- fortify.bathy(getNOAA.bathy(lon1 = min(bbox$x), 
                                   lon2 = max(bbox$x),
                                   lat1 = min(bbox$y), 
                                   lat2 = max(bbox$y), 
                                   resolution = 5))

这很好用：

ggplot() +
  geom_sf(data = ne_countries(scale = \"medium\", 
                              returnclass = \"sf\")) +
  geom_contour(data = nea,
               aes(x = x, 
                   y = y, 
                   z = z),
               breaks = c(-600)) +
  coord_sf(xlim = c(min(bbox$x), 
                    max(bbox$x)),
           ylim = c(min(bbox$y), 
                    max(bbox$y)))

但这没有（只显示国家层，而不是 nea 一个）：

ggplot() +
  geom_sf(data = ne_countries(scale = \"medium\", 
                              returnclass = \"sf\")) +
  geom_contour(data = nea,
               aes(x = x, 
                   y = y, 
                   z = z)) +
  coord_sf(crs = 3035,
           xlim = c(min(bbox$x_laea), 
                    max(bbox$x_laea)),
           ylim = c(min(bbox$y_laea), 
                    max(bbox$y_laea)))

即使我首先将nae 转换为 LAEA：

nea_laea <- nea %>% 
  st_as_sf(coords = c(\"x\", \"y\")) %>%
  st_set_crs(4326) %>% #current CRS is WSG84
  st_transform(3035) %>% #transform CRS to 3035 (Lambert)
  mutate(x = unlist(map(geometry, 1)),
         y = unlist(map(geometry, 2))) %>%
  st_set_geometry(NULL)
ggplot() +
  geom_sf(data = ne_countries(scale = \"medium\", 
                              returnclass = \"sf\")) +
  geom_contour(data = nea_laea,
               aes(x = x, 
                   y = y, 
                   z = z)) +
  coord_sf(crs = 3035,
           xlim = c(min(bbox$x_laea), 
                    max(bbox$x_laea)),
           ylim = c(min(bbox$y_laea), 
                    max(bbox$y_laea)))

我一直在寻找解决方案，并认为一个不错的方法是从非重新投影的轮廓中提取基础数据，然后将它们与 LAEA 投影一起绘制为 geom_line：

extracted_data <- ggplot_build(ggplot() +
                                 geom_contour(data = nea,
                                              aes(x = x, y = y, z = z),
                                              breaks = c(-600)))$data[[1]] %>% 
  st_as_sf(coords = c(\"x\", \"y\")) %>%
  st_set_crs(4326) %>% #current CRS is WSG84
  st_transform(3035) %>% #transform CRS to 3035 (Lambert)
  mutate(x = unlist(map(geometry, 1)),
         y = unlist(map(geometry, 2))) %>%
  st_set_geometry(NULL)

这个几乎工作，除了似乎点没有很好地排序，所以结果图都被打乱了：

ggplot() +
  geom_sf(data = ne_countries(scale = \"medium\", 
                              returnclass = \"sf\")) +
  geom_line(data = extracted_data,
            aes(x = x, 
                y = y, 
                group = group)) +
  coord_sf(crs = 3035,
           xlim = c(min(bbox$x_laea), 
                    max(bbox$x_laea)),
           ylim = c(min(bbox$y_laea), 
                    max(bbox$y_laea)))

知道如何解决这个问题吗？

非常感谢！

“这不[工作得很好]”在这里是什么意思？错误信息？你能显示它产生的输出吗？
我编辑了我的帖子。我的意思是测深层没有显示出来，考虑到 Spacedman 在下面的帖子，这似乎是有道理的......

标签： r tidyverse sf

【解决方案1】：

一旦你将你的点投影到常规的经纬度网格以外的东西上，你就不能使用geom_contour 做一个轮廓而不重新插值到一个常规的网格上。 help(geom_contour) 解释：

Description:

     ggplot2 can not draw true 3D surfaces, but you can use
     ‘geom_contour()’, ‘geom_contour_filled()’, and ‘geom_tile()’ to
     visualise 3D surfaces in 2D. To specify a valid surface, the data
     must contain ‘x’, ‘y’, and ‘z’ coordinates, and each unique
     combination of ‘x’ and ‘y’ can appear at most once. Contouring
     requires that the points can be rearranged so that the ‘z’ values
     form a matrix, with rows corresponding to unique ‘x’ values, and
     columns corresponding to unique ‘y’ values. Missing entries are
     allowed, but contouring will only be done on cells of the grid
     with all four ‘z’ values present. If your data is irregular, you
     can interpolate to a grid before visualising using the
     ‘interp::interp()’ function from the ‘interp’ package (or one of
     the interpolating functions from the ‘akima’ package.)

我很难理解您的代码在某些地方试图做什么——您可以通过删除一些 tidyverse 的东西来使其更整洁。这个：

> bbox %>% select(y_laea) %>% filter(y_laea == min(y_laea)) %>% pull()
[1] 2904046

相当于：

> min(bbox$y_laea)
[1] 2904046

对于 max 和 X 坐标也是如此。

【讨论】：