Chapter 2 bigmemory and ff
2.1 Regular data example
For testing, load the 1\(^\circ\) World Ocean Atlas (Levitus) sea surface temperature (SST) data from Locarnini et al. (2013) and sea surface salinity (SSS) data from Zweng et al. (2013) data via the oceadata package (D. Kelley 2019).
library(ocedata)
data(levitus, package="ocedata")
str(levitus) # the data is filled with NA at land
## List of 4
## $ longitude: num [1:360] -180 -178 -178 -176 -176 ...
## $ latitude : num [1:180] -89.5 -88.5 -87.5 -86.5 -85.5 -84.5 -83.5 -82.5 -81.5 -80.5 ...
## $ SSS : num [1:360, 1:180] NA NA NA NA NA NA NA NA NA NA ...
## $ SST : num [1:360, 1:180] NA NA NA NA NA NA NA NA NA NA ...
Based on this SST data, create bigmemory (Kane et al. 2018) and ff (Adler et al. 2018) objects:
sst_base <- levitus$SST # the default format from base package
library(bigmemory)
sst_bm <- bigmemory::as.big.matrix(levitus$SST) # bigmemory format
library(ff)
sst_ff <- ff::as.ff(levitus$SST) # ff format
c(identical(sst_base, sst_bm[,]), # note the different index syntax [,]
identical(sst_base, sst_ff[,]))
## [1] TRUE TRUE
str(sst_base)
## num [1:360, 1:180] NA NA NA NA NA NA NA NA NA NA ...
str(describe(sst_bm))
## Formal class 'big.matrix.descriptor' [package "bigmemory"] with 1 slot
## ..@ description:List of 12
## .. ..$ sharedType: chr "SharedMemory"
## .. ..$ sharedName: chr "0f2964c1-361c-47b3-8765-225d4e98647b"
## .. ..$ totalRows : int 360
## .. ..$ totalCols : int 180
## .. ..$ rowOffset : num [1:2] 0 360
## .. ..$ colOffset : num [1:2] 0 180
## .. ..$ nrow : num 360
## .. ..$ ncol : num 180
## .. ..$ rowNames : NULL
## .. ..$ colNames : NULL
## .. ..$ type : chr "double"
## .. ..$ separated : logi FALSE
str(virtual(sst_ff))
## List of 4
## $ Length : int 64800
## $ Dim : int [1:2] 360 180
## $ Dimorder : int [1:2] 1 2
## $ Symmetric: logi FALSE
The sizes of the three different SST objects highlight the large difference between the standard base compared to the bigmemory and ff packages in terms of memory usage (see my adapted form of Dirk Eddelbuettel’s lsos() function here):
df <- myls(pattern="sst_(base|ff|bm)", show=F)
library(formattable)
ft <- formattable::formattable(df$ls, list(
RelSize=formattable::formatter("span",
style=function(x) {
formattable::style(display="inline-block",
direction="ltr",
"border-radius"="4px",
"background-color"="#b8dff3",
width=formattable::percent(x))
})))
formattable::as.datatable(ft, options=list(dom='t'), class="display nowrap")
Manipulations of the bigmemory and ff objects can be obtained as with base. For example, the area-averaged SST can be calculated for all data points within or on the edge of an arbitrary polygon (see Fig. 2.1).
library(sp) # for point.in.polygon()
poly <- list(x=c(-30, -75, -30, -15, -30), # deg longitude
y=c(60, 30, 0, 30, 60)) # deg latitude
lon <- levitus$longitude
lat <- levitus$latitude
xy <- expand.grid(lon, lat, KEEP.OUT.ATTRS=F)
colnames(xy) <- c("x", "y")
inds <- sp::point.in.polygon(xy$x, xy$y, poly$x, poly$y) # returns indices as vector, not matrix
inds <- which(inds == 1 | inds == 2) # inside or on edge of polygon1800 (length(inds)) out of 64800 (dim(xy)[1]) 1\(^\circ\) data points are located within or on the edge of poly (note that in this example the indices would be the same if only the ones completely within were considered, i.e. inds == 1). Fig. 2.1 shows the SST data and the obtained points in orthographic projection and was realized with the oce package (D. Kelley and Richards 2018). More information on available projections can be found by vignette("oce"), ?oce::mapPlot or in Appendix C of D. E. Kelley (2018). The plot was generated like this:
library(oce) # for mapPlot() etc.
library(fields) # for colorbar via image.plot()
# make colorbar
zlim <- range(sst_base, na.rm=T)
breaks <- pretty(zlim)
cm <- oce::colormap(breaks=breaks,
col=oce::oce.colorsPalette(length(breaks)-1))
# define projection
p <- "+proj=ortho +lat_0=30 +lon_0=-45" # orthographic projection
# open plot device and set margins
#par(mar=c(0, 0, 0, 5), oma=c(0, 0, 0, 0)) # larger right margin for colorbar
par(mar=c(0, 0, 0, 0), oma=c(0, 0, 0, 0)) # larger right margin for colorbar
# set coordinate system
oce::mapPlot(xy$x, xy$y, projection=p,
grid=T,
#type="n",
longitudelim=c(-95, -5),
latitudelim=c(20, 90),
axes=F, drawBox=F)
# color everything in gray for land
oce::mapImage(levitus$longitude, levitus$latitude,
array(1, dim(sst_base)), col="gray66")
# add the data with colors
oce::mapImage(levitus$longitude, levitus$latitude,
sst_base, breaks=cm$breaks,
filledContour=F, gridder=NA)
# add breaks as contour lines
oce::mapContour(levitus$longitude, levitus$latitude,
sst_base, levels=cm$breaks, lwd=0.75)
# add grid lines
oce::mapGrid(dlongitude=15, dlatitude=10, polarCircle=5,
col="black", lwd=0.5)
# add our arbitrary polygon (project piece-wise)
for (i in 1:(length(poly$x) - 1)) {
tmp <- oce::lonlat2map(seq(poly$x[i], poly$x[i+1], l=100),
seq(poly$y[i], poly$y[i+1], l=100),
projection=p)
lines(tmp, lwd=1.5)
} # for i points
# add points within/on the edge of our arbitrary polygon
tmp <- oce::lonlat2map(xy$x[inds], xy$y[inds], projection=p)
points(tmp, pch=".")
# add colorbar
fields::image.plot(zlim=zlim, legend.only=T, lwd=0.5,
breaks=cm$breaks, col=cm$col,
legend.mar=par("mar")[4] + 1)
mtext(text="SST [°C]", side=4, line=par("mar")[4] - 2)![Average (1955-2012) SST [in °C; colors; no data at gray areas; @locarnini_etal_2013] plotted in orthographic projection with the `oce` package [@R-oce]. Black dots mark the center of all data matrix elements within or on the edge of the arbitrary polygon `poly` shown by the thick black line.](figs/sst-ortho.png)
Fig. 2.1: Average (1955-2012) SST [in °C; colors; no data at gray areas; Locarnini et al. (2013)] plotted in orthographic projection with the oce package (D. Kelley and Richards 2018). Black dots mark the center of all data matrix elements within or on the edge of the arbitrary polygon poly shown by the thick black line.
A subset of the data within the polygon can be obtained from the bigmemory and ff objects by integer indexing as for the default base object. Remember that the returned indices by sp::point.in.polygon() are a vector of nlon*nlat length, not a nlon \(\times\) nlat dimensioned matrix. That means the data needs to be converted from a 2D matrix to a 1D vector before selecting elements via inds.
sst_vec_base <- as.vector(sst_base)
sst_sub_base <- sst_vec_base[inds]
sst_vec_bm <- bigmemory::as.big.matrix(sst_vec_base)
sst_sub_bm <- bigmemory::as.big.matrix(sst_vec_bm[][inds])
sst_vec_ff <- ff::as.ff(sst_vec_base)
sst_sub_ff <- sst_vec_ff[ff::as.ff(inds)]
c(identical(sst_sub_base, sst_sub_bm[]),
identical(sst_sub_base, sst_sub_ff[]))
## [1] TRUE TRUE
Now, the average SST within the polygon can be calculated by, for example, the area \(A\) -weighted arithmetic mean \(\mu_\text{SST} = (\sum_i A_i)^{-1} \, \sum_i (A_i \cdot \text{SST}_i)\), with \(i = (1, \dots, 1800)\):
library(SDMTools) # for grid.info()
area_m2 <- grid.info(lats=lat, cellsize=diff(lon)[1])$area
# bring on same format as the y dimension of the xy vector
area_m2 <- rep(area_m2, e=length(lon))
# select within/on the edge of our arbitrary polygon
area_m2 <- area_m2[inds]
# set potential NA values also in the area_m2 vector
if (any(is.na(sst_sub_base))) {
area_m2[is.na(sst_sub_base)] <- NA
}
sst_mean_base <- sum(sst_sub_base*area_m2, na.rm=T)/sum(area_m2, na.rm=T)
sst_mean_bm <- sum(sst_sub_bm[]*area_m2, na.rm=T)/sum(area_m2, na.rm=T)
sst_mean_ff <- sum(sst_sub_ff[]*area_m2, na.rm=T)/sum(area_m2, na.rm=T)
# todo: esd::aggregate.area(x,is=NULL,it=NULL,FUN=’sum’,na.rm=TRUE,smallx=FALSE)
c(sst_mean_base, sst_mean_bm, sst_mean_ff)
## [1] 22.03309 22.03309 22.03309
c(identical(sst_mean_base, sst_mean_bm),
identical(sst_mean_base, sst_mean_ff))
## [1] TRUE TRUECompare this result to cdo fldmean:
# save the sst field as nc file
lon_dim <- ncdim_def(name="lon", units="degrees_east", vals=lon)
lat_dim <- ncdim_def(name="lat", units="degrees_north", vals=lat)
sst_var <- ncvar_def(name="sst", units="degC", dim=list(lon_dim, lat_dim),
missval=NA, prec="double")
outnc <- nc_create("data/sst.nc", vars=sst_var, force_v4=T)
ncvar_put(outnc, sst_var, sst_base)
nc_close(outnc)
# save our arbitrary polygon
write.table(poly, file="data/poly.txt", row.names=F, col.names=F)
cat(system("cat data/poly.txt", intern=T), sep="\n")
## -30 60
## -75 30
## -30 0
## -15 30
## -30 60
# apply the polygon and calculate the cdo field mean
system("cdo -maskregion,data/poly.txt data/sst.nc data/sst_mask.nc")
system("rm data/sst_mean_cdo.nc; cdo -fldmean data/sst_mask.nc data/sst_mean_cdo.nc")
# load the result in R workspace
ncin <- nc_open("data/sst_mean_cdo.nc")
sst_mean_cdo <- ncvar_get(ncin, "sst")
sst_mean_cdo
## [1] 22.03309
identical(sst_mean_base, sst_mean_cdo)
## [1] FALSEWhops, the R and cdo results are not the same? They are, but differ after the 12th decimal postion:
print(sst_mean_base, digits=20)
## [1] 22.033094656022573332
print(sst_mean_cdo, digits=20)
## [1] 22.033094656022786495
sst_mean_base - sst_mean_cdo
## [1] -2.131628e-13In conclusion, for the tested base, bigmemory and ff objects, the area-averaged SST within the polygon (total area \(\sim\) 19 \(\times\) 106 km2; Fig. 2.1) \(\mu_\text{SST}\) = 22.0330947 \(^\circ\)C. The same result is obtained with cdo fldmean.