我有一個縣邊界的矢量空間資料和同一縣的柵格格式的地形資料。

我需要使用 ggplot 創建一個地圖，以便它僅以柵格格式顯示縣邊界內的那些資料（這些資料又位于矢量空間檔案中）。

換句話說，我需要洗掉縣輪廓之外的柵格資料。可以用ggplot做到這一點嗎？

可重現的例子：

# load packages

library(elevatr) 
library(terra) 
library(geobr)

# get the municipality shapefile (vectorized spatial data)

municipality_shape <- read_municipality(code_muni = 3305802)

plot(municipality_shape$geom) 



# get the raster topographical data

prj_dd <- "EPSG:4674" 

t <- elevatr::get_elev_raster(locations = municipality_shape, 
                              z = 10,
                              prj = prj_dd) 

obj_raster <- rast(t) 

plot(obj_raster) 



# create the ggplot map

df_tere_topo <- obj_raster %>%
  as.data.frame(xy = TRUE) %>%
  rename(altitude = file40ac737835de)

ggplot() 
  geom_raster(data = df_tere_topo, aes(x = x, y = y, fill = `altitude`)) 
  geom_sf(municipality_shape, mapping = aes(), color = 'red', fill = NA)

uj5u.com熱心網友回復：

已編輯

請參閱評論，使用terra::crop()andterra::mask()代替：

# load packages
library(elevatr)
library(terra)
library(geobr)
library(dplyr)
library(ggplot2)
# Use tidyterra
library(tidyterra)

# get the municipality shapefile (vectorized spatial data)

municipality_shape <- read_municipality(code_muni = 3305802)



# get the raster topographical data

prj_dd <- "EPSG:4674"

t <- elevatr::get_elev_raster(
  locations = municipality_shape,
  z = 10,
  prj = prj_dd
)

obj_raster <- rast(t)

# Crop   Mask

obj_raster_mask <- crop(obj_raster, vect(municipality_shape)) %>%
  mask(vect(municipality_shape))


# create the ggplot map
# using tidyterra
ggplot()  
  geom_spatraster(data = obj_raster_mask)  
  geom_sf(municipality_shape, mapping = aes(), color = "white", fill = NA)  
  # Not plotting NAs of the raster
  scale_fill_continuous(na.value = NA)  
  labs(fill="altitude")

原始答案

我認為最有效的方法是將crop您的 SpatRaster 擴展到矢量資料的范圍內。使用這種方法，繪圖效率更高，因為您沒有使用您不想繪制的資料。

另一種選擇是在coord_sf().

在這個 reprex 中，我也使用了tidyterra包，它具有一些與ggplot2 terra一起使用的功能（我是 tidyterra 的開發人員）：

# load packages
library(elevatr)
library(terra)
library(geobr)
library(dplyr)
library(ggplot2)
# Use tidyterra
library(tidyterra)

# get the municipality shapefile (vectorized spatial data)

municipality_shape <- read_municipality(code_muni = 3305802)



# get the raster topographical data

prj_dd <- "EPSG:4674"

t <- elevatr::get_elev_raster(
  locations = municipality_shape,
  z = 10,
  prj = prj_dd
)

obj_raster <- rast(t)

# Option1: Crop first

obj_raster_crop <- crop(obj_raster, vect(municipality_shape))



# create the ggplot map
# using tidyterra
ggplot()  
  geom_spatraster(data = obj_raster_crop)  
  geom_sf(municipality_shape, mapping = aes(), color = "white", fill = NA)  
  coord_sf(expand = FALSE)  
  labs(fill="altitude")

 # Option 2: Use limits and no crop
lims <- sf::st_bbox(municipality_shape)

ggplot()  
  geom_spatraster(
    data = obj_raster,
    # Avoid resampling
    maxcell = ncell(obj_raster)
  )  
  geom_sf(municipality_shape, mapping = aes(), color = "white", fill = NA)  
  coord_sf(
    expand = FALSE,
    xlim = lims[c(1, 3)],
    ylim = lims[c(2, 4)]
  )  
  labs(fill="altitude")

標籤：rggplot2光栅空间的科幻

上一篇：st_interpolate_aw回傳錯誤-替換有x行，資料有y

下一篇：將緩沖區與向量相交以提取資訊