I was thrilled to be clued into the package fasterize
, which speeds up the rasterization process in R tremendously.
However, I cannot figure out why the fasterize()
function outputs two different rasters when generated from a multipolygon with internally resolved boundaries and a multipolygon without internally resolved boundaries. I have tried specifying fun = "any"
in the fasterize parameters, as well as fun = "sum"
. The former should (I believe) take any area where there's a polygon; the latter should take the sum of any overlapping polygons. In both cases, the coverage should be the same. However, the raster from the multipolygon generated from st_combine()
drops overlapping areas whereas the raster from the multipolygon generated from st_union()
keeps any/all areas.
Am I missing something? Why are there holes in the middle plot (rasterized from st_combine()
?? There's reproducible code below the plot, with sample shapefile here (of Yellowstone and Grand Teton NPs -- provided as separate shapefile components and as zipped folder).
# Read in shapefile.
YNPGT <- st_read("YNPGT.shp")
bbox <- st_bbox(YNPGT) # Keep bounding box for plotting b/c global extent
# Rasterize. "Sum" so any spot with poly is retained (e.g., not just 'first' poly).
YNPGT.r <- YNPGT %>% fasterize(raster = ref.ras, fun = "sum")
# Combine that to flatten all polygons into one -- resolves internal boundaries
combo <- st_combine(YNPGT)
# Rasterize
combo <- st_sf(combo) # combine outputs sfc; need sf for fasterize.
combo.r <- combo %>% fasterize(raster = ref.ras, fun = "sum")
# Union that to flatten all polygons into one -- drops internal boundaries.
union <- st_union(YNPGT)
# Rasterize
union <- st_sf(union) # combine outputs sfc; need sf for fasterize.
union.r <- union %>% fasterize(raster = ref.ras, fun = "sum")
par(mfrow=c(2,3))
plot(st_geometry(YNPGT),
xlim=c(bbox[1], bbox[3]), ylim=c(bbox[2], bbox[4]),
col = c("green", "pink"))
title("YNP-GT orig")
plot(st_geometry(combo),
xlim=c(bbox[1], bbox[3]), ylim=c(bbox[2], bbox[4]),
col = c("green", "pink"))
title("YNP-GT st_combine")
plot(st_geometry(union),
xlim=c(bbox[1], bbox[3]), ylim=c(bbox[2], bbox[4]),
col = c("green", "pink"))
title("YNP-GT st_union")
plot(YNPGT.r,
xlim=c(bbox[1], bbox[3]), ylim=c(bbox[2], bbox[4]),
col = c("green", "pink"))
title("YNP-GT raster from orig")
plot(combo.r,
xlim=c(bbox[1], bbox[3]), ylim=c(bbox[2], bbox[4]),
col = c("green", "pink"))
title("YNP-GT raster from st_combine")
plot(union.r,
xlim=c(bbox[1], bbox[3]), ylim=c(bbox[2], bbox[4]),
col = c("green", "pink"))
title("YNP-GT raster from st_union")