Smoothing isobands

[michelk]

Does anyone has an idea how to smooth isobands, while preserving a valid geometry, where class-boundaries share the same polygon vertices?

@mstrimas's smoothr::smooth function smooths each polygon separately. The downside of this approach is, that this leads to gaps and overlaps with neighbouring polygons (see below the image).

Would it be possible to solve this within the isoband package?

Would it be possible to generate isolines, smooth them and convert them afterwards back to isobands?

The following illustrates the issue:

library(ggplot2)
library(isoband)
library(sf)
m <- volcano
b <-
  isobands(x = (1:ncol(m))/(ncol(m)+1)
         , y = (nrow(m):1)/(nrow(m)+1)
         , z = m
         , levels_low = 10*(9:19)
         , levels_high = 10*(10:20))
bands <- iso_to_sfg(b)
volcano_contour <-
  st_sf(
    level = 1:length(bands),
    geometry = st_sfc(bands)
  )
volcano_contour_smooth <-
  smoothr::smooth(volcano_contour, method = "ksmooth", smoothness = 2)

p1 <-
  ggplot() +
  geom_sf(data = volcano_contour, aes(fill = level), color = NA, alpha = 1.0) +
  scale_fill_viridis_c(guide = "none") +
  coord_sf(expand = FALSE)
p2 <-
  ggplot() +
  geom_sf(data = volcano_contour_smooth, aes(fill = level), color = NA, alpha = 1.0) +
  scale_fill_viridis_c(guide = "none") +
  coord_sf(expand = FALSE)

p <- cowplot::plot_grid(p1, p2, labels = c('Raw', 'Smooth'), label_size = 12)

ggsave('volcano.png', p, width = 210, height = 148, units = 'mm')

[clauswilke]

I think you'll be better off smoothing the input raster data. Increase the resolution of the grid and interpolate and then run the isobanding algorithm on that input.

[michelk]

Thanks. I will try that.

[michelk]

The Disaggregation is quite tricky. Is there an algorithm you would suggest? @mstrimas uses in his smoothr::smooth function ksmooth() from base. So I guess, we need a method to apply that function to a raster object.

library(ggplot2)                                                                                                                                                                                                                                                                  
library(isoband)                                                                                                                                                                                                                                                                  
library(sf)                                                                                                                                                                                                                                                                       
library(raster)                                                                                                                                                                                                                                                                   
library(rasterVis)                                                                                                                                                                                                                                                                

toIsobandSf <-                                                                                                                                                                                                                                                                    
  function                                                                                                                                                                                                                                                                        
(                                                                                                                                                                                                                                                                                 
  r  # ^ raster::RasterLayer                                                                                                                                                                                                                                                      
)                                                                   
{                                                                   

  b <-                                                              
    isobands(x = xFromCol(r)                                        
           , y = yFromRow(r)                                        
           , z = as.matrix(r[[1]])                                  
           , levels_low = 10*(9:19)                                 
           , levels_high = 10*(10:20))                              
  bands <- iso_to_sfg(b)                                            
  volcano_contour <-                                                
    st_sf(                                                          
      level = 1:length(bands),                                      
      geometry = st_sfc(bands)                                      
    )                                                               
} # ^ sf::sf                                                        

plotContour <-                                                      
  function                                                          
(                                                                   
  dd # ^ sf::sf                                                     
)                                                                   
{                                                                   
  ggplot() +                                                        
  geom_sf(data = dd, aes(fill = level), color = NA, alpha = 1.0) +                                                                       
  scale_fill_viridis_c(guide = "none") +                            
  coord_sf(expand = FALSE)                                          
} # ^ ggplot2::gg                                                   

v <- raster(volcano)                                                

## Raw Contour                                                      
volcano_contour <- toIsobandSf(v)                                   

## Contour smoothed with smoothr                                    
volcano_contour_smooth <-                                           
  smoothr::smooth(volcano_contour, method = "ksmooth", smoothness = 5)                                                                   

## Increase resolution                                              
v_res_hi <- raster::disaggregate(v, 3, method = "bilinear")                                                                              
volcano_contour_res_hi <- toIsobandSf(v_res_hi)                     

## Blur raster                                                      
v_res_hi_blur <- focal(v_res_hi, w=matrix(1, 5, 5), mean)                                                                                
volcano_contour_res_hi_blur <- toIsobandSf(v_res_hi_blur)                                                                                

## Plot data                                                        
p1 <- plotContour(volcano_contour)                                  
p2 <- plotContour(volcano_contour_smooth)                           
p3 <- plotContour(volcano_contour_res_hi)                           
p4 <- plotContour(volcano_contour_res_hi_blur)                      

p <- cowplot::plot_grid(p1, p2, p3, p4, labels = c('Raw', 'smoothr::smooth', 'Raw Res. High', 'Raw Res. High, Raster blur'))                                                                                                                                                      

ggsave('volcano.png', p, width = 210, height = 148, units = 'mm')                                                                        

[clauswilke]

Looks to me like it's working. You'd probably have to increase the resolution or the blur even more to get smoother curves.

Btw., if you set hjust = 0 in the plot_grid() function the labels will stay in place. https://stackoverflow.com/a/47724512/4975218

[michelk]

Thanks. With Gauss, it gets a little better, although I still prefer the aesthetics of smoothr::smooth

library(ggplot2)                                                                                                                                                                                                                                                                                                               
library(isoband)                                                                                                                                                                                                                                                                                                               
library(sf)                                                                                                                                                                                                                                                                                                                    
library(raster)                                                                                                                                                                                                                                                                                                                
library(rasterVis)                                                                                                                                                                                                                                                                                                             

toIsobandSf <-                                                                                                                                                                                                                                                                                                                 
  function                                                                                                                                                                                                                                                                                                                     
(                                                                                                                                                                                                                                                                                                                              
  r  # ^ raster::RasterLayer                                                                                                                                                                                                                                                                                                   
)                                                                                                                                                                                                                                                                                                                              
{                                                                                                                                                                                                                                                                                                                              

  b <-                                                                                                                                                                                                                                                                                                                         
    isobands(x = xFromCol(r)                                                                                                                                                                                                                                                                                                   
           , y = yFromRow(r)                                                                                                                                                                                                                                                                                                   
           , z = as.matrix(r[[1]])                                                                                                                                                                                                                                                                                             
           , levels_low = 10*(9:19)                                                                                                                                                                                                                                                                                            
           , levels_high = 10*(10:20))                                                                                                                                                                                                                                                                                         
  bands <- iso_to_sfg(b)                                                                                                                                                                                                                                                                                                       
  volcano_contour <-                                                                                                                                                                                                                                                                                                           
    st_sf(                                                                                                                                                                                                                                                                                                                     
      level = 1:length(bands),                                                                                                                                                                                                                                                                                                 
      geometry = st_sfc(bands)                                                 
    )                                                                          
} # ^ sf::sf                                                                   

plotContour <-                                                                 
  function                                                                     
(                                                                              
  dd # ^ sf::sf                                                                
)                                                                              
{                                                                              
  ggplot() +                                                                   
  geom_sf(data = dd, aes(fill = level), color = NA, alpha = 1.0) +                                                                                             
  scale_fill_viridis_c(guide = "none") +                                       
  coord_sf(expand = FALSE)                                                     
} # ^ ggplot2::gg                                                              

v <- raster(volcano)                                                           

## Raw Contour                                                                 
volcano_contour <- toIsobandSf(v)                                              

## Contour smoothed with smoothr                                               
volcano_contour_smooth <-                                                      
  smoothr::smooth(volcano_contour, method = "ksmooth", smoothness = 5)                                                                                         

## Increase resolution                                                         
v_res_hi <- raster::disaggregate(v, 3, method = "bilinear")                                                                                                    
volcano_contour_res_hi <- toIsobandSf(v_res_hi)                                

## Blur raster                                                                 
v_res_hi_blur <-                                                               
  focal(v_res_hi, w=matrix(1, 5, 5), mean)                                     
volcano_contour_res_hi_blur <- toIsobandSf(v_res_hi_blur)                                                                                                      

## Blur raster gauss                                                           
gf <- focalWeight(v_res_hi, c(1, 1), "Gauss")                                  
v_res_hi_blur_gauss <-                                                         
  focal(v_res_hi, w = gf)                                                      
volcano_contour_res_hi_blur_gauss <- toIsobandSf(v_res_hi_blur_gauss)                                                                                          

## Plot data                                                                   
p1 <- plotContour(volcano_contour)                                             
p2 <- plotContour(volcano_contour_smooth)                                      
p3 <- plotContour(volcano_contour_res_hi)                                      
p4 <- plotContour(volcano_contour_res_hi_blur)                                 
p5 <- plotContour(volcano_contour_res_hi_blur_gauss)                           

p <-                                                                           
  cowplot::plot_grid(p1, p2, p3, p4, p5                                        
                   , labels = c('Raw', 'smoothr::smooth', 'Raw Res. High'                                                                                      
                              , 'Raw Res. High, Raster blur'                                                                                                   
                              , 'Raw Res. High, Raster blur Gauss'                                                                                             
                                )                                              
                   , hjust = 0)                                                

ggsave('volcano.png', p, width = 210, height = 148, units = 'mm')