floating forest with stars & sf

gistfile1.txt

# from: https://github.com/Irosenthal/Floating_Forests/blob/master/scripts/demo_of_merge/demo_merge.R
# see https://twitter.com/jebyrnes/status/963279659868749825
# alt: indicates a stars+sf workflow alternative to sp+raster
###########
# Libraries
###########

library(sp)
library(rgdal)
library(rgeos)
library(raster)
library(lubridate)
library(tidyverse)
library(spdplyr)

a_spdf <- readRDS("./subject_demo.rd")

# alt:
library(sf)
a_sf <- st_as_sf(a_spdf)

#----------------
#For one subject
#----------------

# Make a blank raster layer
combnR <- raster(crs=a_spdf@proj4string, ext=extent(a_spdf))

# alt:
library(stars)
combnR_stars <- st_stars(st_bbox(a_sf))

# Get the consensus raster
arast <- rasterize(SpatialPolygons(a_spdf@polygons), combnR, 
                            fun="count")  

# alt:
a <- st_as_stars(st_geometry(a_sf)) %>% # counts; use mutate for 0 -> NA conv
 mutate(values = if_else(values == 0, NA_integer_, values))

# What are the user thresholds possible in this subject?
thresholds <- na.exclude(unique(raster::values(arast)))

# alt:
thresholds_sf <- na.exclude(unique(as.vector(a[[1]])))

# Make a list of spdfs, 1 for each threshold
# This takes... a long time.
spdf_list <- lapply(thresholds, function(x){
  temp_rast <- arast
  raster::values(temp_rast)[raster::values(temp_rast)<x] <- NA
  raster::values(temp_rast)[!is.na(raster::values(temp_rast))] <- x
  rasterToPolygons(temp_rast, dissolve = TRUE)
})

# st_as_sf+st_union polygonizes, and is still much slower than raster::rasterToPolygons: W.I.P.!
sf_list <- lapply(thresholds_sf, function(x) {
  a %>%
  mutate(values = if_else(values < x, NA_integer_, x)) %>%
  st_as_sf(as_points = FALSE) %>%
  st_union
})

# Now make it all into one spdf
merged_spdf <- do.call(rbind, spdf_list) %>%
  arrange(layer)

# alt:
## st_combine merges single POLYGONs into MULTIPOLYGON, for each layer
## c combines the 9 MULTIPOLYGON into an sfc, combining layers
## st_sf makes it a single spatial data.frame
merged_sf <- lapply(sf_list, st_combine) %>% 
	do.call(c, .) %>%
	st_sf(data.frame(layer = 1:9))

# Plot them
par(mfrow=c(3,4))
arast2 <- crop(arast, extent(merged_spdf))
plot(arast2)
for(i in 1:10)
  plot(merged_spdf%>%filter(layer==i), main=paste0("Threshold = ", i)) #not sure what is up with 9 and 10...
  # Answer: control xlim and ylim, see below;
par(mfrow=c(1,1))

# alt: plot them
par(mfrow=c(3,4))
xlim = c(st_bbox(merged_sf)[c("xmin", "xmax")])
ylim = c(st_bbox(merged_sf)[c("ymin", "ymax")])
image(a, xlim = xlim, ylim = ylim) # crops
for (i in 1:9)
  merged_sf %>% filter(layer==i) %>% st_geometry %>% 
  plot(main=paste0("Threshold = ", i), xlim = xlim, ylim = ylim) 
par(mfrow=c(1,1))