RingGrids

RingGrids is a submodule that has been developed for SpeedyWeather.jl which is technically independent (SpeedyWeather.jl however imports it and so does SpeedyTransforms) and can also be used without running simulations. It is just not put into its own respective repository.

RingGrids defines several iso-latitude grids, which are mathematically described in the section on Grids. In brief, they include the regular latitude-longitude grids (here called FullClenshawGrid) as well as grids which latitudes are shifted to the Gaussian latitudes and reduced grids, meaning that they have a decreasing number of longitudinal points towards the poles to be more equal-area than full grids.

RingGrids defines and exports the following grids:

full grids: FullClenshawGrid, FullGaussianGrid, FullHEALPix, and FullOctaHEALPix
reduced grids: OctahedralGaussianGrid, OctahedralClenshawGrid, OctaHEALPixGrid and HEALPixGrid

The following explanation of how to use these can be mostly applied to any of them, however, there are certain functions that are not defined, e.g. the full grids can be trivially converted to a Matrix (i.e. they are rectangular grids) but not the OctahedralGaussianGrid.

What is a ring?

We use the term ring, short for iso-latitude ring, to refer to a sequence of grid points that all share the same latitude. A latitude-longitude grid is a ring grid, as it organises its grid-points into rings. However, other grids, like the cubed-sphere are not based on iso-latitude rings. SpeedyWeather.jl only works with ring grids because its a requirement for the Spherical Harmonic Transform to be efficient. See Grids.

Creating data on a RingGrid

Every grid in RingGrids has a grid.data field, which is a vector containing the data on the grid. The grid points are unravelled west to east then north to south, meaning that it starts at 90˚N and 0˚E then walks eastward for 360˚ before jumping on the next latitude ring further south, this way circling around the sphere till reaching the south pole. This may also be called ring order.

Data in a Matrix which follows this ring order can be put on a FullGaussianGrid like so

using SpeedyWeather.RingGrids
map = randn(Float32, 8, 4)

8×4 Matrix{Float32}:
  0.810838   0.906474  -1.7733     0.272695
 -0.180423  -0.122792  -0.793858  -0.794579
 -0.585897   0.472033  -1.37577   -0.223613
  0.391809   0.587566  -0.306473   2.8628
  1.03392   -0.611295  -1.74545    0.37196
  0.836483   0.539314  -0.162308  -0.273887
  0.663684   1.40548    0.798213  -0.22473
 -0.563424   1.49567   -0.258994  -1.55334

grid = FullGaussianGrid(map, input_as=Matrix)

32-element, 4-ring FullGaussianGrid{Float32}:
  0.81083834
 -0.18042298
 -0.58589685
  0.39180875
  1.0339204
  0.8364826
  0.6636845
 -0.56342363
  0.9064742
 -0.122791916
  ⋮
 -0.25899416
  0.27269548
 -0.7945794
 -0.22361264
  2.862804
  0.3719597
 -0.27388656
 -0.2247304
 -1.5533442

Note that input_as=Matrix is necessary as, RingGrids have a flattened horizontal dimension into a vector. To distinguish the 2nd horizontal dimension from a possible vertical dimension the keyword argument here is required.

A full Gaussian grid has always $2N$ x $N$ grid points, but a FullClenshawGrid has $2N$ x $N-1$, if those dimensions don't match, the creation will throw an error. To reobtain the data from a grid, you can access its data field which returns a normal Vector

grid.data

32-element Vector{Float32}:
  0.81083834
 -0.18042298
 -0.58589685
  0.39180875
  1.0339204
  0.8364826
  0.6636845
 -0.56342363
  0.9064742
 -0.122791916
  ⋮
 -0.25899416
  0.27269548
 -0.7945794
 -0.22361264
  2.862804
  0.3719597
 -0.27388656
 -0.2247304
 -1.5533442

Which can be reshaped to reobtain map from above. Alternatively you can Matrix(grid) to do this in one step

map == Matrix(FullGaussianGrid(map))

false

You can also use zeros, ones, rand, randn to create a grid, whereby nlat_half, i.e. the number of latitude rings on one hemisphere, Equator included, is used as a resolution parameter and here as a second argument.

nlat_half = 4
grid = randn(OctahedralGaussianGrid{Float16}, nlat_half)

208-element, 8-ring OctahedralGaussianGrid{Float16}:
 -0.1133
  0.1583
  0.7266
  0.529
  0.4053
 -1.796
 -0.4343
  0.9526
 -0.1082
 -0.853
  ⋮
 -0.6943
  1.536
 -0.4985
  0.3433
  0.534
  1.05
 -0.2954
 -0.2107
 -0.6245

and any element type T can be used for OctahedralGaussianGrid{T} and similar for other grid types.

Visualising RingGrid data

As only the full grids can be reshaped into a matrix, the underlying data structure of any AbstractGrid is a vector. As shown in the examples above, one can therefore inspect the data as if it was a vector. But as that data has, through its <:AbstractGrid type, all the geometric information available to plot it on a map, RingGrids also exports plot function, based on UnicodePlots' heatmap.

nlat_half = 24
grid = randn(OctahedralGaussianGrid, nlat_half)
RingGrids.plot(grid)

                   48-ring OctahedralGaussianGrid{Float64}                
       ┌────────────────────────────────────────────────────────────┐  3  
    90 │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ ┌──┐
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
    ˚N │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
   -90 │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ └──┘
       └────────────────────────────────────────────────────────────┘ -4  
        0                           ˚E                           360

(Note that to skip the RingGrids. in the last line you can do import SpeedyWeather.RingGrids: plot, import SpeedyWeather: plot or simply using SpeedyWeather. It's just that LowerTriangularMatrices also defines plot which otherwise causes naming conflicts.)

Indexing RingGrids

All RingGrids have a single index ij which follows the ring order. While this is obviously not super exciting here are some examples how to make better use of the information that the data sits on a grid.

We obtain the latitudes of the rings of a grid by calling get_latd (get_lond is only defined for full grids, or use get_latdlonds for latitudes, longitudes per grid point not per ring)

grid = randn(OctahedralClenshawGrid, 5)
latd = get_latd(grid)

9-element Vector{Float64}:
  72.0
  54.0
  36.0
  18.0
   0.0
 -18.0
 -36.0
 -54.0
 -72.0

Now we could calculate Coriolis and add it on the grid as follows

rotation = 7.29e-5                  # angular frequency of Earth's rotation [rad/s]
coriolis = 2rotation*sind.(latd)    # vector of coriolis parameters per latitude ring

rings = eachring(grid)
for (j, ring) in enumerate(rings)
    f = coriolis[j]
    for ij in ring
        grid[ij] += f
    end
end

eachring creates a vector of UnitRange indices, such that we can loop over the ring index j (j=1 being closest to the North pole) pull the coriolis parameter at that latitude and then loop over all in-ring indices i (changing longitudes) to do something on the grid. Something similar can be done to scale/unscale with the cosine of latitude for example. We can always loop over all grid-points like so

for ij in eachgridpoint(grid)
    grid[ij]
end

or use eachindex instead.

Interpolation on RingGrids

In most cases we will want to use RingGrids so that our data directly comes with the geometric information of where the grid-point is one the sphere. We have seen how to use get_latd, get_lond, ... for that above. This information generally can also be used to interpolate our data from grid to another or to request an interpolated value on some coordinates. Using our data on grid which is an OctahedralGaussianGrid from above we can use the interpolate function to get it onto a FullGaussianGrid (or any other grid for purpose)

grid = randn(OctahedralGaussianGrid{Float32}, 4)

208-element, 8-ring OctahedralGaussianGrid{Float32}:
  1.0485916
 -0.20203464
 -1.2415137
  0.95038855
  0.10498307
  0.05388951
 -1.1637807
 -0.83020276
 -1.0859487
 -0.8832507
  ⋮
  0.061498024
 -1.9674281
 -0.73826927
 -1.2760131
  0.064602114
  0.8781293
  0.96678364
  0.4221278
 -0.6415451

interpolate(FullGaussianGrid, grid)

128-element, 8-ring FullGaussianGrid{Float64}:
  1.0485916137695312
 -0.4619044214487076
 -0.14556264877319336
  0.31633445620536804
  0.05388950929045677
 -1.0803862065076828
 -0.9580757319927216
 -0.933925211429596
 -0.9254169464111328
 -0.13810716569423676
  ⋮
  0.9270212054252625
 -0.8849471211433411
 -0.4457335099577904
 -1.3528486490249634
 -1.1415771692991248
  0.06460211426019669
  0.9002928882837296
  0.6944557428359985
 -0.3756268620491028

By default this will linearly interpolate (it's an Anvil interpolator, see below) onto a grid with the same nlat_half, but we can also coarse-grain or fine-grain by specifying nlat_half directly as 2nd argument

interpolate(FullGaussianGrid, 6, grid)

288-element, 12-ring FullGaussianGrid{Float64}:
  0.6521090668277228
 -0.05679980851672897
 -0.6699589011841325
 -0.1601687033800131
  0.393625277453864
  0.0044630330805663365
 -0.024499000313265296
 -0.7147269910631263
 -0.7015034994691319
 -0.7128496560877943
  ⋮
 -0.8752726973535553
 -0.5791551440837058
 -0.6710245660234825
  0.08889397817203225
  0.5500362650989591
  0.6824672203521456
  0.5173278221273485
  0.09091289118318141
 -0.3786162483306045

So we got from an 8-ring OctahedralGaussianGrid{Float16} to a 12-ring FullGaussianGrid{Float64}, so it did a conversion from Float16 to Float64 on the fly too, because the default precision is Float64 unless specified. interpolate(FullGaussianGrid{Float16}, 6, grid) would have interpolated onto a grid with element type Float16.

One can also interpolate onto a given coordinate ˚N, ˚E like so

interpolate(30.0, 10.0, grid)

-0.42760348f0

we interpolated the data from grid onto 30˚N, 10˚E. To do this simultaneously for many coordinates they can be packed into a vector too

interpolate([30.0, 40.0, 50.0], [10.0, 10.0, 10.0], grid)

3-element Vector{Float32}:
 -0.42760348
 -0.59949327
 -0.7839959

which returns the data on grid at 30˚N, 40˚N, 50˚N, and 10˚E respectively. Note how the interpolation here retains the element type of grid.

Performance for RingGrid interpolation

Every time an interpolation like interpolate(30.0, 10.0, grid) is called, several things happen, which are important to understand to know how to get the fastest interpolation out of this module in a given situation. Under the hood an interpolation takes three arguments

output vector
input grid
interpolator

The output vector is just an array into which the interpolated data is written, providing this prevents unnecessary allocation of memory in case the destination array of the interpolation already exists. The input grid contains the data which is subject to interpolation, it must come on a ring grid, however, its coordinate information is actually already in the interpolator. The interpolator knows about the geometry of the grid the data is coming on and the coordinates it is supposed to interpolate onto. It has therefore precalculated the indices that are needed to access the right data on the input grid and the weights it needs to apply in the actual interpolation operation. The only thing it does not know is the actual data values of that grid. So in the case you want to interpolate from grid A to grid B many times, you can just reuse the same interpolator. If you want to change the coordinates of the output grid but its total number of points remain constants then you can update the locator inside the interpolator and only else you will need to create a new interpolator. Let's look at this in practice. Say we have two grids an want to interpolate between them

grid_in = rand(HEALPixGrid, 4)
grid_out = zeros(FullClenshawGrid, 6)
interp = RingGrids.interpolator(grid_out, grid_in)

AnvilInterpolator{Float64, HEALPixGrid}
├ from: 7-ring HEALPixGrid, 48 grid points
└ onto: 264 points

Now we have created an interpolator interp which knows about the geometry where to interpolate from and the coordinates there to interpolate to. It is also initialized, meaning it has precomputed the indices to of grid_in that are supposed to be used. It just does not know about the data of grid_in (and neither of grid_out which will be overwritten anyway). We can now do

interpolate!(grid_out, grid_in, interp)
grid_out

264-element, 11-ring FullClenshawGrid{Float64}:
 0.7583144976728629
 0.7313216572248218
 0.7043288167767806
 0.6773359763287395
 0.6563328140159135
 0.6353296517030875
 0.6143264893902615
 0.5933233270774354
 0.5723201647646095
 0.5513170024517834
 ⋮
 0.42806040036510895
 0.38842002361625394
 0.3487796468673989
 0.30913927011854375
 0.2694988933696887
 0.22985851662083379
 0.19021813987197855
 0.2647823090394976
 0.3393464782070165

which is identical to interpolate(grid_out, grid_in) but you can reuse interp for other data. The interpolation can also handle various element types (the interpolator interp does not have to be updated for this either)

grid_out = zeros(FullClenshawGrid{Float16}, 6);
interpolate!(grid_out, grid_in, interp)
grid_out

264-element, 11-ring FullClenshawGrid{Float16}:
 0.7583
 0.7314
 0.704
 0.6772
 0.6562
 0.6353
 0.6143
 0.5933
 0.5723
 0.5513
 ⋮
 0.428
 0.3884
 0.3489
 0.309
 0.2695
 0.2299
 0.1902
 0.265
 0.3394

and we have converted data from a HEALPixGrid{Float64} (Float64 is always default if not specified) to a FullClenshawGrid{Float16} including the type conversion Float64-Float16 on the fly. Technically there are three data types and their combinations possible: The input data will come with a type, the output array has an element type and the interpolator has precomputed weights with a given type. Say we want to go from Float16 data on an OctahedralGaussianGrid to Float16 on a FullClenshawGrid but using Float32 precision for the interpolation itself, we would do this by

grid_in = randn(OctahedralGaussianGrid{Float16}, 24)
grid_out = zeros(FullClenshawGrid{Float16}, 24)
interp = RingGrids.interpolator(Float32, grid_out, grid_in)
interpolate!(grid_out, grid_in, interp)
grid_out

4512-element, 47-ring FullClenshawGrid{Float16}:
 -0.3645
 -0.61
 -0.8555
 -1.101
 -1.346
 -1.38
 -1.282
 -1.185
 -1.087
 -0.9463
  ⋮
  0.6226
  0.8354
  0.9062
  0.977
  1.048
  0.976
  0.926
  0.876
  0.8257

As a last example we want to illustrate a situation where we would always want to interpolate onto 10 coordinates, but their locations may change. In order to avoid recreating an interpolator object we would do (AnvilInterpolator is described in Anvil interpolator)

npoints = 10    # number of coordinates to interpolate onto
interp = AnvilInterpolator(Float32, HEALPixGrid, 24, npoints)

AnvilInterpolator{Float32, HEALPixGrid}
├ from: 47-ring HEALPixGrid, 1728 grid points
└ onto: 10 points

with the first argument being the number format used during interpolation, then the input grid type, its resolution in terms of nlat_half and then the number of points to interpolate onto. However, interp is not yet initialized as it does not know about the destination coordinates yet. Let's define them, but note that we already decided there's only 10 of them above.

latds = collect(0.0:5.0:45.0)
londs = collect(-10.0:2.0:8.0)

now we can update the locator inside our interpolator as follows

RingGrids.update_locator!(interp, latds, londs)

With data matching the input from above, a nlat_half=24 HEALPixGrid, and allocate 10-element output vector

output_vec = zeros(10)
grid_input = rand(HEALPixGrid, 24)

we can use the interpolator as follows

interpolate!(output_vec, grid_input, interp)

10-element Vector{Float64}:
 0.5817838849073665
 0.46326375754275684
 0.6071189329363699
 0.40447637533898256
 0.540618570781595
 0.4233420259692797
 0.7233977804893904
 0.453832092539283
 0.7982269864525035
 0.6040032217908099

which is the approximately the same as doing it directly without creating an interpolator first and updating its locator

interpolate(latds, londs, grid_input)

10-element Vector{Float64}:
 0.5817838867202452
 0.46326375746429205
 0.6071189349086645
 0.40447639308795535
 0.5406185704545488
 0.4233420243447585
 0.7233977887502223
 0.45383209586060913
 0.7982269874875757
 0.6040032224998209

but allows for a reuse of the interpolator. Note that the two output arrays are not exactly identical because we manually set our interpolator interp to use Float32 for the interpolation whereas the default is Float64.

Anvil interpolator

Currently the only interpolator implemented is a 4-point bilinear interpolator, which schematically works as follows. Anvil interpolation is the bilinear average of a, b, c, d which are values at grid points in an anvil-shaped configuration at location x, which is denoted by Δab, Δcd, Δy, the fraction of distances between a-b, c-d, and ab-cd, respectively. Note that a, c and b, d do not necessarily share the same longitude/x-coordinate.

        0..............1    # fraction of distance Δab between a, b
        |<  Δab   >|

0^      a -------- o - b    # anvil-shaped average of a, b, c, d at location x
.Δy                |
.                  |
.v                 x 
.                  |
1         c ------ o ---- d

          |<  Δcd >|
          0...............1 # fraction of distance Δcd between c, d

^ fraction of distance Δy between a-b and c-d.

This interpolation is chosen as by definition of the ring grids, a and b share the same latitude, so do c and d, but the longitudes can be different for all four, a, b, c, d.

Function index

Core.Type — Method

() Initialize an instance of the grid from an Array. For keyword argument input_as=Vector (default) the leading dimension is interpreted as a flat vector of all horizontal entries in one layer. For input_as==Matrx the first two leading dimensions are interpreted as longitute and latitude. This is only possible for full grids that are a subtype of AbstractFullGridArray.