Display two-dimensional color gradient for two different quantities

Question

I have two overlapping raster layers. One contains quantity A from 0 to 100, and the other contains quantity B from 0 to 100. I'd like to create a map displaying both sets of data at once based on a color box with two different axes.

So, for an area where:

• A is 100 and B is 100, the map would be red
• A is 100 and B is 0, it would be yellow
• A is 0 and B is 100, it would be green
• A is 0 and B is 0, it would be blue

For intermediate values the hue would change accordingly, so A=100, B=50 would be orange; A=50, B=0 would be blueish-green.

I have seen maps display this type of data before, and I have a couple of things I'd like to display this way. What is the easiest way to do this? I can use ArcGIS, QGIS and R.

Answer 1

What you're trying to create is known as a bivariate map. There's a couple of ways to do this, but since you've already got raster data that leans toward certain methods. The big thing is going to be your color selections, and I'm not sure exactly how to get the blending you desire.

Typically this is done with two colors, one for each variable. So A 0 = white, A 100 = blue, B 0 = white, and B 100 = red. You would then merge the two rasters and get shades of purple where both are high, white where both are low, and a color leaning toward red or blue where one is greater than the other.

This can be adapted to two color axis. A 0 is yellow, A 100 is red and B 0 is green with B 100 blue. But then you can start to get some muddy values in the middle, depending on your colors. In this case, A and B 100 is purple, while A and B 0 is a yellow green. A 0 and B 100 is going to be gray while A 100 and B 0 is olive.

With just the three primaries, the blends are fairly clean. That fourth complicates things though. According to your question, both A and B 0 would be blue, A 100 would be yellow, B 100 is green, and both A and B 100 is red. This doesn't really work since yellow plus green isn't red. And yellow plus the midpoint of green and blue isn't orange.

I think you'll have better results if you choose just two colors. You might even get away with three (both are the same primary color low or high, but then each has a different primary for the other end of the range). But the basic principle is the same. You'll symbolize each of your rasters with the desired color ramp, export/save that to a colorized raster, then merge/blend the two colorized rasters together to get your final bivariate raster. There are different approaches to merge/blend as well, and some software gives more options than others. For example it could be a transparency blend or an additive color blend.

There's an example specific method at this question using gdal.

Answer 2

If you would like to do this in R, you might want to have a look at the blog post by Oscar Perpiñán Lamigueiro. Its not exactly what you are looking for, but there might be a way to modify it to work for your purposes.

In my last post I described how to produce a multivariate choropleth map with R. Now I will show how to create a map from raster files. One of them is a factor which will group the values of the other one. Thus, once again, I will superpose several groups in the same map.

Addition:

It may be a bit late, but I thought I would add this as it seems like a relatively simple solution to your problem.

In R, when you use plot() to display rasters, there is an alpha argument that sets the transparency. If you plotted your two variables as different color rasters with color scales ranging from light to dark and each alpha argument set to something like 0.5, the colors would combine where they overlap.

I have done this with spatialPolygons using spplot as well. I imagine it would work the same for rasters. Here is an example with spatialPolygons:

library(sp)
library(rgeos)
library(rworldmap)

box <- readWKT("POLYGON((-180 90, 180 90, 180 -90, -180 -90, -180 90))")
proj4string(box) <- CRS("+proj=cea +datum=WGS84")
set.seed(1)
pts <- spsample(box, n=2000, type="random")
pols <- gBuffer(pts, byid=TRUE, width=50) # create circle polys around each point
merge = sample(1:40, 100, replace = T) # create vector of rand #s between 0-100 to merge pols on

Sp.df <- gUnionCascaded(pols, id = merge) # combine polygons with the same 'merge' value
# create SPDF using polygons and randomly assigning 1 or 2 to each in the @data df
Sp.df <- SpatialPolygonsDataFrame(Sp.df, data.frame(z = factor(sample(1:2, length(Sp.df), replace = TRUE)),
                                                    row.names= unique(merge)))
Sp.df <- crop(Sp.df, box)
colors <- c(rgb(r=0, g=0, blue=220, alpha=50, max=255), rgb(r=220, g=0, b=0, alpha=50, max=255))

land <- getMap()

overlay.map <- spplot(Sp.df, zcol = "z", col.regions = colors, col = NA, alpha = 0.5, breaks=c(0,1)) +
  layer(sp.polygons(land, fill = "transparent",  col = "grey50"))

Obviously, the legend is not very helpful for this map. To create a legend, you then need to plot the gradients in a square with one increasing left to right and the other increasing bottom to top.

Some code to do this in base graphics is here: https://stackoverflow.com/a/11103414/3897439.

I also did it for two colors using ggplot as I needed grid graphics. Here is a simple example:

Variable_A <- 100 # max of variable
Variable_B <- 100

x <- melt(outer(1:Variable_A, 1:Variable_B)) # set up the data frame to plot from

p <- ggplot(x) + theme_classic() + scale_alpha(range=c(0,0.5), guide="none") +
  geom_tile(aes(x=Var1, y=Var2, fill="Variable_A", col.regions="red", alpha=Var1)) +
  geom_tile(aes(x=Var1, y=Var2, fill="Variable_B", col.regions="blue", alpha=Var2)) +
  scale_x_continuous(limits = c(0, Variable_A), expand = c(0, 0)) +
  scale_y_continuous(limits = c(0, Variable_B), expand = c(0, 0)) +
  xlab("Variable_A") + ylab("Variable_B") +
  guides(fill=FALSE)
p

It's not perfect, but may work for your purposes.

Another source that might have some interesting information relevant to your task is: http://andrewpwheeler.wordpress.com/2012/08/24/making-value-by-alpha-maps-with-arcmap/.