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

Question

我正在转换一些区域形状文件数据到光栅使用光栅功能在光栅和地球软件包。

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

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

## download shapefile of NL coropgebied regions
geo_nam <- "coropgebied"

## define year
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)
  
  # 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。此外，在您的问题中，您提到了shapefiles，但是在您的示例中没有shapefile(一种特定的文件格式)。你有一个多边形的sf (简单功能)对象。我用了一个多边形的SpatVector。