Handling outliers in a choropleth map plotted with ggplot2

Question

I want to do some plotting of map data in ggplot where the polygons are filled according to values of a continuous variable. One problem I have is that there are just a handful of extreme outliers (~10 out of 3000) in the data that prevent most of the other observations from being displayed with enough meaningful color variation.

Basically I'd like to be able to specify a color scheme for the bulk of my data, then color all positive outliers bright blue and all negative outliers bright red.

Example Code - For the sake of example, let's say that most of the data falls between -0.01 and 0.01, but there are a few values that fall way below/above these levels...

library(ggplot2)

DAT <- structure(list(long = c(848025.138769486, 827715.400155344, 819783.06692123, 
    811273.46685221, 807112.24435793, 825218.389316822, 824653.263676394, 
    851268.022958262, 854876.20206792, 848025.138769486, 824710.033713034, 
    808060.98832434, 791388.306466439, 791011.02400712, 790256.857765555, 
    790077.133379578, 792915.8857683, 803799.11479841, 812042.67207171, 
    811086.452951819, 818801.811870496, 824623.483334553, 829923.458416495, 
    844732.373710393, 843509.962977133, 840793.468791028, 844661.128104187, 
    841054.158477243, 824318.536510595, 824710.033713034, 957744.153625377, 
    964378.49410997, 963723.372709918, 963536.949007208, 950656.36657876, 
    923479.701493794, 918733.432704199, 950416.579446804, 951806.914646236, 
    952020.610154412, 957744.153625377, 838474.166485858, 839637.34900613, 
    846161.553188089, 854105.657535852, 850749.837928768, 851268.022958262, 
    824653.263676394, 825218.389316822, 807112.24435793), lat = c(-755063.209671518, 
    -751075.227428769, -758486.024742793, -742209.818502709, -725584.096412241, 
    -723927.958082718, -719335.285579824, -716264.904081879, -747370.371610989, 
    -755063.209671518, -607299.311749675, -601135.889250199, -605563.568360578, 
    -602482.006717006, -602560.163831817, -601169.088819494, -595490.118653201, 
    -589823.673980013, -577642.515107293, -574207.831414412, -570499.93978401, 
    -555676.769699399, -554574.027849141, -569849.400696679, -581914.546950335, 
    -587110.213020993, -586640.550459734, -595326.416960721, -603469.770252399, 
    -607299.311749675, -686468.891810573, -696059.793425604, -696994.182842627, 
    -698894.537935268, -702169.738050773, -705710.496557967, -666147.504830981, 
    -662203.856856129, -666852.009638632, -667571.668902733, -686468.891810573, 
    -684524.444323129, -694047.01193517, -705551.617269712, -711060.296586236, 
    -711747.897688833, -716264.904081879, -719335.285579824, -723927.958082718, 
    -725584.096412241), group = c("0.1", "0.1", "0.1", "0.1", "0.1", 
    "0.1", "0.1", "0.1", "0.1", "0.1", "1.1", "1.1", "1.1", "1.1", 
    "1.1", "1.1", "1.1", "1.1", "1.1", "1.1", "1.1", "1.1", "1.1", 
    "1.1", "1.1", "1.1", "1.1", "1.1", "1.1", "1.1", "2.1", "2.1", 
    "2.1", "2.1", "2.1", "2.1", "2.1", "2.1", "2.1", "2.1", "2.1", 
    "3.1", "3.1", "3.1", "3.1", "3.1", "3.1", "3.1", "3.1", "3.1"
    ), VALUE = c(0.0197178619295337, 0.0197178619295337, 0.0197178619295337, 
    0.0197178619295337, 0.0197178619295337, 0.0197178619295337, 0.0197178619295337, 
    0.0197178619295337, 0.0197178619295337, 0.0197178619295337, 0.0144402919365254, 
    0.0144402919365254, 0.0144402919365254, 0.0144402919365254, 0.0144402919365254, 
    0.0144402919365254, 0.0144402919365254, 0.0144402919365254, 0.0144402919365254, 
    0.0144402919365254, 0.0144402919365254, 0.0144402919365254, 0.0144402919365254, 
    0.0144402919365254, 0.0144402919365254, 0.0144402919365254, 0.0144402919365254, 
    0.0144402919365254, 0.0144402919365254, 0.0144402919365254, -0.00812118892018265, 
    -0.00812118892018265, -0.00812118892018265, -0.00812118892018265, 
    -0.00812118892018265, -0.00812118892018265, -0.00812118892018265, 
    -0.00812118892018265, -0.00812118892018265, -0.00812118892018265, 
    -0.00812118892018265, 0.00751936235807205, 0.00751936235807205, 
    0.00751936235807205, 0.00751936235807205, 0.00751936235807205, 
    0.00751936235807205, 0.00751936235807205, 0.00751936235807205, 
    0.00751936235807205)), .Names = c("long", "lat", "group", "VALUE"
    ), row.names = c(6L, 5L, 4L, 10L, 7L, 8L, 9L, 2L, 3L, 1L, 151L, 
    150L, 156L, 162L, 168L, 163L, 159L, 154L, 160L, 149L, 161L, 158L, 
    164L, 152L, 153L, 155L, 165L, 167L, 166L, 157L, 226L, 223L, 233L, 
    232L, 225L, 227L, 224L, 230L, 228L, 229L, 231L, 258L, 253L, 257L, 
    262L, 254L, 256L, 261L, 260L, 269L), class = "data.frame")

#Construct plot, placing limits on color scale. How to make outliers bright red or blue?
gmp <- ggplot(data=DAT, aes(x=long,y=lat,group=group)) +
scale_fill_gradient2(low=rgb(0.8,0.2,0.2),high=rgb(0.2,0.3,0.8),mid=rgb(0.9,0.9,0.9),limits=c(-0.01,0.01)) + 
geom_polygon(aes(fill=VALUE,group=group),colour=NA, size = .3) + 
coord_fixed()
plot(gmp)

Answer 1

Here's one option:

Create subsets for each criteria (minimum extremes, maximum extremes, and mid ranges). I've added a minimum extreme to your data for illustration.

DAT$VALUE[DAT$group == 0.1] <- -0.019
DAT.mid <- DAT[abs(DAT$VALUE) < 0.01,]
DAT.max <- DAT[DAT$VALUE >= 0.01,]
DAT.min <- DAT[DAT$VALUE <= -0.01,]

Then build a base plot with just the mid data:

base <- ggplot(DAT.mid, aes(x=long,y=lat,group=group)) +
  geom_polygon(aes(fill=VALUE,group=group),colour=NA, size = .3) + 
  scale_fill_gradient2(low=rgb(0.8,0.2,0.2),high=rgb(0.2,0.3,0.8),
  mid=rgb(0.9,0.9,0.9),limits=c(-0.01,0.01)) + 
  coord_fixed()
base 

Then you can add layers with your extremes with whatever colours you want:

ext <- base +
  geom_polygon(data = DAT.max, fill = "blue") +
  geom_polygon(data = DAT.min, fill = "red") 
ext

The tricky bit is to get a legend representing your extremes. Following on from here, you can add invisible geoms mapped to non-existent aes, and then ab/use the legend to get it how you want. Here, I've used to geom_points. For the first, I've mapped size to ">= 0.01", which will give us a legend with that text, and similar for the second, except I've used alpha to get another legend entry. You can then use guides() to override the legend attributes, and get a square with the right fill. It's not perfect, but works in most circumstances.

gmp <- ext +
  geom_point(data = DAT.max, aes(size = ">= 0.01", shape = NA), colour = "blue") +
  geom_point(data = DAT.min, aes(alpha = "<= -0.01", shape = NA), colour = "red") +
  guides(size = guide_legend(title = "", override.aes = list(shape = 15, size = 10)),
         alpha = guide_legend(title = "", override.aes = list(shape = 15, size = 10, alpha = 1)))
gmp