Gridpp in R

[Under development]

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) 

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)