1

I know this variants of question already exist, but I'm still stuck: How do I plot specific categorical levels of a raster in R, using sppplot or ggplot?

Right now I have a raster layer holding soil zinc values, called zn. Here is the information:

    class       : RasterLayer 
    dimensions  : 1308, 3188, 4169904  (nrow, ncol, ncell)
    resolution  : 250, 250  (x, y)
    extent      : -178002.4, 618997.6, 2914810, 3241810  (xmin, xmax, ymin, ymax)
    coord. ref. : +proj=utm +zone=45 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 

I can plot that variable, along with administrative outlines held in the npadmin1 shapefile, using this code:

   spplot(zn, scales = list(draw = FALSE), 
   col.regions = terrain.colors(100)) + 
   layer(sp.polygons(npadmin1, lwd = 1))

But I also want to plot 4 levels of the raster value, in 4 colors: 0-.5, .5-1, 1-1.5, >1.5. And then I want the legend to say "low" "medium", "high", "very high". No color bar. I know there are similar questions out there, but (a) I can't seem to get gplot to work at all, and (b) I can't figure out how to do this in spplot. To be honest, I don't really know enough about spplot, ggplot and/or levelsplot to know when each is more appropriate to use.

I'm already aware of the 2 links below... still can't figure it out.

Legend of a raster map with categorical data

Plot continuous raster data in binned classes with ggplot2 in R

Leah Bevis
  • 317
  • 1
  • 11
  • But, why you can achieve this? Just reclass, add Raster Attribute Table and plot with `levelplot()` – aldo_tapia Jul 24 '18 at 21:01
  • I actually did reclassify, but (a) it did something odd to the raster, and (b) I still couldn't figure out how to plot the reclassified raster. I tried, couldn't get the legends to work. Like I said, I'm not great with any of these 3 (levelplot,spplot,ggplot) so if you could give an example of the code, that would be great. – Leah Bevis Jul 25 '18 at 04:09

2 Answers2

2

As I said in commentaries: reclassify, add raster attribute table and plot with levelplot():

library(raster)

# Reproducible example
r <- raster()
r[] <- runif(ncell(r), min = 0, max = 2)

# Reclassify
r <- reclassify(r, c(0, 0.5, 1,
                     0.5, 1, 2,
                     1, 1.5, 3,
                     1.5,Inf,4))

# View
plot(r)

enter image description here

# Values as factor
r <- as.factor(r)
# Extract attribute table
rat <- levels(r)[[1]]
# Set custom breaks
rat[["zn"]] <- c("low", "medium", "high", "very high")
# Add back RAT
levels(r) <- rat

# Plot
rasterVis::levelplot(r)

enter image description here

aldo_tapia
  • 1,153
  • 16
  • 27
  • Thanks for the example! However, this doesn't work for me. First of all, you're not specifying the break values, so I'm guessing they are set automatically at something like 25th, 50th, or 75th quantiles? I want, instead, to set mine to specific values (.5, 1, 1.5). Also, when I do run the code anyway, allowing automatic breaks, I get the following error: "Error in .local(x, y, ...) : duplicate "by" values not allowed" – Leah Bevis Jul 25 '18 at 14:15
  • There are no automatic breaks. Breaks are defined by `reclassify()` function. If you want to ser quantiles breaks, you need to compute those quantiles and apply a reclassification. – aldo_tapia Jul 25 '18 at 14:45
  • Oh sorry! I had missed the code at the top, my apologies. Now it works perfectly, thanks you so much! One last question, if you'll allow me. I'm successfully mapping the raster on a defined color scale, but when I try to overlay a polygon called npadmin1 (code: "+ layer(sp.polygons(npadmin1, lwd = 1)") it give me this error: "Only one argument supplied to binary operator + which requires two." I see people using this code online w/ levelplot, and I've successfully used it with spplot. Any idea why it's not working here? I can't interpret the error. – Leah Bevis Jul 25 '18 at 15:17
  • Check https://oscarperpinan.github.io/rastervis/FAQ.html, I have no problems with a similar code using `layer(sp.polygons())` – aldo_tapia Jul 25 '18 at 18:25
1

Here is a ggplot2 solution:

library(raster)
library(ggplot2)
library(dplyr)

# Reproducible example
r <- raster()
r[] <- runif(ncell(r), min = 0, max = 2)

# Reclassify
r <- reclassify(r, c(0, 0.5, 1,
                     0.5, 1, 2,
                     1, 1.5, 3,
                     1.5,Inf,4))

r_df <- as.data.frame(r, xy = TRUE) %>%
  mutate(layer = factor(layer))

# View
ggplot() + 
  geom_raster(data = r_df, aes(x = x, y = y, fill = layer)) +
  scale_fill_manual(values = viridis::viridis(4), breaks = 1:4,
                    labels = c("low", "medium", "high", "very high"))

enter image description here

jsta
  • 3,216
  • 25
  • 35