Next Previous

Introduction

I posted an oce issue recently [1] that indicated my displeasure with the present defaults for the hues of Landsat plots that use band="terralook". Others have reported similar problems, and this was enough to motivate some experimentation.

I choose a region with which I am very familiar, namely the city in which I live, Halifax Nova Scotia. I happen to have two Landsat scenes on my computer, one from summer and one from winter, so I used them for my tests.

Readers who are unfamiliar with Halifax might want to consult a map to understand place names in the following discussion of my reference points in the Landsat images.

First, note that the chosen region contains a fair bit of ocean, but I lack ground-truth on its colour, so I have no way to know if the hue on a given day should be blueish, greenish, or grayish.

The scene also does not contain much in the way of red reference material, nor yellow. I suppose other Landsat images show fields or deserts that would be good for such colours, and I will on other people working on such domains to let me know whether plot.landsat will need further adjustment for those hues.

In my reference images, I have three main elements: green areas, snow-covered areas, which should have a light and neutral hue, and asphalt-covered areas that should be dark and neutral.

Year-round, the region of Point Pleasant Park should be a green hue, because the trees there are mostly coniferous. In the rest of the city, trees are mostly deciduous, and so comparing summer and winter images will be helpful there. The southern dockyard region and the train-tracks leading to it are either covered by asphalt or gravel, and so should be neutral in tone.

Procedure

The first step in the code is to create local files. This saves time, because reading raw Landsat files is time consuming. (Note that the files are checked into the git repository holding this blog, so that savvy readers can reproduce the other steps of this posting without downloading 2Gb worth of Landsat files.)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
library(oce)
## 1. Create local-view datasets
Hc <- c(-63.57, 44.63)                 # south end of Halifax, NS
Hd <- 0.025 * c(1/cos(pi*Hc[2]/180), 1) # approx 5 km box
band <- c("red", "green", "nir") # lets us plot 'terralook'
## Wintertime
if (0==length(list.files(pattern="^Hw.rda$"))) {
    w <- read.landsat("/data/archive/landsat/LC80080292014065LGN00", band=band)
    Hw <- landsatTrim(w, ll=Hc-Hd, ur=Hc+Hd)
    save(Hw, file="Hw.rda")
    rm(w)
} else {
    load("Hw.rda")
}

## Summertime
if (0==length(list.files(pattern="^Hs.rda$"))) {
    s <- read.landsat("/data/archive/landsat/LC80070292014170LGN00", band=band)
    Hs <- landsatTrim(s, ll=Hc-Hd, ur=Hc+Hd)
    save(Hs, file="Hs.rda")
    rm(s)
} else {
    load("Hs.rda")
}

At this stage, we have landsat-class objects named Hw and Hs to work with. It is convenient to set up a demo function that will plot a four-panel plot.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
demo <- function(red.f, green.f, blue.f, name="")
{
    if (!missing(name)) png(name, unit="in", width=6, height=3.1, res=150)
    par(mfrow=c(1,2))
    mar <- c(0.25, 0.25, 1, 0.25)
    mar <- c(2, 2, 1.5, 0.5)
    axes <- TRUE
    cex <- 0.8
    plot(Hw, band='terralook', red.f=red.f, green.f=green.f, blue.f=blue.f,
         mar=mar, axes=axes, cex=cex)
    mtext(sprintf("red.f=%g, green.f=%g, blue.f=%g", red.f, green.f, blue.f),
          side=3, line=0, adj=0, cex=cex)
    mtext(format(Hw[['time']], "%B %Y"), side=3, line=0, adj=1, cex=cex)
    plot(Hs, band='terralook', red.f=red.f, green.f=green.f, blue.f=blue.f,
         mar=mar, axes=axes, cex=cex)
    mtext(sprintf("red.f=%g, green.f=%g, blue.f=%g", red.f, green.f, blue.f),
          side=3, line=0, adj=0, cex=cex)
    mtext(format(Hs[['time']], "%B %Y"), side=3, line=0, adj=1, cex=cex)
    if (!missing(name)) dev.off()
}

Now, we can try some tests. Leaving out name argument will produce interactive plots, and indeed that is how I selected the values I show below. Note that I have hard-coded the plot sizes, rather than use Jekyll and Rmarkdown to do that.

Case 1.

Present-day default, much too green overall. Look especially at the snow-covered citadel/commons region near 63.58W and 44.65N.

1
demo(2, 2, 4, "2016-02-20-landsat-1.png")

center

Case 2.

The greens are more reasonable now, e.g. still green in the evergreen trees of Point Pleasant Park, and along the railway cut line that crosses the city a little south of 44.63N. However, there a slight yellow-reddish hue overall (see especially the snow-covered citadel).

1
demo(2, 1.7, 4, "2016-02-20-landsat-2.png")

center

Case 3.

Decrease red and increase blue, which neutralize the winter hue (especially of the snow-covered citadel) and of summer paved (especially Pier 21 and railway funnel). The asphalt regions (e.g. the dockyard, year-round) are now a neutral gray. A blue hue covers water in both seasons. The summertime greens are rather bright, but this is what this tree-covered region looks like.

The lakes on the south side of the Northwest Arm are probably covered with smooth and rather clean snow, and so it is agreeable that the colour is both brighter there than on the citadel, which is, after all, a fort surrounded by a steep hill.

Overall, I am quite satisfied by the neutrality of both the ice/snow and darker asphalt/gravel regions. The greens also look good, and close inspection of the winter image will show greens in the parts of Halifax that have coniferous trees.

The main remaining question is whether the coefficients will handle red properly, because there are no prominent red items in this view. My experiments have shown that a neutralization of the gray/white regions can be achieved with a variety of red-blue pairings, and so I would like to have another test image that will help me to set that balance in a way that works for both the green-gray domain of Halifax and some other domain. If any reader wants to help with that, please contact me directly (since I do not have comments turned on for this blog).

1
demo(1.7, 1.7, 6, "2016-02-20-landsat-3.png")

center

Conclusions

These tests suggest that the values red.f=1.7, green.f=1.7 and blue.f=6 may be more reasonable as defaults than the present values of 2, 2, and 5, respectively. I will wait a week for any personal advice from colleagues, and consider switching the defaults after that time.

References and resources

issue report on landsat hue
Jekyll source code for this blog entry: 2016-02-20-landsat-hue.Rmd