R rasterize 为不同的形状文件提供相同的值答案

【问题标题】：R rasterize giving same values for different shape filesR rasterize 为不同的形状文件提供相同的值
【发布时间】：2022-08-18 14:46:21
【问题描述】：

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

但是，我发现所有栅格最终都具有与第一个数据集相同的值。下面是一个使用两个相同 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
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)

标签： r-raster rasterizing

【解决方案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。

【讨论】：