CODEWITHSUNDEEP
rrasterr-raster
Pascal
Pascal

Reputation: 1672

raster resolution not displaying properly

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

enter image description here

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

Upvotes: 2

Views: 486

Answers (1)

lbusett
lbusett

Reputation: 5932

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)

enter image description here

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

enter image description here

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)

enter image description here

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.

Upvotes: 1

Related Questions