Agustin.Lobo@ictja.csic.es
20140508
Last update: 20140514
file: /media/alobo/Iomega_HDD/RrasterGuide/rasterVisTutorial/rasterOnGM.R
GM tiles can be downloaded as R objects using either ggmap:get_map() or dismo:gmap().
We prefer dismo:gmap() because, unlike ggmap:get_map(), it returns a raster layer from package raster including complete CRS information
require(utils)
require(colorRamps)
require(rgdal)
require(raster)
require(dismo)
require(rasterVis)
As an example, we will be using a field of the probability that a given pixel has got > 1kg/m2 of ash at a given eventual eruption of the Cerro Blanco volcano (Argentina) during the period of study
The probability field is a raster layer in GeoTif format
download.file("https://dl.dropboxusercontent.com/u/3180464/rprob520.tif", "rprob520.tif",
method = "curl")
rprob <- raster("rprob520.tif")
CAUTION: by some reason, using a remote file, i.e.:
# delme <-
# raster('https://dl.dropboxusercontent.com/u/3180464/rprob520.tif')
CRASHES RStudio
rprob
## class : RasterLayer
## dimensions : 161, 200, 32200 (nrow, ncol, ncell)
## resolution : 10000, 10000 (x, y)
## extent : 418200, 2418200, 6231355, 7841355 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : /media/alobo/Iomega_HDD/RrasterGuide/rasterVisTutorial/rprob520.tif
## names : rprob520
## values : 0.007692, 1 (min, max)
extent(rprob)
## class : Extent
## xmin : 418200
## xmax : 2418200
## ymin : 6231355
## ymax : 7841355
plot(rprob)
projection(rprob)
## [1] "+proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
gmap() requires the extent to be in either longitude/latitude (epsg:4326) or Google PseudoMercator (epsg:3857),
but can also accept an spatial object in any CRS provided it is correctely defined: gmap() will perform the reprojection if needed.
As extent objects are not geospatial objects (they lack the CRS information), in case you provide the extent then you must take care
of the eventual reprojection on your own first.
The simplest way is to ask gmap() to download the GM layer centered on our raster object at a given zoom level:
migmap <- gmap(x = rprob, type = "hybrid", zoom = 5)
plot(migmap)
But note the following would not be what we expect:
delme <- gmap(x = extent(rprob), type = "hybrid", zoom = 5)
plot(delme)
because in case we provide the extent and not the spatial object we must deal with the reprojection on our own:
delme <- gmap(x = projectExtent(rprob, crs = CRS("+init=epsg:3857")), type = "hybrid",
zoom = 5)
plot(delme)
rm(delme)
We deal with projections and reprojections a bit more in depth elsewhere.
In order to overlay rprob on top of migmap, we must make both layers have the same CRS. As this is a large area, we select Google PseudoMercator (epsg:3857) and reproject rprob:
rprobGM <- projectRaster(from = rprob, crs = CRS("+init=epsg:3857"))
The reprojection involves some interpolation, thus we make sure the |0,1| interval is fulfilled:
rprobGM[rprobGM < 0] <- 0
rprobGM[rprobGM > 1] <- 1
plot(migmap)
plot(rprobGM, add = T, legend = F, col = rev(rainbow(10, alpha = 0.35)))
tune the palette:
get pal() from HCL-Based Color Palettes in R
# pal(rev(rainbow(10, alpha=1)))
micolor <- rev(rainbow(12, alpha = 0.35))
transp <- rainbow(12, alpha = 0)
micolor[1:3] <- transp[1]
plot(migmap)
plot(rprobGM, add = T, legend = F, col = micolor)
The layer returned by gmap() is pseudocolor (indexed 256 colors) raster that must be represented with its own colortable, which is included as an slot within slot legend in the raster layer object:
length(migmap@legend@colortable)
## [1] 256
migmap@legend@colortable[1:3]
## [1] "#00081A" "#111A11" "#1A1F16"
We must specify the breaks between colors with at the 256 levels
colors <- migmap@legend@colortable
migmaplv <- levelplot(migmap, maxpixels = ncell(migmap), col.regions = colors,
at = 0:255, panel = panel.levelplot.raster, interpolate = TRUE, colorkey = FALSE,
margin = FALSE)
print(migmaplv)
levelplot of the rprob layer must include an opacity factor (alpha.regions) < 1 to be partially transparent:
rproblv <- levelplot(rprobGM, margin = FALSE, contour = TRUE, par.settings = rasterTheme(region = matlab.like(n = 10)),
alpha.regions = 0.35, at = (0:10)/10, main = "p(dep > 1kg/m2 per eruption event)")
print(rproblv)
and finally we overlay the two graphics:
print(migmaplv + rproblv)
