R 光栅化为不同形状文件提供相同的值

Question

我正在使用 raster 和 terra 包中的 rasterize 函数将一些区域 shapefile 数据转换为栅格。

但是，我发现所有栅格最终都具有与第一个数据集相同的值。下面是使用两个相同的空间矢量多边形的示例，一个的数据在 1 到 100 之间，另一个的数据在 1000 到 2000 之间。即使我重新启动 R 并仅运行使用 geo_dat2 (1000-2000) 生成的 rast_d2，我也会得到 geo_dat1 的输出（1-100）。是否有某种缓存存储在我需要清除的地方？

library(raster)
library(httr)
library(sf)
library(dplyr)
library(mapview)

## download shapefile of NL coropgebied regions
## some geo options: coropgebied, gemeente, provincie, lansdeel
geo_nam <- "coropgebied"

url <- parse_url("https://service.pdok.nl/cbs/gebiedsindelingen/2023/wfs/v1_0")

url$query <- list(service = "WFS",
                  version = "2.0.0",
                  request = "GetFeature",
                  typename = paste0("gebiedsindelingen:", geo_nam, "_gegeneraliseerd"),
                  outputFormat = "application/json")
request <- build_url(url)
  
  # request sf and transform to 4326
  geo_sf <- st_read(request) %>% 
    st_transform(4326)
  
  ## generate desired raster
  r <- raster(nrows=122, ncols=87, xmn=3.15, xmx=7.5, ymn=50.65, ymx=53.7, 
              crs = 4326)
  
  ## crop to extent
  r_crop <- crop(r, geo_sf)

  ## generate some random data for the regions that is significantly different
  ## dat1 is 1-100
  geo_dat1 <- geo_sf %>% 
    mutate(dat = as.numeric(sample(1:100, NROW(geo_sf$statcode)))) %>% 
    select(dat, geometry)
  
  ## dat2 is 1000:2000
  geo_dat2 <- geo_sf %>% 
    mutate(dat = as.numeric(sample(1000:2000, NROW(geo_sf$statcode)))) %>% 
    select(dat, geometry)
  
  ## use raster to rasterize the shape file data
  rast_d1 <- terra::rasterize(geo_dat1, r_crop)
  crs(rast_d1) <- 4326
  
  rast_d2 <- terra::rasterize(geo_dat2, r_crop)
  crs(rast_d2) <- 4326
  
  ## plot both
  plot(rast_d1)
  plot(rast_d2)
  
  ## extract values from raster
  rast_d1@data@max
  rast_d2@data@max
 
  ## however the scale when plotting with mapview is consistent with expected range, the cell values are not.
mapview(rast_d1)  
mapview(rast_d2)

Answer 1

您需要指定要栅格化的变量。我在下面展示了一个使用“terra”的简化脚本。

library(terra)
library(httr)

geo_nam <- "coropgebied"
year <- "2021"
url <- parse_url("https://geodata.nationaalgeoregister.nl/cbsgebiedsindelingen/wfs")
url$query <- list(service = "WFS", version = "2.0.0", request = "GetFeature",
             typename = paste0("cbsgebiedsindelingen:cbs_", geo_nam, "_", year, "_gegeneraliseerd"),
             outputFormat = "application/json")
request <- build_url(url)
  
geo <- vect(request) |> project("EPSG:4326")      
r <- rast(geo, res=0.025) 

geo$dat1 <- sample(100, nrow(geo))
geo$dat2 <- 1000 + sample(100, nrow(geo))
  
rd1 <- rasterize(geo, r, "dat1")
rd2 <- rasterize(geo, r, "dat2")

请注意，在调用

terra::rasterize

时，您提供了

RasterLayer

参数，在这种情况下，您使用的方法将是

raster::rasterize

。另外，在您的问题中，您引用了 shapefile，但您的示例中没有 shapefile（一种特定的文件格式）。您有一个多边形的 sf（简单特征）对象。我使用了多边形的 SpatVector。