将 coord_sf 与 geom_spatraster 一起使用会改变分辨率

问题描述 投票:0回答:1

我正在尝试使用

tidyterra
ggplot2
绘制大型栅格的子区域(例如,覆盖北美的栅格中的亚利桑那州)。我尝试通过使用
coord_sf
将 CRS 设置为经/纬度 WGS 84,然后使用
lim
以及经度/纬度的限制来设置绘图的限制。当我这样做时,生成的图显示的分辨率比栅格最初的分辨率(1 公里)低得多。我意识到我可以在栅格的原始 CRS 中设置限制,但我还想在经纬度 WGS 84 的顶部绘制矢量,而不是原始 CRS。有没有办法可以保留原始分辨率,同时使用纬度/经度 WGS 84 来指定限制和其他几何图形?

library(tidyverse)
library(ggplot2)
library(terra)
library(tidyterra)
library(maps)

### Example data:

# example raster:
ex <- rast(
  extent = c(-4352000, 3336000, -3341000, 4276000),
  resolution = c(1000, 1000),
  crs = "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
)
ex[] <- runif(n = ncell(ex), min = 0, max = 100)

# example polygon
az <- map_data("state") %>%
  filter(region %in% c("arizona")) %>%
  vect(
    x = cbind(
      as.factor(.$region),    # object (state) id
      .$group,                # polygon part id
      .$long,                 # x
      .$lat,                  # y
      0                             # hole flag (all 0)
    ),
    type = "polygons",
    crs = "+proj=longlat +datum=WGS84"
  )

### Plot:

p <- ggplot() +
  geom_spatraster(data = ex) +
  geom_spatvector(data = az, fill = NA, color = "white") +
  coord_sf(crs = 4326) +
  lims(x = c(-116, -108), y = c(30, 38))
print(p)

enter image description here

r ggplot2 terra tidyterra
1个回答
0
投票

为了避免重新采样问题,您可以在绘图之前

terra::crop()
您的数据。鉴于您想要在 WGS84/EPSG:4326 中绘制数据,您必须在裁剪之前投影您的 spatraster。这样做就不需要在
coords_sf()
中定义 CRS,因为所有数据已经具有相同的 CRS。

library(terra)
library(tidyterra)
library(sf)
library(maps)
library(ggplot2)

# Example raster
set.seed(1)
ex <- rast(
  extent = c(-4352000, 3336000, -3341000, 4276000),
  resolution = c(1000, 1000),
  crs = "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
) %>%
  project("EPSG:4326")

ex[] <- runif(n = ncell(ex), min = 0, max = 100)

# Create sf for cropping raster using your extent values
crop_sf <- data.frame(lon = c(-116, -116, -108, -108),
                      lat = c(30, 38, 38, 30)) %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

crop_sf <- do.call(c, st_geometry(crop_sf)) %>%
  st_cast("POLYGON") %>%
  st_sfc() %>%
  st_as_sf(crs = 4326)

# Crop ex spatraster
ex1 <- crop(ex, crop_sf)

# Example polygon
az <- map_data("state") %>%
  filter(region %in% c("arizona")) %>%
  vect(
    x = cbind(
      as.factor(.$region),    # object (state) id
      .$group,                # polygon part id
      .$long,                 # x
      .$lat,                  # y
      0                       # hole flag (all 0)
    ),
    type = "polygons",
    crs = "EPSG:4326"
  )

# Plot
ggplot() +
  geom_spatraster(data = ex1) +
  geom_spatvector(data = az, fill = NA, color = "white") +
  coord_sf(xlim = c(-116, -108), ylim = c(30, 38))

result

© www.soinside.com 2019 - 2024. All rights reserved.