文章/答案/技术大牛

发布

社区首页 >问答首页 >带有Lambert投影的geom_contour

问带有Lambert投影的geom_contour
EN

Stack Overflow用户

提问于 2022-08-03 10:11:40

回答 2查看 127关注 0票数 0

我正在绘制一张东北大西洋的地图，上面突出了特定的等深线。考虑到这个区域，我更愿意使用LAEA投影。然而，这似乎会导致geom_contour失败。

以下是背景代码：

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

#Set boundaries
bbox <- tibble(x = c(-20, 10), y = c(45, 60))
#Add coordinates 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)))

知道怎么解决这个问题吗？

非常感谢!

回答 2

Stack Overflow用户

回答已采纳

发布于 2022-08-03 22:28:22

另一种解决方案是使用以下方法：

sf::st_transform())

Download

在没有投影的情况下设置边框(WGS84 -> crs = 4326)，并在兰伯特(使用sf::st_transform())

plot使用ggplot2和sf.

在Lambert使用测量水深数据和在兰伯特使用marmap和raster raster的国家数据和项目)进行投影(sf::st_transform())

plot使用ggplot2和sf.

)

这种方法的优点是，您的所有数据层都是在Lambert中投影的。因此，您可以使用ggplot2和sf的全部潜力进行绘图。更具体地说，您可以使用geom_tile()，并根据需要使用尽可能多的geom_contour()层。如果要绘制几条等高线，则不必单独计算:只需在breaks参数geom_contour()中添加等值线值即可。

下面是一个坐标和-600米等深线的例子：

# Load useful packages
library(sf)
library(marmap)
library(tidyverse)
library(rnaturalearth)

# Set Lambert projection
proj = 3035

# Set bounding box (WGS 84)
min_lon <- -20
max_lon <- 10
min_lat <- 45
max_lat <- 60

# Project bbox
bbox <- st_sfc(st_point(c(min_lon, min_lat)), 
               st_point(c(max_lon, max_lat)), 
               crs = 4326)
bbox_proj <- st_coordinates(st_transform(bbox, crs = proj))

# Get bathymetric data
bathy <- getNOAA.bathy(lon1 = min_lon, lon2 = max_lon,
                       lat1 = min_lat, lat2 = max_lat,
                       resolution = 5, keep = TRUE)

# Project bathymetric data and transform to plot with ggplot2
proj_bathy <- as.raster(bathy) %>% 
  raster::projectRaster(crs = proj) %>% 
  as.bathy() %>% 
  as.xyz()

# Import country data and project
country <- ne_countries(scale = "medium", returnclass = "sf")
country_proj <- st_transform(country, crs = proj)

# Plot using ggplot and sf
ggplot() + 
  geom_sf(data = country_proj) +
  geom_contour(data = proj_bathy, 
               aes(x = V1, y = V2, z = V3),
               breaks = -600, color = "red") +
  coord_sf(xlim = bbox_proj[,'X'], 
           ylim = bbox_proj[,'Y']) +
  labs(x = "Longitude", y = "Latitude") +
  theme_minimal()

这个阴谋产生了：

票数 1

Stack Overflow用户

发布于 2022-08-03 13:19:36

一旦你把你的点投射到常规的长网格之外，你就不能使用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.)

我很难理解您的代码在某些地方想要做什么--您可以通过删除一些有趣的东西来使它更整洁。这是：

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

相当于：

> min(bbox$y_laea)
[1] 2904046

与之类似的是最大和弦和X弦。

如果您只想获得一个特定的轮廓集，请使用contourLines，而不是尝试从绘图结构中获取值。下面是一个完整的例子，我尝试使它尽可能简单，所以它只使用sf和marmap

使用以下软件包：

library(sf)
library(marmap)

设置限制--我们不会转换它们，所以没有必要用它们来构建数据框架。它还简化了提取限制：

xr=c(-20, 10)
yr=c(45,60)

把水深测量作为一个矩阵--不要“强化”它：

bathy <- getNOAA.bathy(lon1 = min(xr), 
                     lon2 = max(xr),
                     lat1 = min(yr),
                     lat2 = max(yr),
                     resolution = 5)

现在从酒窝里拿出拉特和龙的和弦。有关我从哪里获得这些信息，请参见marmap::as.xyz：

lon = as.numeric(rownames(bathy))
lat = as.numeric(colnames(bathy))

现在我们有了一个矩阵和坐标，这样我们就可以将它提供给contourLines并得到-600级别：

c600_list = contourLines(x=lon, y=lat, z=bathy, levels=-600)

这是一个包含$x和$y组件的元素列表，因此我们使用lapply循环它们并构建st_linestring对象，然后使用st_sfc连接该几何学原语列表。

c600 = do.call(st_sfc,
             lapply(c600_list, function(l){st_linestring(cbind(l$x, l$y))}) 
             )

到目前为止，一切都是长时间的epsg:4326，让我们来设定：

st_crs(c600) = 4326

现在我们可以转换为LAEA epsg:3035：

c600_laea = st_transform(c600, 3035)

您可以使用ggplot2来绘制：

library(ggplot2)
ggplot() + geom_sf(data=c600_laea)

或者在交互模式下使用tmap包获取一些上下文，并确保我们在正确的位置：

library(tmap)
tmap_mode("view")
# tmap mode set to interactive viewing
tm_shape(c600_laea) + tm_lines()

票数 3

页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持

原文链接：

https://stackoverflow.com/questions/73219945

复制

相似问题

问带有Lambert投影的geom_contour
EN

回答 2

Stack Overflow用户

Stack Overflow用户

社区

活动

圈层

关于

腾讯云开发者

热门产品

热门推荐

更多推荐

问带有Lambert投影的geom_contourEN

回答 2

Stack Overflow用户

Stack Overflow用户

社区

活动

圈层

关于

腾讯云开发者

热门产品

热门推荐

更多推荐

问带有Lambert投影的geom_contour
EN