裁剪、更改值和合并具有重叠范围的栅格

Question

我正在尝试获取一个州的土壤数据栅格，按县进行裁剪，更改每个县的单元格值（更改为县 fips 代码），然后将县栅格重新合并回州栅格.

在这里，我读取了州土壤栅格（默认情况下，它作为与每个土壤类型关联的地图单元键作为单元格值）和美国县的多边形。然后我只选择一个州的多边形，将其转换为与土壤栅格相同的坐标系统，然后选择土壤栅格和多边形两个示例县。

state_soils_raster <- raster("MapunitRaster_IL_10m.tif")
us_county_polygons <- readOGR("cb_2016_us_county_500k/cb_2016_us_county_500k.shp")

IL_county_polygons <- us_county_polygons[us_county_polygons$STATEFP == 17,]
IL_county_polygons  <- spTransform(IL_county_polygons, CRS = crs(state_soils_raster))

county1 <- "Douglas"
county2 <- "Coles"

county1_polygon <- IL_county_polygons[IL_county_polygons$NAME %in% county1,]
county2_polygon <- IL_county_polygons[IL_county_polygons$NAME %in% county2,]

county1_raster <- crop(state_soils_raster, county1_polygon)
county2_raster <- crop(state_soils_raster, county2_polygon)

如果我单独绘制每个县，您可以看到裁剪区域的范围是矩形的，并且超出了县本身的区域。着色很疯狂，因为 mukey 值无处不在（尽管通常按县分组）。 County1 位于 County2 的北部。

plot(county1_raster)
plot(county1_polygon, add = T)

plot(county2_raster)
plot(county2_polygon, add = T)

如果我保留这些值并将两个县栅格重新合并在一起，一切都很好。即使两个栅格的范围确实重叠，但无论 merge 是从哪个栅格中提取的，像元值都是相同的。在这种情况下，我实际上不确定merge 是从哪个光栅中提取的，但这并不重要。一切都很好地重新组合在一起，并且单元格值正确。

both_counties_raster <- merge(county1_raster, county2_raster)
plot(both_counties_raster)
plot(county1_polygon, add = T)
plot(county2_polygon, add = T)

但是，我想做的是在重新组合县栅格之前按县更改像元值。

values(county1_raster) <- 1
values(county2_raster) <- 2
both_counties_raster_new <- merge(county1_raster, county2_raster)

一切都合并得很好，但是当我现在绘制新的组合栅格时，很明显，对于包含在两个县栅格 merge 中的像元，只是从其中一个栅格中获取像元值。显然merge 默认优先考虑第一个输入栅格。

plot(both_counties_raster_new)
plot(county1_polygon, add = T)
plot(county2_polygon, add = T)

我正在寻找的只是更改每个县边界内的单元格值，然后将所有县重新合并在一起。

我知道raster::mask 功能可以将县界以外的任何东西变成NA，分辨率为10m（描述为here），这需要大量时间！

我还尝试了另一种方法，使用 raster::rasterize 函数将县边界多边形转换为与州土壤栅格具有相同像元大小和范围的栅格。同样，对于 10m 的单元分辨率，这需要很长时间。我能够在 1.5 小时内在我的 8 个核心上处理一个县。而且我有整个国家的事情要做！

我不知道有任何 10m 栅格美国县数据集，但如果有人指出我会很神奇。

土壤数据是 gSSURGO 数据 - 我也不知道 gSSURGO 在其许多表中是否包含县属性。如果它在那里，我找不到它。这也是一个简单的解决方案。

Answer 1

它可能不会更快，但您尝试过raster::cellFromPolygon 吗？
这是一个简单的例子：

# Create a raster with zero values
r <- raster(ncols=30, nrows=30, res = 1/3)
values(r) <- 0
# Create polygons
cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
pols <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1), Polygons(list(Polygon(cds2)), 2)))
plot(r)
plot(pols, add = TRUE)

r2 <- r
# Find which cells are in which polygons
cellpol <- cellFromPolygon(r, pols)
# Not a really clean way to attribute values in the global environment...
lapply(1:length(cellpol), function(x) values(r2)[cellpol[[x]]] <<- x)
plot(r2)
plot(pols, add = TRUE)