isoband
Generate Isolines and Isobands from Regularly Spaced Elevation Grids
A fast C++ implementation to generate contour lines (isolines) and contour polygons (isobands) from regularly spaced grids containing elevation data.
Generate contour lines (isolines) and contour polygons (isobands) from regularly spaced grids containing elevation data.
Installation
Install from github with:
devtools::install_github("clauswilke/isoband")
Examples
The two main workhorses of the package are the functions isolines()
and isobands(), respectively. They return a list of isolines/isobands
for each isolevel specified. Each isoline/isoband consists of vectors of
x and y coordinates, as well as a vector of ids specifying which sets of
coordinates should be connected. This format can be handed directly to
grid.polyline()/grid.path() for drawing. However, we can also
convert the output to spatial features and draw with ggplot2 (see
below).
library(isoband)m <- matrix(c(0, 0, 0, 0, 0,0, 1, 2, 1, 0,0, 1, 2, 0, 0,0, 1, 0, 1, 0,0, 0, 0, 0, 0), 5, 5, byrow = TRUE)isolines(1:ncol(m), 1:nrow(m), m, 0.5)#> $`0.5`#> $`0.5`$x#> [1] 4.00 3.50 3.00 2.50 2.00 1.50 1.50 1.50 2.00 3.00 4.00 4.50 4.00 3.75#> [15] 4.00 4.50 4.00#>#> $`0.5`$y#> [1] 4.50 4.00 3.75 4.00 4.50 4.00 3.00 2.00 1.50 1.25 1.50 2.00 2.50 3.00#> [15] 3.50 4.00 4.50#>#> $`0.5`$id#> [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1#>#>#> attr(,"class")#> [1] "isolines" "iso"isobands(1:ncol(m), 1:nrow(m), m, 0.5, 1.5)#> $`0.5:1.5`#> $`0.5:1.5`$x#> [1] 2.50 2.00 1.50 1.50 1.50 2.00 3.00 4.00 4.50 4.00 3.75 4.00 4.50 4.00#> [15] 3.50 3.00 3.00 3.25 3.50 3.00 2.50 2.50#>#> $`0.5:1.5`$y#> [1] 4.00 4.50 4.00 3.00 2.00 1.50 1.25 1.50 2.00 2.50 3.00 3.50 4.00 4.50#> [15] 4.00 3.75 3.25 3.00 2.00 1.75 2.00 3.00#>#> $`0.5:1.5`$id#> [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2#>#>#> attr(,"class")#> [1] "isobands" "iso"
The function plot_iso() is a convenience function for debugging and
testing.
plot_iso(m, 0.5, 1.5)
A few more simple examples. Missing values and disconnected areas are all handled correctly.
m <- matrix(c(0, 0, 0, 0, 0, 0,0, 0, 0, 1, 1, 0,0, 0, 1, 1, 1, 0,0, 1, 1, 0, 0, 0,0, 0, 0, 1, 0, 0,0, 0, 0, 0, 0, 0), 6, 6, byrow = TRUE)plot_iso(m, 0.5, 1.5)
m <- matrix(c(NA, NA, NA, 0, 0, 0,NA, NA, NA, 1, 1, 0,0, 0, 1, 1, 1, 0,0, 1, 1, 0, 0, 0,0, 0, 0, 1, 0, 0,0, 0, 0, 0, 0, 0), 6, 6, byrow = TRUE)plot_iso(m, 0.5, 1.5)
m <- matrix(c(0, 0, 1, 1,0, 1, 1, 1,1, 1, 0, 0,0, 0, 0.8, 0), 4, 4, byrow = TRUE)plot_iso(m, 0.5, 1.5)
The algorithm has no problem with larger datasets. Let’s calculate isolines and isobands for the volcano dataset, convert to sf, and plot with ggplot2.
library(ggplot2)library(sf)#> Linking to GEOS 3.6.1, GDAL 2.1.3, PROJ 4.9.3m <- volcanob <- isobands((1:ncol(m))/(ncol(m)+1), (nrow(m):1)/(nrow(m)+1), m, 10*(9:19), 10*(10:20))l <- isolines((1:ncol(m))/(ncol(m)+1), (nrow(m):1)/(nrow(m)+1), m, 10*(10:19))bands <- iso_to_sfg(b)data_bands <- st_sf(level = 1:length(bands),geometry = st_sfc(bands))lines <- iso_to_sfg(l)data_lines <- st_sf(level = 2:(length(lines)+1),geometry = st_sfc(lines))ggplot() +geom_sf(data = data_bands, aes(fill = level), color = NA, alpha = 0.7) +geom_sf(data = data_lines, color = "black") +scale_fill_viridis_c(guide = "none") +coord_sf(expand = FALSE)
We can also use this approach to convert bitmap images into polygons and plot with ggplot2.
library(magick)#> Linking to ImageMagick 6.9.9.39#> Enabled features: cairo, fontconfig, freetype, lcms, pango, rsvg, webp#> Disabled features: fftw, ghostscript, x11sf_from_image <- function(image) {image_gray <- image %>% image_quantize(colorspace = "gray")image_raster <- as.raster(image_gray)d <- dim(image_raster)m <- matrix(c((255-col2rgb(image_raster)[1,])), nrow = d[1], ncol = d[2], byrow = TRUE)b <- isobands(1:d[2], d[1]:1, m, 20*(0:13), 20*(1:14))bands <- iso_to_sfg(b)data <- st_sf(level = letters[1:length(bands)],geometry = st_sfc(bands))}img <- image_resize(image_read("https://pbs.twimg.com/profile_images/905186381995147264/7zKAG5sY_400x400.jpg"), "200x200")img_sf <- sf_from_image(img)p1 <- ggplot(img_sf) +geom_sf(color = "gray10", fill = NA, size = 0.05) +coord_sf(expand = FALSE) +theme_gray() +theme(axis.ticks = element_blank(),axis.text = element_blank(),axis.title = element_blank(),axis.ticks.length = grid::unit(0, "pt"),plot.margin = margin(0, 0, 0, 0))p2 <- ggplot(img_sf) +geom_sf(aes(fill = level, color = level)) +coord_sf(expand = FALSE) +theme_void()p3 <- ggplot(img_sf) +geom_sf(aes(fill = level), color = "gray30", size = 0.1) +scale_fill_hue(aesthetics = c("color", "fill"), guide = "none", direction = -1) +coord_sf(expand = FALSE) +theme_void()cowplot::plot_grid(p1,p2 + scale_fill_viridis_d(aesthetics = c("color", "fill"), option = "B", guide = "none",direction = -1),p2 + scale_fill_viridis_d(aesthetics = c("color", "fill"), option = "D", guide = "none",direction = -1),p3,scale = 0.9)
Performance
The code is written in C++ and performance is generally good. Isolining
is about as fast as grDevices::contourLines(), isobanding is
approximately 2.5 times slower.
microbenchmark::microbenchmark(grDevices::contourLines(1:ncol(volcano), 1:nrow(volcano), volcano, levels = 10*(10:18)),isolines(1:ncol(volcano), 1:nrow(volcano), volcano, 10*(10:18)),isobands(1:ncol(volcano), 1:nrow(volcano), volcano, 10*(9:17), 10*(10:18)))#> Unit: milliseconds#> expr#> grDevices::contourLines(1:ncol(volcano), 1:nrow(volcano), volcano, levels = 10 * (10:18))#> isolines(1:ncol(volcano), 1:nrow(volcano), volcano, 10 * (10:18))#> isobands(1:ncol(volcano), 1:nrow(volcano), volcano, 10 * (9:17), 10 * (10:18))#> min lq mean median uq max neval#> 1.623163 1.754017 2.316944 1.892529 2.178831 10.729141 100#> 1.750317 1.857125 2.276085 1.947959 2.359221 9.654656 100#> 4.382390 4.625374 5.072554 4.800300 5.340838 9.739161 100