[R-sig-Geo] How to “smooth” a raster map

[Thiago V. dos Santos]

Thank you all for the comments; I really appreciate them.

Julian, interpolation was what I was looking for. Instead of 'resample', I ended up using 'disaggregate' with the argument 'method', which I was not using on my code and made a complete difference.

Roger (and Julian too), I guess this is not the case of hiding the actual resolution or cheating. Had I thought about the repercussion of this ethical issue before, I would have made the final purpose clear from the beginning. For that, I am sorry.

The crop model I am using employs mainly climate data as input to solve the equations. There is no information on the spatial distribution of crops. Therefore, the output data (maps) assume that crops are present everywhere on the simulation domain.

A next step on my code (which I omitted in my question for simplification purposes and was not a good idea) is to mask the model output using another map [1], that shows a finer distribution of the studied crop.

Therefore, instead of for example having a southern Brazil completely covered with rice (over water bodies too? - remember this is derived from climate date) like this: http://i.imgur.com/3pPpRXp.png

we can have a more realistic spatial distribution like this, taking into account a high-resolution, satellite-derived map of rice distribution and showing only the areas where rice is actually grown: http://i.imgur.com/HiKqCQW.png

I should reinforce that, for now, no analysis is being conducted on this map: it is for visualization only. Also, all the information on the input dataset (resolution included) will be provided on the supportive text. However, for future works, I still think that the second map would be more suitable for any kind of quantitative analysis. I wonder if some lies are for good?

Thank you guys for bringing up this issue. I would like to hear feedback from other people.

[1] from Monfreda et al. (2008), "Farming the planet: 2. Geographic distribution of crop areas, yields, physiological types, and net primary production in the year 2000", Global Biogeochemical Cycles, Vol.22
 
Greetings,
--
Thiago V. dos Santos
PhD student
Land and Atmospheric Science
University of Minnesota
http://www.laas.umn.edu/CurrentStudents/MeettheStudents/ThiagodosSantos/index.htm
Phone: (612) 323 9898



On Monday, June 29, 2015 12:49 AM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
Julian,

A measured response, with an answer but crucially your comments further 
down. I really hope Thiago reads them, as they are of central importance.

Even for display, nobody should ever hide the actual resolution of the 
information being displayed (even if others do this in error). Of course 
the authors of spatial packages in R will never provide default displays 
that intentionally deceive the viewer.

If the half degree output for the model was what it was designed to 
provide, nobody can know the fitted values at higher resolution without 
re-running the model itself at higher resolution. The input data to the 
model may not be available at this resolution.

More important, the model output almost certainly is accompanied by 
measures of the prediction uncertainly, so that each half-degree cell 
value is actually a summary of the predicted distribution.

Trying to smooth only a central tendency measure of these distributions 
deterministically is creating complete chaos - you will not know what the 
resampled cell distributions are. That is why creating a "nice map" is a 
really bad idea, see also:

http://www.markmonmonier.com/how_to_lie_with_maps_14880.htm

The pixelation is actually your friend, because it is communicating the 
support of the model fitted values visually.

Again, thanks to Julian for a measured and very rapid response.

Roger


On Mon, 29 Jun 2015, Julian Burgos wrote:

> Hi Thiago,
> If the output of your model has a resolution of 0.5 degrees, you will have
> to do some kind of interpolation to get the "smooth" look that you are
> looking for.  If you are only doing this for visualization purposes, you
> can use the resample function and do a simple bilinear interpolation.  The
> function goes something like this:
>
> new.raster <- b[[2]] # Create a new raster (with same extent, etc. as your
> original raster)
> res(new.raster) <- 0.25 # Change the resolution.. select whatever value
> you want... small values require more time
> resample(b[[2]], new.raster, method="bilinear")
> levelplot(new.raster)
>
> Now, remember that when you do this you are in a way cheating.  You are
> showing a model output at much higher resolution that the output really
> is.  But again, if it is only to have a pretty picture then it is fine. 
> On the other hand, if you are going to use the new.raster for other
> analysis or as input for other models, then things get complicated.
>
> All the best,
>
> Julian
> -- 
> Julian Mariano Burgos, PhD
> Hafrannsóknastofnun/Marine Research Institute
> Skúlagata 4, 121 Reykjavík, Iceland
> Sími/Telephone : +354-5752037
> Bréfsími/Telefax:  +354-5752001
> Netfang/Email: julian at hafro.is
>
>> Dear all,
>>
>> 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: http://i.stack.imgur.com/WssPy.png
>>
>> 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:
>> http://i.stack.imgur.com/jXYI7.png
>>
>> 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)
>>
>> --------------------------------------------------------------------------
>> Thanks in advance for any input,
>> --
>> Thiago V. dos Santos
>> PhD student
>> Land and Atmospheric Science
>> University of Minnesota
>> http://www.laas.umn.edu/CurrentStudents/MeettheStudents/ThiagodosSantos/index.htm
>> Phone: (612) 323 9898
>>

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
e-mail: Roger.Bivand at nhh.no