The fldgen package allows you to ingest temperature output from an earth system model (ESM) and generate randomized temperature fields that have the same space and time correlation properties as the original ESM data. This tutorial focuses on how to use the functions in the package to generate and analyze temperature fields. The details of how the method works are covered in a companion paper. The statistics and analysis from that paper are used as the case study in this tutorial.

All of the functions used here are documented in R’s help system. Since our purpose here is to outline what functions have to be called, and in what sequence, to perform the analysis, we haven’t repeated material from the help files. If you’re confused about how a function is supposed to work, consult the help files. For example, help(read.temperatures) will print the docs for the function that reads the netCDF temperature fields.

Software Version

These calculations were run with fldgen-v1.0.0

Setup

These are the parameters we will be using for the calculations below.

## parameters for the code below.
ngen <- 9             # number of fields to generate
nplt <- 4             # number of fields to plot
exfld <- 20           # example field to plot from the time series
set.seed(8675309)     # Set RNG seed so results will be reproducible
siglvl <- 0.05        # significance level for statistical tests
method <- 1           # method to use in mkcorrts

Reading the ESM data

The ESM temperature field should be in netCDF files. If the files are all in a single directory (and there are no extraneous files in that directory), then all we need is the directory name.

datadir <- './data'
emulator <- train(datadir, latvar='lat_2', lonvar='lon_2')
griddata_l <- splitGrids(emulator$griddata)

Generating fields

The total power in each EOF mode is proportional to its variance. The modes with higher power are more important to determining the field structure; the ones with lower power are less important.

reof <- emulator$reof
rawpwr <- reof$sdev[-1]   # Exclude PC0, which has almost no power in it.
relpwr <- (rawpwr / max(rawpwr))^2
pltpwr <- ggplot(data=data.frame(total.power=relpwr, EOF=seq_along(relpwr)-1), 
                 aes(x=EOF, y=total.power)) + geom_col() + ylab('Relative Power') + theme_minimal()
pltlogpwr <- ggplot(data=data.frame(log.power=log(relpwr)-log(min(relpwr)), EOF=seq_along(relpwr)-1),
                    aes(x=EOF, y=log.power)) + geom_col()

print(pltpwr)

#print(pltlogpwr)

p1_10 <- round(sum(relpwr[2:11])/sum(relpwr), 2) * 100
p1_50 <- round(sum(relpwr[2:51])/sum(relpwr), 2) * 100
ptail <- round(sum(relpwr[400:length(relpwr)])/sum(relpwr), 2) * 100
#p1_10
#p1_50
#ptail
#min(relpwr[-1])

Two things are apparent. First, most of the power is in the first few dozen components. The first 10 components contain 49% of the total power, while the first 50 contain 72%. However, there is a long tail of components that have low power individually, but which together contribute nontrivially to the total. Even though the contributions of the components past number 400 are all but invisible on the plot, together they contribute 1% of the total power.

The last thing we need to generate our fields is the temporal structure of the EOF coefficients. We get this from the power spectral density (PSD) of the coordinates of the residuals in the coordinate system defined by the EOF basis vectors. In order to do this, we have to split the eof structure, so that we can analyze the time series from each ESM run separately.

Fxmag <- emulator$fx$mag

Now we’re ready to generate fields. We’ll generate 9 fields in this example, to get valid correlation statistics, but we’ll only plot 4 of them. We need tgav to be the same length as the time series we used to generate the power spectrum. In practical terms, this means taking the tgav values from one of the input ESM simulations (or those generated from another process, such as an SCM.) For this example we’ll use one of the RCP6 runs

esmscen <- 'tas_annual_CESM1-CAM5_rcp60_r2i1p1_200601-210012.nc'   # ESM scenario to use
griddata <- emulator$griddata
tgav <- emulator$tgav
pscl <- emulator$pscl

strt <- griddata$tags[[esmscen]][1]
end <- griddata$tags[[esmscen]][2]
tgreconst <- tgav[strt:end]
meanfield <- pscl_apply(pscl, tgreconst)
tempgrids <- lapply(1:ngen, 
                    function(x) {
                        bcoord <- mkcorrts(emulator, method=method)   # Because this uses a RNG, it will be different every call.
                        reconst_fields(reof$rotation, bcoord, meanfield)
                    })

## Subtract off the mean field.  Save these because we will want to use them later
residgrids <- lapply(tempgrids, function(g) {
    g - meanfield
})

## We will also want a grand matrix of all of the generated residual fields to use for statistics.  
allresidgrids <- do.call(rbind, residgrids)

Plotting the fields

We can extract a single field from each time series and plot them all for comparision. We will be able to see a lot more detail if we plot the residuals instead of the total temperature values.

plotglobalfields <- function(fieldlist, sequencenum, minval=-3.5, maxval=3.5, legendstr) {
    ## Extract a single example field from each series (sequencenum determines which one) and create a plot
    lapply(fieldlist, function(g) {
        suppressWarnings(
            plot_field(g[sequencenum,], griddata, 14,  # 14 is the number of color levels in the plot
                       minval, maxval) + 
                guides(fill=ggplot2::guide_colorbar(barwidth=ggplot2::unit(4,'in'), barheight=ggplot2::unit(0.125,'in'),
                                                    title=legendstr, title.position='top', title.hjust=0.5))
        )
    })
}

residfieldplots <- plotglobalfields(residgrids[1:nplt], exfld, legendstr='delta-T (K)')
#> Loading required namespace: gcammaptools

## Display the plots
for (plt in residfieldplots) {
    print(plt)
}

For completeness, here are the full-temperature plots at the same time.


tempfieldplots <- plotglobalfields(tempgrids[1:nplt], exfld, 220, 307, legendstr = 'Temperature (K)')

## Display the plots
for (plt in tempfieldplots) {
    print(plt)
}

Differences among these fields are subtler because the residual field’s range of +/- a few degrees is much smaller than the range of temperatures represented in the mean field.

Analysis

Global mean temperatures

The global mean temperature should be close to zero for all of the residual fields, but it probably won’t be exact. The differences here tell us something about the performance of the mean response algorithm. If the global mean of the mean response field for input temperature \(T_g\) is something other than \(T_g\), then the algorithm isn’t consistent.

meantemps <- lapply(seq_along(residgrids), function(i) {
    residgrids[[i]] %*% griddata$tgop
    }) 
for(mt in meantemps)
    print(as.vector(summary(mt)))
#> [1] "Min.   :-3.968e-12  " "1st Qu.:-1.051e-12  " "Median : 1.317e-13  "
#> [4] "Mean   : 1.470e-13  " "3rd Qu.: 1.501e-12  " "Max.   : 4.032e-12  "
#> [1] "Min.   :-4.943e-12  " "1st Qu.:-1.229e-12  " "Median : 1.142e-13  "
#> [4] "Mean   :-1.470e-13  " "3rd Qu.: 1.047e-12  " "Max.   : 3.626e-12  "
#> [1] "Min.   :-4.097e-12  " "1st Qu.:-1.097e-12  " "Median : 1.095e-13  "
#> [4] "Mean   : 1.471e-13  " "3rd Qu.: 1.293e-12  " "Max.   : 4.320e-12  "
#> [1] "Min.   :-3.480e-12  " "1st Qu.:-1.518e-12  " "Median : 1.548e-13  "
#> [4] "Mean   : 1.471e-13  " "3rd Qu.: 1.418e-12  " "Max.   : 5.192e-12  "
#> [1] "Min.   :-3.367e-12  " "1st Qu.:-1.133e-12  " "Median :-7.228e-14  "
#> [4] "Mean   : 1.471e-13  " "3rd Qu.: 1.437e-12  " "Max.   : 6.120e-12  "
#> [1] "Min.   :-5.576e-12  " "1st Qu.:-1.257e-12  " "Median :-4.689e-13  "
#> [4] "Mean   :-1.470e-13  " "3rd Qu.: 9.232e-13  " "Max.   : 6.370e-12  "
#> [1] "Min.   :-6.623e-12  " "1st Qu.:-9.121e-13  " "Median : 3.937e-13  "
#> [4] "Mean   : 1.470e-13  " "3rd Qu.: 1.411e-12  " "Max.   : 3.751e-12  "
#> [1] "Min.   :-4.395e-12  " "1st Qu.:-1.256e-12  " "Median : 7.038e-14  "
#> [4] "Mean   : 1.470e-13  " "3rd Qu.: 1.607e-12  " "Max.   : 4.809e-12  "
#> [1] "Min.   :-4.972e-12  " "1st Qu.:-1.329e-12  " "Median : 5.387e-14  "
#> [4] "Mean   :-1.470e-13  " "3rd Qu.: 1.205e-12  " "Max.   : 3.748e-12  "

Grid cell variance

We also expect the variance in each grid cell to be the same as the variance in the original ESM data. We can test this by calculating the \(F\) statistic. We don’t want to do an F-test on each grid cell, and some of them are expected to fail anyhow, since there are grid cells. Instead, we’ll calculate all the \(F\) values and see how many would fail at the 0.05 level.

## variance in each cell in the observations
obscellvar <- apply(pscl$r, 2, var)

## variance in each cell in the reconstructions
cellvar <- lapply(residgrids, function(g) {
    apply(g, 2, var)   # calculate variance in each grid cell across time slices.
    })

fstats <- lapply(cellvar, function(cvg) {
    cvg / obscellvar     # F = var(x) / var(y)
    })

pv <- c(siglvl/2, 1-siglvl/2)      # This is a two-tailed test, so the p-vals are half the sig level
df1 <- nrow(residgrids[[1]]) - 1   # degrees of freedom is N-1; all grids have the same N
df2 <- nrow(pscl$r) - 1
fq <- qf(pv, df1, df2)
fracfail <- sapply(fstats, function(f) {
    sum(f<fq[1] | f>fq[2]) / length(f)
})
print(fracfail)
#> [1] 0.03776042 0.04788773 0.05720124 0.06913701 0.05311415 0.04163050
#> [7] 0.02718099 0.04673032 0.03020110

In other words the number of grid cells that appear to show significant differences in variance is about what we would expect by chance alone (i.e., roughly 1% of each grid).

An alternate version of this calculation is to look at the variances of the grid cells in all of the field realizations, treating the variance as an operation across both time and alternate realizations. We’ll do that here and see what if any difference it makes.

fstats.all <- apply(allresidgrids, 2, var) / obscellvar
df.all <- nrow(allresidgrids) - 1    # More degrees of freedom in the merged version.
fq.all <- qf(pv, df.all, df2)          # Original dataset still has same number of df.
fracfail.all <- sum(fstats.all < fq.all[1] | fstats.all > fq.all[2]) / length(fstats.all)
fflo <- sum(fstats.all < fq.all[1]) / length(fstats.all)
ffhi <- sum(fstats.all > fq.all[2]) / length(fstats.all)
print(fracfail.all)
#> [1] 0.003291377
print(fflo)
#> [1] 0.003273293
print(ffhi)
#> [1] 1.808449e-05

## Calculate the power of this test
pwr <- power.f(df.all, df2, 0.01, siglvl)
print(pwr)
#> [1] 0.05241955

The fraction that fail is even lower in this version. I’m still trying to decide if there is any important information here, or if it’s just telling us that our method is producing the appropriate variance more consistently than a purely random process would.

Grid cell normality

We also need to establish that the data are normally distributed. (Technically we should do this before we apply the F-test.) We’ll use the Shapiro-Wilk statistic for this test. We don’t have a separate quantile function for null hypothesis distribution of this statistic, so we’ll just use the p-values produced by shapiro.test.

swtestpvals <- lapply(residgrids, function(g) {
    apply(g, 2, function(x){
        shapiro.test(x)$p.value
    })
})

fracfailsw <- sapply(swtestpvals, function(pv) {
    sum(pv < siglvl) / length(pv)
})
print(fracfailsw)
#> [1] 0.06421803 0.09045862 0.02846499 0.05253545 0.02839265 0.04633247
#> [7] 0.04255281 0.04486762 0.08416522

Looking at the results of the Shapiro-Wilk test, we find about the expected number of failures, indicating that the data are consistent with a normal distribution.

This is what we would expect, even if the original data were not normally distributed because the the process is essentially summing up a collection of random variables, and the Central Limit Theorem tells us that such sums approach a normal distribution as the number of terms in the sum increases.

As with the variance, we can do this test on the merged dataset.

swpvals.all <- apply(allresidgrids, 2, function(x) {
    shapiro.test(x)$p.value
})

fracfailsw.all <- sum(swpvals.all < siglvl) / length(swpvals.all)
print(fracfailsw.all)
#> [1] 0.06976997

Roughly 7% of the grid cells fail the test. This is a little on the high side, but not worryingly so. I believe that what we’re seeing here is that the Central Limit Theorem says we approach a normal distribution in the limit of an infinite number of summed values. We’re only summing values in batches of 95, so it’s not surprising that if we collect enough of them we can start to detect some deviations from a strict normal distribution. I’m ok with this.

Correlation of projection coefficients

We expect the projection coefficients for each EOF to be uncorrelated with the projection coefficients for other EOFs. With 854 EOFs, this is not actually possible in a single realization, so this test needs to be averaged across all of the realizations.

phi <- allresidgrids %*% emulator$reof$rotation
iscorr <- function(i,j) {
    cor.test(phi[,i], phi[,j], method='spearman', alternative='two.sided')$p.value < siglvl
}

## skip i=1 in this test because it is EOF-0 and isn't expected to follow this rule
neof <- ncol(phi)-1
ncor <- sum(
    sapply(seq(2, ncol(phi)-1),        
           function(i) {
               sum(simplify2array(Map(iscorr, i, seq(i+1, ncol(phi)))))
         })
)
## fraction correlated
fcor <- ncor / (neof*(neof-1)/2)
fcor
#> [1] 0.05056718

The fraction of the EOF pairs that test positive for correlation is 5.1%, which is right about what we expect for the chosen significance level of 0.05 for our test. I looked up the statistical power of this test in a table, and for the 855 samples we have here, the power of the test is about 0.64 for a presumed correlation of 0.1, and it is 0.13 for a presumed correlation of 0.05. Therefore, from this test we can be fairly confident that the correlation between projection coefficients for distinct EOFs is less than 0.05.

EOF Characteristics

Time behavior

Start by constructing a heatmap of the power spectrum for the EOFs.

nt <- length(griddata_l[[1]]$time)
## There is no additional information in the negative frequencies, so keep only
## the positive ones.
np <- 
    if(nt %% 2 == 1) { 
        (nt+1)/2
    } else {
        nt/2 + 1
    }

psd <- emulator$fx$mag^2
hmdata <- reshape2::melt(psd[1:np,])
hmdata$freq <- (hmdata$Var1 - 1) / nt    # The - 1 is due to R's unit-indexed arrays.
hmdata$EOF <- as.integer(substring(hmdata$Var2, 3))
## discretize the power so we can isolate structure more easily
nbrk <- 10
hmdata$discval <- findInterval(hmdata$value / max(hmdata$value), seq(0.01, 0.99, length.out=nbrk)) / nbrk
hmdata <- dplyr::select(hmdata, EOF, freq, value, discval) %>% as_tibble()

hmplt <- ggplot(hmdata, aes(x=freq, y=EOF, fill=discval)) + geom_raster() + scale_fill_distiller(palette='YlOrRd', direction=1, name='Relative Power')

hmdata50 <- filter(hmdata, EOF <= 50)
hmplt50 <- ggplot(hmdata50, aes(x=freq, y=EOF, fill=discval)) + geom_raster() + scale_fill_distiller(palette='YlOrRd', direction=1, name='Relative Power')

print(hmplt)

print(hmplt50)

Just as we saw above, the power drops off quite a bit after the first few EOFs.

More importantly, the power spectrum whitens after the first few EOFs, so those early EOFs represent periodic signals, while the later ones don’t. Here are the smoothed power spectra for the first 9 EOFs.

eofspectra <- filter(hmdata, EOF>0, EOF<10) %>% mutate(EOF=factor(EOF), value=value/max(value))
eofspectraplt <- ggplot(eofspectra, aes(x=freq, y=value, color=EOF)) + geom_smooth(se=FALSE) + scale_color_brewer(palette='Set1') + xlab('freq (1/yr)') + ylab('Relative Spectral Power Density') + theme_minimal()
print(eofspectraplt)
#> `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Spatial behavior

We will make spatial plots of the 9 EOFs shown in the plot above. We’ll also plot EOFs 25 and 50, just to get an idea of what’s happening in those lower power modes.

### Plotting global maps is still a little slow, so expect this to take some time.
## The EOFs are in reof$rotation.  Each column is the grid cell values for an EOF.
## Also, the EOFs are scaled to unit norm.  We'll rescale them to unit max value

eofcols <- c(2:10, 26, 51, 201, 401, 601, 854)  # EOF numbering starts at 0, but array numbering starts at 1
eofvis <- t(reof$rotation[,eofcols])   # EOFs are now in rows, not columns
eofvis <- eofvis / max(abs(eofvis))

eofplts <- lapply(seq_along(eofcols), function(i) {
    title <- paste0('EOF-', eofcols[i]-1)
    suppressWarnings(plot_field(eofvis[i,], griddata, 14, -0.5, 0.5)) + 
        ggtitle(title) +
        guides(fill=ggplot2::guide_colorbar(barwidth=ggplot2::unit(4,'in'), barheight=ggplot2::unit(0.125,'in'),
                                                  title="delta-T (K)", title.position='top', title.hjust=0.5))
})
    
for(i in seq_along(eofcols)) {
    print(eofplts[[i]])
}