Gridpp in R

[Under development]

Overview
Structure Functions Overview
Setup for OI examples
OI example: Combining Multiple Structure Functions

Overview

The package can be loaded as follows (from the build/swig/R folder):

dyn.load("gridpp.so")
source("gridpp.R")
cacheMetaData(1)

NOTE: there is a bug in SWIG Versions 4.0.2 and 4.2.1 (we have not tested other versions), it creates a build/swig/R/gridpp.R that is not correct. A quick fix is:

cd build/swig/R
for i in {1..50}; do sed -i "s/all(sapply(argv\[\[$i\]\] , is.integer) || sapply(argv\[\[$i\]\], is.numeric))/all(sapply(argv\[\[$i\]\] , is.integer) | sapply(argv\[\[$i\]\], is.numeric))/g" gridpp.R; done

Functions can be called in the same way as in python. Objects can also be created in the same way. However, class member functions are called in a different way.

structure <- BarnesStructure(10000, 0)
p1 <- Point(0, 0)
p2 <- Point(0, 0.1)
correlation <- structure$corr(p1, p2)
print(correlation) 

Alternatively, for Cartesian coordinates:

structure <- BarnesStructure(10000, 0)
p1 <- Point(0, 0, 0, 0, "Cartesian")
p2 <- Point(30000, 0, 0, 0, "Cartesian")
correlation <- structure$corr(p1, p2)
print(correlation) 

Structure Functions Overview

There are five main analytical forms for defining structure functions in gridpp, all of which depend on the distance between two points ($d$) and a scaling length parameter ($L$): Barnes, SOAR, TOAR, Powerlaw, and Cressman.

• Barnes (Barnes, 1973): This function is of the form $e^{-x^2}$, where $x = d / L$. It resembles a Gaussian function. For $L = 1$, a correlation value of 0.0013 is reached at $d = 3.65$.
• SOAR (Second-order autoregressive function, Gasparri and Cohn, 1999): Defined as $(1 + x)⋅e^{-x}$, with $x = d / L$. For $L = 1$, a correlation value of 0.0013 is reached at $d = 8.54$.
• TOAR (Third-order autoregressive function, Gasparri and Cohn, 1999): Has the form $(1 + x + 1/3 ⋅ x^2)⋅e^{-x}$, with $x = d / L$. For $L = 1$, a correlation value of 0.0013 is reached at $d = 9.49$.
• Powerlaw (Gasparri and Cohn, 1999): This function is defined as $(1 + 1/2 ⋅ x^2)^{-1}$, with $x = d / L$. For $L = 1$, a correlation value of 0.0013 is reached at $d = 39.20$.
• Cressman: The function takes the form $(L^2 - d^2) / (L^2 + d^2)$. For $L = 1$, a correlation value of 0.0013 is reached at $d = 1.002$.

The sixth structure function in GridPP is a special linear function, designed for variables with values strictly between 0 and 1, inclusive. A common use case is for variables like land area fraction, which represents the fraction of land surrounding a point. Unlike typical correlation functions, this function returns a damping factor (a value smaller or equal to 1) rather than a correlation, allowing the correlation between two points to be adjusted based on differences in their land area fractions, for example. In this context, $d$ represents the difference in land area fraction between two points, and the function is defined as: $1 − (1−min)⋅∣d∣$ where $min$ is the minimum value of the damping coefficient. If two points have identical land area fractions, the damping coefficient is 1. If the difference is maximal (e.g., one point has a fraction of 0 and the other 1), the damping coefficient is set to $min$. This structure function allows users to reduce the correlation between points with differing land area fractions while maintaining some level of correlation, particularly if the points are spatially close.

The R-code used to obtain the figure above is reported in the following:

coord  <- seq(0,10,by=0.1)
ncoord <- length(coord)
structure_Barnes   <- BarnesStructure(1, 0, 0)  # Example: horizontal structure
structure_Soar     <- SoarStructure(1, 0, 0)  # Example: horizontal structure
structure_Toar     <- ToarStructure(1, 0, 0)  # Example: horizontal structure
structure_Powerlaw <- PowerlawStructure(1, 0, 0)  # Example: horizontal structure
structure_Cressman <- CressmanStructure(1, 0, 0)  # Example: horizontal structure
correlation_Barnes   <- vector( mode="numeric", length=ncoord)
correlation_Soar     <- vector( mode="numeric", length=ncoord)
correlation_Toar     <- vector( mode="numeric", length=ncoord)
correlation_Powerlaw <- vector( mode="numeric", length=ncoord)
correlation_Cressman <- vector( mode="numeric", length=ncoord)
p1 <- Point(0, 0, 0, 0, "Cartesian")
for (i in 1:ncoord) {
  x <- coord[i]
  p2 <- Point(x, 0, 0, 0, "Cartesian")
  correlation_Barnes[i]   <- structure_Barnes$corr(p1, p2)
  correlation_Soar[i]     <- structure_Soar$corr(p1, p2)
  correlation_Toar[i]     <- structure_Toar$corr(p1, p2)
  correlation_Powerlaw[i] <- structure_Powerlaw$corr(p1, p2)
  correlation_Cressman[i] <- structure_Cressman$corr(p1, p2)
}
png(file="correlations.png",width=800,height=600)
par(mar=c(5,5,1,1))
plot(coord,correlation_Barnes, ylim=c(0,1),xlim=c(0,10),col="white",cex.axis=2,xlab="Spatial 1D Coordinate", ylab="Correlation", cex.lab=2)
abline(v=0:100,h=seq(0,1,by=0.1),lty=2,col="gray")
abline(v=0,h=c(0,1),lty=1,col="darkgray")
lines(coord,correlation_Barnes,lwd=15,col="Blue")
lines(coord,correlation_Soar,lwd=12,col="Sienna")
lines(coord,correlation_Toar,lwd=10,col="Tan")
lines(coord,correlation_Powerlaw,lwd=10,col="Plum2")
lines(coord,correlation_Cressman,lwd=10,col="Chartreuse4")
legend(x="topright",col=c("Blue","Sienna","Tan","Plum2","Chartreuse4"),lwd=15,legend=c("Barnes","SOAR","TOAR","Powerlaw","Cressman"),cex=2)
dev.off()

Setup for OI examples

To prepare for the following examples, run the code below to set up the necessary variables. This code generates a grid with Easting and Northing coordinates, where elevation increases with the Easting coordinate. The grid is split into land and sea: points where $x>y$ are classified as land, and all others as sea. Additionally, 30 observations are generated at random locations on the grid. The background values at grid points are set to 0, and the observations are all set to 1. In this setup, the analysis corresponds to the Integral Data Influence (IDI, Uboldi et al. 2008).

dyn.load("/home/cristianl/projects/gridpp/build/swig/R/gridpp.so")
source("/home/cristianl/projects/gridpp/build/swig/R/gridpp.R")
cacheMetaData(1)
# set grid parameters
min_x <- 1
max_x <- 10
min_y <- 1
max_y <- 10
min_elev <- 0
max_elev <- 20
by_x <- .1
by_y <- .1
# elevation is assumed to be a function of x coord
m <- (max_elev - min_elev) / (max_x - min_x)
q <- min_elev - m * min_x
# generate the regular 2D grid
grid_x_coord <- seq( min_x, max_x, by=by_x)
grid_y_coord <- seq( min_y, max_y, by=by_y)
ngrid_x <- length(grid_x_coord)
ngrid_y <- length(grid_y_coord)
grid_points <- expand.grid( grid_x_coord, grid_y_coord)
grid_x <- as.numeric(grid_points[,1])
grid_y <- as.numeric(grid_points[,2])
ngrid <- length( grid_x)
grid_z <- m * grid_x + q
grid_laf <- rep(0, ngrid)
grid_laf[which(grid_x>grid_y)] <- 1
pgrid  <- Points( grid_x, grid_y, grid_z, grid_laf, "Cartesian")
# generate observations
nobs <- 30
set.seed(1)
obs_x <- runif( nobs, min=min_x, max=max_y)
set.seed(2)
obs_y <- runif( nobs, min=min_y, max=max_y)
set.seed(3)
obs_z <- m * obs_x + q + runif( nobs, min=-3, max=3)
obs_laf <- rep(0, nobs)
obs_laf[which(obs_x>obs_y)] <- 1
points <- Points( obs_x, obs_y, obs_z, obs_laf, "Cartesian")
# set background to 0 everywhere
background <- rep(0, ngrid)
pbackground <- rep(0, nobs)
# set all observations to 1 
obs <- rep(1, nobs)
# OI parameters
ratios <- rep( 0.1, nobs)
max_points <- 5

OI example: Combining Multiple Structure Functions

In this example, we illustrate how to use the optimal_interpolation() function and define a multiple structure function that allows the user to combine three different structure functions into a single one.

To define a multiple structure function, you can use one of the following approaches. Typically, a Barnes structure function is a good choice for describing correlations based on vertical distances (elevation differences) between points. The linear structure function can be used to adjust correlations between points in different surroundings, such as between inland and island locations. Finally, different structure functions (e.g., Barnes, TOAR, Cressman) can be applied to describe correlations based on the radial or horizontal distance between two points.

# Create multiple structure functions
multiple_structure_Barnes <- MultipleStructure( BarnesStructure(0.2, 0, 0), BarnesStructure(0, 1, 0), LinearStructure(0, 0, 0.5))
multiple_structure_Soar <- MultipleStructure( SoarStructure(0.2, 0, 0), BarnesStructure(0, 1, 0), LinearStructure(0, 0, 0.5))
multiple_structure_Toar <- MultipleStructure( ToarStructure(0.2, 0, 0), BarnesStructure(0, 1, 0), LinearStructure(0, 0, 0.5))
multiple_structure_Powerlaw <- MultipleStructure( PowerlawStructure(0.2, 0, 0), BarnesStructure(0, 1, 0), LinearStructure(0, 0, 0.5))
multiple_structure_Cressman <- MultipleStructure( CressmanStructure(0.2, 0, 0), BarnesStructure(0, 1, 0), LinearStructure(0, 0, 0.5))

Gridpp allows you to perform OI by calling the optimal_interpolation() function.

# analysis, in this case is equal to the Integral Data Influence IDI 
analysis_Barnes <- optimal_interpolation(pgrid, background, points, obs, ratios, pbackground, multiple_structure_Barnes, max_points)
analysis_Soar <- optimal_interpolation(pgrid, background, points, obs, ratios, pbackground, multiple_structure_Soar, max_points)
analysis_Toar <- optimal_interpolation(pgrid, background, points, obs, ratios, pbackground, multiple_structure_Toar, max_points)
analysis_Powerlaw <- optimal_interpolation(pgrid, background, points, obs, ratios, pbackground, multiple_structure_Powerlaw, max_points)
analysis_Cressman <- optimal_interpolation(pgrid, background, points, obs, ratios, pbackground, multiple_structure_Cressman, max_points)

Once defined, you can plot the analysis fields obtained from the interpolation. In the analysis, values will be close to 1 near observation locations (indicated by warm colors in the figures) and approach 0 further from the observations (shown by yellow-ish colors, with beige indicating values smaller than 0.1).

The analysis fields effectively show the correlation between grid points near observations and their surrounding grid points. The spatial patterns of different horizontal correlations are clearly marked in the figures. Given the same value of $L$, Cressman provides the smallest neighborhoods influenced by observations, while TOAR or Powerlaw extend the influence over larger areas. The right panels of the figures show how elevation differences affect the correlations. The dashed line in the left panels marks the sharp boundary between land and sea, and you can see how correlations diminish after crossing this boundary.

idi_Barnes <- array(data=analysis_Barnes, dim=c(ngrid_x,ngrid_y))
idi_Soar <- array(data=analysis_Soar, dim=c(ngrid_x,ngrid_y))
idi_Toar <- array(data=analysis_Toar, dim=c(ngrid_x,ngrid_y))
idi_Powerlaw <- array(data=analysis_Powerlaw, dim=c(ngrid_x,ngrid_y))
idi_Cressman <- array(data=analysis_Cressman, dim=c(ngrid_x,ngrid_y))
# Figures
breaks <- c( -1, seq( 0,1, by=0.1), 2)
col <- c( "beige", rev( heat.colors(length(breaks)-2)))
for (i in 1:5) {
  if (i == 1) { str <- "Barnes"; idi <- idi_Barnes; analysis <- analysis_Barnes } else 
  if (i == 2) { str <- "SOAR"; idi <- idi_Soar; analysis <- analysis_Soar } else 
  if (i == 3) { str <- "TOAR"; idi <- idi_Toar; analysis <- analysis_Toar } else 
  if (i == 4) { str <- "Powerlaw"; idi <- idi_Powerlaw; analysis <- analysis_Powerlaw } else 
  if (i == 5) { str <- "Cressman"; idi <- idi_Cressman; analysis <- analysis_Cressman }
  png( file=paste0("fa_",str,".png"), width=800, height=800)
  par( mar = c( 5, 5, 1, 1))
  image( z=idi, x=grid_x_coord, y=grid_y_coord, breaks=breaks, col=col,
         xlab="Easting Coordinate", ylab="Northing Coordinate", cex.lab=2, cex.axis=2)
  points( obs_x, obs_y, pch=21, bg="cornflowerblue", cex=2)
  lines(-1000:1000,-1000:1000,lty=2)
  text(x=8.8,y=8,labels="Land",cex=2.5)
  text(x=8,y=8.7,labels="Sea",cex=2.5)
  dev.off()
  #
  png( file=paste0("fb_",str,".png"), width=800, height=800)
  par( mar = c( 5, 5, 1, 1))
  plot( grid_x, grid_z,
        xlab="Easting Coordinate", ylab="Elevation", cex.lab=2, cex.axis=2)
  for (i in 1:length(col)) {
    if (length(ix <- which(analysis >= breaks[i] & analysis < breaks[i+1]))>0) {
      points( grid_x[ix], grid_z[ix], col=col[i], bg=col[i], pch=21, cex=4)
    }
  }
  text(x=3,y=18,labels=str,cex=4.5)
  points(obs_x,obs_z,pch=21, bg="cornflowerblue", cex=3)
  dev.off()
  system( paste0("convert +append fa_",str,".png fb_",str,".png fig_",str,".png"))
  system( paste0("rm fa_",str,".png fb_",str,".png"))
}

Required setup for examples

To get ready for the examples in the next sections, run the code below to set up necessary variables (you will also need the test datasets). This retrieves air temperature and precipitation from the observation, analysis, and forecast files as well as metadata about the grids.

dyn.load("R/gridpp.so")
source("R/gridpp.R")
cacheMetaData(1)

require(ncdf4)

array2list = function(ar) {
    q = lapply(seq(dim(ar)[2]), function(x) ar[ , x])
    return(q)
}

nc = nc_open("analysis.nc")
ilats = ncvar_get(nc, 'latitude')
ilons = ncvar_get(nc, 'longitude')
ielevs = ncvar_get(nc, 'surface_geopotential')
igrid = Grid(array2list(ilats), array2list(ilons), array2list(ielevs))
temp_analysis = ncvar_get(nc, 'air_temperature_2m')
nc_close(nc)

nc = nc_open("output.nc")
olats = ncvar_get(nc, 'latitude')
olons = ncvar_get(nc, 'longitude')
oelevs = ncvar_get(nc, 'altitude')
ogrid = Grid(array2list(olats), array2list(olons), array2list(oelevs))
nc_close(nc)

nc = nc_open("obs.nc")
plats = ncvar_get(nc, 'latitude')
plons = ncvar_get(nc, 'longitude')
pelevs = ncvar_get(nc, 'altitude')
points = Points(plats, plons, pelevs)
temp_obs = ncvar_get(nc, 'air_temperature_2m')
precip_obs = ncvar_get(nc, 'precipitation_amount')
nc_close(nc)

Optimal interpolation example

We have verified that the following example works with swig version 4.2.1 while it does not work with version 4.0.2 (our best guess: there is a problem with the way cpp shared pointers are handled in 4.0.2, which has been fixed after version 4.1.1).

structure = BarnesStructure(10000, 0)
pbackground = nearest(igrid, points, array2list(temp_analysis[,,1]))
ratios = rep(0.1 , length(pbackground))
max_points = 5
ovalues = optimal_interpolation(igrid, array2list(temp_analysis[,,1]), points, temp_obs, ratios, pbackground, structure, max_points)