Why does fasterize gives different output for multipolygons with and without internally resolved boundaries?

Question

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")

Answer 1

There's a couple of things going on here, I think. First, fasterize doesn't have ... enabled so you can't pass na.rm = TRUE through your summarising function, and all of those default to FALSE.

The other issue is that as the documentation warns, the output of st_combine() is often not valid under the simple features spec so may perform unexpectedly. If you run combo through lwgeom::st_make_valid() %>% st_collection_extract(., 'POLYGON'), you'll see the combination resolves to a valid sfc_MULTIPOLYGON but treats some overlapping polygons from YNPGT[2, ] as holes, because they're fully contained by YNPGT[1, ].

I don't think the results of st_combine() with polygons are generally safe for rasterizing. I get the sense that its most useful for point and line geometries.