AustralianAntarcticDivision / SOmap

Southern Ocean round maps
https://australianantarcticdivision.github.io/SOmap/
24 stars 6 forks source link

raster to polygons with isoband #16

Closed mdsumner closed 4 years ago

mdsumner commented 5 years ago

stuffing around with the new isoband package

Function to take a raster and create polygons .

isoband_raster <- function(x, lo, hi, auto = FALSE) {
  if (auto) {
    breaks <- pretty(values(x))
    lo <- head(breaks, -1)
    hi <- tail(breaks, -1)
  }
  ## OMG: note the [[1]] to avoid the case of as.matrix(brick) which is not helpful ...
  b <- isoband::isobands(xFromCol(x), yFromRow(x), as.matrix(x[[1]]), levels_low = lo, levels_hi = hi)
  sf::st_sf(lo = lo, hi = hi, geometry  = sf::st_sfc(isoband::iso_to_sfg(b), crs = projection(x)))
}

For the sea ice I have to smooth it, with focal() but sst and other fields work without problems.

ice <- focal(raadtools::readice()[[1]], matrix(1, 3, 3), fun = mean)
nbrks <- 20
bks <- seq(cellStats(ice, "min"), cellStats(ice, "max"), length.out = nbrks)
polys <- isoband_raster(d, head(bks, -1), tail(bks, -1))

To create bathymetry polygons, I first reduce the resolution (because faster)

d <- aggregate(SOmap::Bathy, fact = 4)
bks <- c(cellStats(d, "min"), c(-5000, -4000, -3000, -2500, -2000, -1500, -1000, -500, 0))
polys <- isoband_raster(d, head(bks, -1), tail(bks, -1))
ramp2 <- grDevices::colorRampPalette(c("#54A3D1","#60B3EB","#78C8F0","#98D1F5","#B5DCFF","#BDE1F0","#CDEBFA","#D6EFFF","#EBFAFF","grey92","grey94","grey96", "white"))
bluepal<-ramp2(nrow(polys))
plot(polys[1], border = NA, col = bluepal)
bluepal<- head(ramp2(nrow(polys) + 3), -3)
plot(polys[1], border = NA, col = bluepal)

image

Cool eh?

sq <- seq(1, nrow(polys), by = 2)
plot(polys[sq, 1], col = bluepal[sq])

image

Maschette commented 4 years ago

Is there a simpler way to do this yet @mdsumner I refer to this issue all the time.

mdsumner commented 4 years ago

Not as far as I know - one possibility is stars::st_polygonize but I think that's not really functional.

There's no natural place for this function and I just search too ¯_(ツ)_/¯

mdsumner commented 2 years ago

this looks promising https://github.com/rCarto/mapiso