raster resolution not displaying properly

Question

I'm seeing some strange behavior when I plot a certain raster. This is a shaded relief raster I acquired from naturalearthdata.com.You can download it here.

I find that depending on what spatial scale I plot the raster at, it plots a different resolution, despite the resolution of the raster not being changed.

library(raster)
relief <- raster('GRAY_50M_SR_W.tif')

# let's use Mexico as an example:
library(maptools)
data(wrld_simpl)
mx <- wrld_simpl[which(wrld_simpl@data$NAME == 'Mexico'),]

# Here I create a cropped version of the raster
reliefMX <- crop(relief, mx)

To illustrate the problem, I plot Mexico to get the map extent, then I plot the full extent of the raster, and then the cropped raster on top.

You can see that the rasters show very different resolutions, but they really have identical resolutions.

plot(mx)
plot(relief, col = grey(0:100/100), legend = FALSE, axes = F, box = F, add=T)
plot(reliefMX, col = grey(0:100/100), legend = FALSE, axes = F, box = F, add=T)

> res(relief)
[1] 0.03333333 0.03333333
> res(reliefMX)
[1] 0.03333333 0.03333333

Any ideas? How can I get these rasters to display properly?

> sessionInfo()
R version 3.4.0 (2017-04-21)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.4

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] maptools_0.9-2 raster_2.5-8   sp_1.2-4      

loaded via a namespace (and not attached):
[1] compiler_3.4.0  rgdal_1.2-7     tools_3.4.0     foreign_0.8-68  Rcpp_0.12.10   
[6] grid_3.4.0      lattice_0.20-35

lbusett lbusett · Accepted Answer · 2017-05-16T08:20:53

This depends on the maxpixels argument in the call to raster::plot:

maxpixels integer > 0. Maximum number of cells to use for the plot. If maxpixels < ncell(x), sampleRegular is used before plotting. If gridded=TRUE maxpixels may be ignored to get a larger sample

When you plot the "full" map, the image is automatically downsampled to save you memory and reduce the rendering time. You can change the value of `maxpixels' to get the desired level of "detail". see for example:

plot(relief, col = grey(0:100/100), legend = FALSE, axes = F, box = F, add=F)

plot(relief, col = grey(0:100/100), legend = FALSE, axes = F, box = F, add=F, maxpixels = 5000000)

Although at this "zoom level" is not evident, the second image is more detailed. You can appreciate this by "zooming in" on your cropped area:

plot(mx)
plot(relief, col = grey(0:100/100), legend = FALSE, axes = F, box = F, add = T, maxpixels = 5000000)

Still not "good" as the a-priori cropped one (since I'm still not using "all" pixels), but already better.

In practice, it's all a compromise between rendering time/memory and quality of the output. Obviously, if you just need to plot a part of the area, it's much more efficient to crop the image beforehand.

HTH.

raster resolution not displaying properly

1 Answers