英文:
Combine sf geometry with transformed geometry
问题
我正在使用R复制xkcd的“Bad Map Projection: ABS(Longitude)”。我已经有了一个不错的开端，但我想用“正确”的方式来完成。
在这里，我使用一个世界地图，并创建一个水平镜像的副本。然后我绘制两个地图，并将视窗裁剪到正经度范围内。
library(ggplot2)library(sf)library(rnaturalearth)library(rnaturalearthdata)world <- ne_countries(scale = "small", returnclass = "sf") %>%sf::st_as_sf() %>%sf::st_set_crs(value = 4326)# 矩阵乘法反射经度world_rev <- st_geometry(world) * matrix(c(-1,0,0,1), 2, 2)world2 <- world |st_set_geometry(world_rev) |sf::st_as_sf() |sf::st_set_crs(value = 4326)ggplot() +geom_sf(data = world, fill = "#336666", alpha = 0.5) +geom_sf(data = world2, fill = "#663366", alpha = 0.5) +coord_sf(xlim = c(0, 180), expand = 0)
这种方法的不足之处在于我似乎无法应用不同的投影方式。例如，如果我最后使用以下代码：
coord_sf(crs= "+proj=robin", xlim = c(0, 180), expand = 0)
我得到一个空白的图。是否有更可靠的方法来创建反转的对象？
英文:
I am replicating xkcd's "Bad Map Projection: ABS(Longitude)" with R. I have a decent start but I'd like to do it the "right" ways.
Here, I take a world map and make a copy with a horizontal reflection. Then I plot both and crop the viewing window to positive longitude.
library(ggplot2)library(sf)library(rnaturalearth)library(rnaturalearthdata)world <- ne_countries(scale = "small", returnclass = "sf") |>sf::st_as_sf() %>%sf::st_set_crs(value = 4326)# matrix multiplication to reflect longitudeworld_rev <- st_geometry(world) * matrix(c(-1,0,0,1), 2, 2)world2 <- world |>st_set_geometry(world_rev) |>sf::st_as_sf() |>sf::st_set_crs(value = 4326)ggplot() +geom_sf(data = world, fill = "#336666", alpha = 0.5) +geom_sf(data = world2, fill = "#663366", alpha = 0.5) +coord_sf(xlim = c(0, 180), expand = 0)
A shortfall of this is that I can't seem to apply different projections. For instance, if I end with
coord_sf(crs= "+proj=robin", xlim = c(0, 180), expand = 0)
I get a blank plot. Is there a more robust way to create the reversed object?
答案1
得分: 2
略微接近，但仍不理想：
- 将 Robinson 投影中心的经度从默认值 0 移动到 90
- 为了处理“超出边缘”的几何图形，将其裁剪到 0 ... 180 经度范围内
- 为了使其工作，禁用 s2
library(ggplot2)library(sf)#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUElibrary(rnaturalearth)library(rnaturalearthdata)world <- ne_countries(scale = "small", returnclass = "sf") |>sf::st_as_sf() %>%sf::st_set_crs(value = 4326)# 矩阵乘法反射经度world_rev <- st_geometry(world) * matrix(c(-1,0,0,1), 2, 2)world2 <- world |>st_set_geometry(world_rev) |>sf::st_as_sf() |>sf::st_set_crs(value = 4326)sf_use_s2(use_s2 = FALSE)#> Spherical geometry (s2) switched offbbox_0_180 <- st_bbox(c(xmin = 0, ymin = -90, xmax = 180, ymax = 90), crs = 4326)world_crp <- st_crop(world, bbox_0_180)world2_crp <- st_crop(world2, bbox_0_180)ggplot() +geom_sf(data = world_crp, fill = "#336666", alpha = 0.5) +geom_sf(data = world2_crp, fill = "#663366", alpha = 0.5) +coord_sf(crs= "+proj=robin +lon_0=90", expand = FALSE) +scale_x_continuous(breaks = seq(0,180,30))
<!-- -->
<sup>Created on 2023-07-27 with reprex v2.0.2</sup>
你使用
coord_sf(crs= "+proj=robin", xlim = c(0, 180), expand = 0)
得到了空白图，因为 Robinson 投影是投影的，单位是米，xlim 也适用于此，你可以尝试使用以下代码：
coord_sf(crs= "+proj=robin", xlim = c(-18e6, 18e6), expand = 0)
使用 world 和 world2 会产生类似以下的结果：
英文:
Somewhat closer, though still not ideal:
- move longitude of Robinson projection center from default 0 to 90
- to handle "over the edge" geometries, crop those to 0 ... 180 longitude
- for this to work, disable s2
library(ggplot2)library(sf)#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUElibrary(rnaturalearth)library(rnaturalearthdata)world <- ne_countries(scale = "small", returnclass = "sf") |>sf::st_as_sf() %>%sf::st_set_crs(value = 4326)# matrix multiplication to reflect longitudeworld_rev <- st_geometry(world) * matrix(c(-1,0,0,1), 2, 2)world2 <- world |>st_set_geometry(world_rev) |>sf::st_as_sf() |>sf::st_set_crs(value = 4326)sf_use_s2(use_s2 = FALSE)#> Spherical geometry (s2) switched offbbox_0_180 <- st_bbox(c(xmin = 0, ymin = -90, xmax = 180, ymax = 90), crs = 4326)world_crp <- st_crop(world, bbox_0_180)world2_crp <- st_crop(world2, bbox_0_180)ggplot() +geom_sf(data = world_crp, fill = "#336666", alpha = 0.5) +geom_sf(data = world2_crp, fill = "#663366", alpha = 0.5) +coord_sf(crs= "+proj=robin +lon_0=90", expand = FALSE) +scale_x_continuous(breaks = seq(0,180,30))
<!-- -->
<sup>Created on 2023-07-27 with reprex v2.0.2</sup>
You got blank plot with
coord_sf(crs= "+proj=robin", xlim = c(0, 180), expand = 0)
because Robinson is projected and units are meters, applies for xlim too, you could try this instead:
coord_sf(crs= "+proj=robin", xlim = c(-18e6, 18e6), expand = 0)
, with world and world2 would produce something like:
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论