我正在尝试使用
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)
为了避免重新采样问题,您可以在绘图之前
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))