2

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

comparison plots

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

  • Thanks for these thoughts -- very helpful. I'll be cautious when using st_combine() in the future. – ltlf653 Apr 22 at 13:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.