[GIS] Smoothing raster map using R

mapsrrastersmoothing

I am trying to create a map from raster data. The file came from a crop model, with resolution of 0.5 degree. Even when I disaggregate it (i.e. increase spatial resolution), the map looks really pixelated. I am trying to make it look better.

My current code produces this image:

generated with my code

where I would like to "smooth" the data, by supressing the pixelated look. Some other visualization programs do this automatically, so I guess it should not be hard to reproduce using R.

For example, this is the same file plotted using Panoply:

generated with Panoply

It doesn't look absolutely smooth, but at least it doesn't have the pixelated look neither. How to achieve a similar result in R?

This is the code to reproduce my problem:

library(RCurl) 
library(rasterVis)

# Go to temp dir and download file - approx. 1.7M
old <- setwd(tempdir())

# download raster and shapefile
download.file('https://dl.dropboxusercontent.com/u/27700634/yield.nc', 'yield.nc', method='curl')
download.file('https://dl.dropboxusercontent.com/u/27700634/southern.zip', 'southern.zip', method='curl')
unzip('southern.zip', exdir='.')

# load southern Brazil shapefile
mapaSHP <- shapefile('southern.shp')

# load brick
b <- brick('yield.nc', level=16)

# create color scheme
mycols <- rasterTheme(region=colorRampPalette(brewer.pal(9,'Greens'))(100))

# use second brick layer to plot map
levelplot(b[[2]], margin = FALSE, main = "Rice yield in tons/ha", par.settings = mycols) +
  layer(sp.lines(mapaSHP, lwd=0.8, col='darkgray'))

# return to your old dir
setwd(old)

Best Answer

Here are some ideas.

With base plot you can do

 plot(x, interpolate=TRUE)

You can also resample your data

 y <- disaggregate(x, 5, method='bilinear')

Or indeed smooth it using a focal operation

 y <- focal(x, w=matrix(1, 5, 5), mean)

Or a combination

 y <- disaggregate(x, 5)
 y <- focal(y, w=matrix(1, 5, 5), mean)

The question whether doing this is a good idea or not is another matter, that I'll leave to you to decide

Related Solutions

[GIS] Performing Raster Noise Reduction and Edge Smoothing

The morphological operations Expand and Shrink were created for this kind of processing. Use ArcGIS (or GRASS or Mathematica) because R's "raster" library is too slow.

Often it helps to experiment a little with the parameters: you have to decide how much expanding and shrinking is needed to clean an image; and usually you want to do as little as possible, because each operation tends to smooth out some of the sharp details. Here is a sequence that works well to eliminate much of the apparent "noise" while maintaining most of the detail in the "entities". "Expand" and "shrink" are both with reference to the white cells, so that expanding causes them to grow outwards and shrinking causes the black cells to encroach into white regions.

Table of images

The "difference" column uses color to highlight differences between the start and end image at each step: blue for black that turned to white, and orange for white that turned to black.

If the larger remaining pieces need to be removed, that might best be done with RegionGroup to identify them, after which they can be obliterated through reclassification. This was an option at the outset, but a little initial cleaning with Expand and Shrink reduces the work and provides the desired smoothing.

Incidentally, I chose to make the eight images in this illustration with Mathematica commands because they are so simple, easy, and fast to execute:

i = Import["http://i.stack.imgur.com/umDg7.png"];
l = Dilation[k = Erosion[j = Dilation[i, 2], 3], 1];      (* This does all the work *)
delta = ColorCombine /@ {{i, j}, {j, k}, {k, l}, {i, l}}; (* Compares images *)

The workflow in ArcGIS is the same but the syntax will be lengthier. If you really want to use R, load the "raster" library and exploit focalFilter to create functions to do the expanding and shrinking. Then wait about a minute each to execute the operations... .

[GIS] Raster time series smoothing package using R

Rather than looking for a specific package for raster time series you could look for functions for smoothing, and then use these with the calc function in the raster package. Here is an example for Savitzky-Golay:

https://stackoverflow.com/questions/37843942/smoothen-rasterstack-using-the-savitzky-golay-sgolayfilt-signal-in-r/37846229#37846229

Best Answer

Related Solutions

[GIS] Performing Raster Noise Reduction and Edge Smoothing

[GIS] Raster time series smoothing package using R

Related Question