I am trying to take a raster of soils data for one state, crop it by county, change the cell values in each county (to the county fips code), and then re-merge the county rasters back into the state raster.
Here I read in the state soils raster (which by default as the map unit key associated with each soil type as the cell value) and a polygon of US counties. I then select out the polygon just of one state, transform it to the same coordiante system as the soils raster, and then select out the soils rasters and polygons two example counties.
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)
If I plot each county by itself, you can see that the extent of the cropped region is rectangular and extends beyond the area of the county itself. The coloring is crazy because the mukey values are all over the place (although typically grouped by county). County1 lies just to the north of County2.
plot(county1_raster)
plot(county1_polygon, add = T)
plot(county2_raster)
plot(county2_polygon, add = T)
If I leave the values as is and merge the two county rasters back together, everything is fine. Even though the extents of the two rasters do overlap, the cell values are identical regardless of which raster merge is pulling from. I'm not actually sure which raster merge pulls from in this case, but it doesn't really matter. Everything fits back together nicely and the cell values are correct.
both_counties_raster <- merge(county1_raster, county2_raster)
plot(both_counties_raster)
plot(county1_polygon, add = T)
plot(county2_polygon, add = T)
However, what I want to do is to change the cell values by county prior to recombining the county rasters.
values(county1_raster) <- 1
values(county2_raster) <- 2
both_counties_raster_new <- merge(county1_raster, county2_raster)
Everything merges just fine, but when I now plot the new combined raster it is clear that for cells that were contained in both county rasters merge just took the cell values from one of the rasters. Clearly merge prioritizes the the first input raster by default.
plot(both_counties_raster_new)
plot(county1_polygon, add = T)
plot(county2_polygon, add = T)
What I'm looking for is to just change the cell values within the boundaries of each county and then merge all the counties back together again.
I am aware of the raster::mask function that can turn anything outside of the county boundary to NA, with a 10m cell resolution (described here), this takes an insane amount of time!
I have also tried an alternative approach using the raster::rasterize function to turn the county boundary polygons into a raster with the same cell size and extent of the state soils raster. Again, with a 10m cell resolution this takes forever. I was able to process one county on each of my 8 cores in 1.5 hours. And I've got a whole country to do!
I am not aware of any 10m raster US county dataset, although that would be amazing if someone pointed me to that.
The soils data is gSSURGO data - I'm also not aware if gSSURGO has within its many tables a county attribute. If it's there, I can't find it. that would also be an easy solution.
