使用 ggplot 进行地图投影

Question

我使用 Chatgpt 根据以下内容绘制了一个图： 使用 R，我想创建一个包含 50x50 网格单元的数组，其东西间距（x 坐标）为 10 度，南北间距（y 坐标）为 10 度，从西经 20 度延伸到东经 40 度，从北纬 30 度延伸到 70 度北纬度。我想要一个用指示值的颜色填充网格单元的图。我希望将其绘制在显示国家和沿海边界的地图上。该图可以保存为 png 文件。我希望颜色是透明的（alpha=0.4）。 我收到的以下代码运行得相当不错，但是当要求另一个不会将北部地区放大太多的地图投影（兰伯特）时，它完全失败了。是否不可能使用这种绘图方式，例如兰伯特项目？

# Install and load necessary packages. Are already installed
#install.packages(c("ggplot2", "sf", "rnaturalearth", "rnaturalearthdata"))
library(ggplot2)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)

# Define the grid parameters
x_range <- seq(-20, 40, length.out = 50)
y_range <- seq(30, 70, length.out = 50)

# Create a data frame with grid coordinates
df <- expand.grid(lon = x_range, lat = y_range)
df$value <- runif(nrow(df), min = 0, max = 100)  # Assign random values to each grid cell

# Prepare the world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Plot the data
plot <- ggplot() +
   geom_tile(data = df, aes(x = lon, y = lat, fill = value), alpha = 0.4) +    
scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 50, 
                   name = "Value", guide = guide_colorbar(barwidth = 10, barheight = 0.5)) +
  geom_sf(data = world, fill = NA, color = "black") +  # Add country boundaries
  coord_sf(xlim = c(-20, 40), ylim = c(30, 70), expand = FALSE) +
  theme_minimal(base_family = "sans") +
  theme(
    plot.background = element_rect(fill = "white"),
    text = element_text(color = "black"),
    axis.title = element_blank(),  # Remove axis titles
    axis.text = element_text(color = "black"),
    legend.position = "bottom"  # Position the legend at the bottom
  ) +
  labs(title = "Geographical Map with Grid Cell Values")

 # Save the plot as a PNG file with custom width and height
 ggsave("map_plot.png", plot = plot, width = 10, height = 8, dpi = 300)

Answer 1

问题在于

geom_tile

使用表示度数的原始数字作为其坐标。当您在投影之间进行转换时，坐标系不再以度为单位，因此

geom_tile

不再与地图对齐。 处理这个问题的最全面的方法是将图块创建为

sf

数据框。这需要一些计算，但看起来像这样：

x_range <- seq(-20, 40, length.out = 50)
y_range <- seq(30, 70, length.out = 50)

df <- expand.grid(lon = x_range, lat = y_range)

df <- lapply(seq(nrow(df) - 1), function(i) {
  if(df$lon[i + 1] < df$lon[i]) return(NULL)
  st_polygon(list(cbind(
    c(seq(df$lon[i], df$lon[i + 1], length = 10),
      rep(df$lon[i + 1], 10),
      seq(df$lon[i + 1], df$lon[i], length = 10),
      rep(df$lon[i], 10)),
    c(rep(df$lat[i], 10),
      seq(df$lat[i], df$lat[i] + 40/49, length = 10),
      rep(df$lat[i] + 40/49, 10),
      seq(df$lat[i] + 40/49, df$lat[i], length = 10))
      )
  ))}) |> 
  st_sfc(crs = 'WGS84') |>
  st_sf() |>
  setNames('geometry') |>
  `st_geometry<-`("geometry")
  
df$value <- runif(nrow(df), min = 0, max = 100)

现在我们可以通过将世界地图裁剪为预先指定的“矩形”来创建具有我们想要的精确区域的绘图：

world <- ne_countries(scale = "medium", returnclass = "sf")

sf_use_s2(FALSE)

plot <- st_polygon(x = list(cbind(c(seq(-20, 40, length = 100),
                            rep(40, 100),
                            seq(40, -20, length = 100),
                            rep(-20, 100)), 
                          c(rep(30, 100),
                            seq(30, 70, length = 100),
                            rep(70, 100), seq(70, 30, length = 100))))) |>
  st_sfc(crs = 'WGS84') |>
  st_intersection(world) |>
  ggplot() +
  geom_sf(data = df, aes(fill = value), alpha = 0.4, color = NA) +    
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 50, 
                       name = "Value", guide = guide_colorbar(barwidth = 10, barheight = 0.5)) +
  geom_sf(fill = NA, color = "black", linewidth = 0.3) + 
  theme_minimal(base_family = "sans") +
  theme(
    plot.background = element_rect(fill = "white"),
    text = element_text(color = "black"),
    axis.title = element_blank(), 
    axis.text = element_text(color = "black"),
    legend.position = "bottom"  
  ) +
  labs(title = "Geographical Map with Grid Cell Values")

所有这些只是给我们提供了与最初相同的情节：

plot

除了我们现在可以使用

coord_sf

将其转换为我们喜欢的任何投影。例如，我们得到一个很好的兰伯特投影

crs = 3416

:

plot + coord_sf(expand = FALSE, crs = 3416) 

或者像这样的 Mollweide 投影：

plot + coord_sf(expand = FALSE, crs = "+proj=moll")