A short script to visualize and explore the updated Land Surface Climatology algorithm that 'lowers the standards' in some areas to increase the number of available observations.
download = F
if (download) system("wget -e robots=off -L -r -np -nd -p 20140304_LST -nc -A tif http://ecocast.arc.nasa.gov/data/pub/climateLayers/LST_new/")
## organize file names
f = data.frame(full = T, path = list.files("20140304_LST", pattern = "tif$",
full = T), stringsAsFactors = F)
f$month = as.numeric(do.call(rbind, strsplit(basename(f$path), "_|[.]"))[, 7])
f$type = do.call(rbind, strsplit(basename(f$path), "_|[.]"))[, 1]
f = f[order(f$month), ]
f$mn = month.name[f$month]
## create raster stacks
lst_mean = stack(f$path[f$type == "mean"])
names(lst_mean) = f$mn[f$type == "mean"]
NAvalue(lst_mean) = 0
lst_nobs = stack(f$path[f$type == "nobs"])
names(lst_nobs) = f$mn[f$type == "nobs"]
lst_qa = stack(f$path[f$type == "qa"])
names(lst_qa) = f$mn[f$type == "qa"]
## define a function to summarize data
fst = function(x, na.rm = T) c(mean = mean(x, na.rm = T), min = min(x, na.rm = T),
max = max(x, na.rm = T))
rfst = function(r) cellStats(r, fst)
Map of LST by month (with white indicating missing data). Note that many inland regions have missing data (white) in some months (mostly winter).
colramp = colorRampPalette(c("blue", "orange", "red"))
dt_mean = rfst(lst_mean)
levelplot(lst_mean, col.regions = c(colramp(99)), at = seq(0, 65, len = 99),
main = "Mean Land Surface Temperature", sub = "Tile H08v05 (California and Northern Mexico)")
Table of mean, min, and maximum LST for this tile by month.
print(xtable(dt_mean), type = "html")
| January | February | March | April | May | June | July | August | September | October | November | December | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| mean | 15.01 | 18.41 | 25.81 | 32.05 | 39.49 | 44.18 | 44.70 | 41.96 | 38.59 | 30.84 | 21.77 | 14.33 |
| min | 1.43 | 1.98 | 1.08 | 1.28 | 2.37 | 5.03 | 12.04 | 12.96 | 12.72 | 6.90 | 2.77 | 2.76 |
| max | 31.31 | 36.74 | 45.88 | 52.29 | 58.66 | 60.88 | 61.15 | 60.99 | 57.07 | 48.82 | 39.17 | 29.42 |
lst_tmean = melt(unlist(as.matrix(lst_mean)))
colnames(lst_tmean) = c("cell", "month", "value")
lst_tmean = lst_tmean[!is.na(lst_tmean$value), ]
lst_tmean$month = factor(lst_tmean$month, levels = month.name, ordered = T)
bwplot(value ~ month, data = lst_tmean, ylab = "Mean LST (c)", xlab = "Month")
This section details the spatial and temporal distribution of the number of LST observations that were not masked by quality control (see section below). Note that the regions with no data in the map above have missing data (nobs=0) here as well, but also the areas surrounding those regions have low numbers of observations in some months (blue colors).
dt_nobs = rfst(lst_nobs)
levelplot(lst_nobs, col.regions = c("grey", colramp(99)), at = c(-0.5, 0.5,
seq(1, 325, len = 99)), main = "Sum Available Observations", sub = "Tile H08v05 (California and Northern Mexico) \n Grey indicates zero observations")
Table of mean, min, and maximum number of observations for this tile by month.
| January | February | March | April | May | June | July | August | September | October | November | December | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| mean | 107.97 | 100.06 | 140.35 | 152.19 | 177.02 | 178.22 | 156.86 | 169.87 | 180.29 | 167.68 | 136.12 | 102.86 |
| min | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |
| max | 272.00 | 238.00 | 281.00 | 293.00 | 319.00 | 304.00 | 319.00 | 324.00 | 310.00 | 305.00 | 271.00 | 260.00 |
The seasonal cycle of missing data is quite noisy, though there tend to be fewer observations in winter months (DJF).
lst_tnobs = melt(unlist(as.matrix(lst_nobs)))
colnames(lst_tnobs) = c("cell", "month", "value")
lst_tnobs = lst_tnobs[!is.na(lst_tnobs$value), ]
lst_tnobs$month = factor(lst_tnobs$month, levels = month.name, ordered = T)
bwplot(value ~ month, data = lst_tnobs, ylab = "Number of availble observations",
xlab = "Month")
Map of the Quality Assessment (QA) level used to fill the pixels. It goes from 0 (highest quality) to 3(low). For h08v05 all pixels are filled with either 0 or 1. So red indicates areas with the lower quality data (most of the tile).
levelplot(lst_qa, col.regions = c("grey", "red"), at = c(-0.5, 0.5, 1.5), cuts = 2,
main = "Quality Assessment Filter", sub = "Tile H08v05 (California and Northern Mexico)")
Proportion of cells in each month with QA=1 (including cells in the Pacific Ocean)
| January | February | March | April | May | June | July | August | September | October | November | December | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| ProportionQA_1 | 0.47 | 0.47 | 0.42 | 0.40 | 0.38 | 0.37 | 0.42 | 0.40 | 0.38 | 0.38 | 0.43 | 0.47 |
A few open questions/comments (in my mind):