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.

Defined grids

RingGrids defines and exports the following grids, see also Implemented grids for a more mathematical description.

Full grids

using SpeedyWeather.RingGrids
subtypes(RingGrids.AbstractFullGrid)

4-element Vector{Any}:
 FullClenshawGrid
 FullGaussianGrid
 FullHEALPixGrid
 FullOctaHEALPixGrid

and reduced grids

subtypes(RingGrids.AbstractReducedGrid)

5-element Vector{Any}:
 HEALPixGrid
 OctaHEALPixGrid
 OctahedralClenshawGrid
 OctahedralGaussianGrid
 OctaminimalGaussianGrid

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.

Grid versus Field

With "grid" we mean the discretization of space. Also called tesselation given that we are tiling a space with polygons, we subdivide the sphere into grid cells, with vertices, faces and centres. In that sense, a grid does not contain any data it purely describes the location of grid cells. Grids in RingGrids are identified by their name, e.g. FullGaussianGrid, and a resolution parameter where we use nlat_half (the number of latitudes on one hemisphere, Equator included) for all grids. This is because some grids have an even number some an odd number number of latitudes so not all nlat are valid, but all nlat_half are. While an instance of a grid stores some precomputed arrays to facilitate faster indexing it does not store coordinates and similar grid information, these can be recomputed on the fly whenever needed given a grid. All grids are considered to be two-dimensional, so they do not contain information about the vertical or time, for example. The horizontal grid points are unravelled into a vector, starting at 0˚E at the north pole, going first east, then ring by ring to the south pole.

Data on a grid is called a Field, many variables, like temperature are a field. Surface temperature would be a 2D field (though represented as a vector), temperature on several vertical layers of the atmosphere would be 3D (data represented as a matrix, horizontal x vertical), including time would make it 4D. Several fields can share the same grid. Given that the grid is always two-dimensional, a 2D and 3D field can also share the same grid, leaving the 3rd dimension not further specified for flexibility.

Creating a grid

All grids are specified by name and the resolution parameter nlat_half::Integer (number of latitude rings on one hemisphere, Equator included). An instance of a grid is simply created by

grid = FullGaussianGrid(24)

48-ring FullGaussianGrid
├ nlat_half=24 (4608 points, ~2.99˚, full)
└ architecture: SpeedyWeather.Architectures.CPU{KernelAbstractions.CPU}

As a second argument architecture can be provided which helps to share information on the computing architecture (CPU/GPU) but this will not be further explained here.

Accessing coordinates

With a grid you can get the coordinates through get_lat, get_latd, get_colat, get_lond but note that the latter is only defined for full grids as the reduced grids do not share the same longitudes across rings. To obtain the longitudes and latitudes for all grid points use get_londlatds, get_lonlats, get_loncolats, e.g.

grid = OctaminimalGaussianGrid(2)   # a tiny grid with 4 latitudes only
get_latd(grid)

4-element Vector{Float64}:
  59.44440828916677
  19.87571914744092
 -19.875719147440904
 -59.44440828916677

Creating a Field

Creating Field, that means data on a grid can be done in many ways, for example using zeros, ones,rand, orrandn`

grid = HEALPixGrid(2)               # smallest HEALPix
field = zeros(grid)                 # 2D field
field = rand(grid, 10)              # 3D field with 10 vertical layers or time steps
field = randn(Float32, grid, 2, 2)  # 4D using Float32 as element type

12×2×2, 3-ring HEALPixField{Float32, 3} on Array on SpeedyWeather.Architectures.CPU{KernelAbstractions.CPU}:
[:, :, 1] =
 -0.524171   0.104508
  0.917835  -0.400061
  0.396045   0.597318
  0.418773  -1.59478
  1.45856    0.214578
 -0.553765  -0.0401893
 -0.431075  -0.266735
 -1.09957   -1.95032
 -1.20999    0.121229
 -0.592127  -0.544872
 -0.377673  -0.719442
 -2.02139   -0.394944

[:, :, 2] =
 -0.850132  -0.022275
  0.903572   1.66832
 -1.75085   -0.485193
 -1.59219    0.102127
 -0.76742    0.884069
 -1.00373    0.465761
  1.24238   -0.822543
  1.08514   -1.90568
  0.338575   1.12088
  0.715901  -0.723818
 -2.41422    0.724279
  0.254568   0.0360136

Note that in this case all fields share the same grid. Or the grid can be created on the fly when the grid type is specified, followed by nlat_half. Every ?Grid also has a corresponding ?Field type, e.g.

field = zeros(OctaminimalGaussianGrid, 2)   # nlat_half=2
field = HEALPixField(undef, 2)              # using undef initializor
field = HEALPixField{Float16}(undef, 2, 3)  # using Float16 as eltype

12×3, 3-ring HEALPixField{Float16, 2} on Array on SpeedyWeather.Architectures.CPU{KernelAbstractions.CPU}:
 6.0e-8  0.0  0.0
 0.0     0.0  0.0
 0.0     0.0  0.0
 0.0     0.0  0.0
 6.6e-7  0.0  0.0
 0.0     0.0  0.0
 0.0     0.0  0.0
 0.0     0.0  0.0
 0.0     0.0  0.000839
 0.0     0.0  0.0
 0.0     0.0  4.005e-5
 0.0     0.0  0.0

Creating a Field from data

A field has field.data (some AbstractArray{T, N}) and field.grid (some AbstractGrid as described above). The first dimension of data describes the horizontal as the grid points on every grid (full and reduced) 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

data = randn(Float32, 8, 4)
field = FullGaussianGrid(data, input_as=Matrix)
field = FullGaussianField(data, input_as=Matrix)    # equivalent

32-element, 4-ring FullGaussianField{Float32, 1} on Array on SpeedyWeather.Architectures.CPU{KernelAbstractions.CPU}:
 -0.65970683
  1.5565485
 -1.3304058
 -0.2069787
  0.8236686
  0.37665534
 -0.49068186
 -1.4700775
  0.40112665
  0.12839806
  ⋮
 -0.025904706
  0.51760197
  0.031388912
 -0.9150743
 -0.939553
 -1.7625887
  0.30175877
 -0.1052116
 -1.8279105

Both return a field (there is no data in the grid itself) and create an instance of FullGaussianGrid on the fly. The input_as=Matrix is required to denote that the horizontal dimension is not unravelled into a vector, triggering a reshape internally, input_as=Vector is the default. 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.

If you have the data and know which grid it comes one you can also create a Field by providing both

data = randn(Float32, 8, 4)         # data of some shape
grid = OctaminimalGaussianGrid(1)   # you need to know the nlat_half (here 1) of that grid!
field = Field(data, grid)

8×4, 2-ring OctaminimalGaussianField{Float32, 2} on Array on SpeedyWeather.Architectures.CPU{KernelAbstractions.CPU}:
 -1.48815   -0.0706487  -1.86668    -0.595802
 -1.40278    1.06431    -0.620261   -1.4903
 -0.779348   0.414487   -0.0192575  -2.01905
 -0.467307  -1.34075    -0.2235      0.145223
  0.808506  -0.292958    0.325626    1.47968
  0.579691   0.30291    -1.69801     1.33374
  1.00979   -0.474029    1.47142    -0.582973
  0.643891  -0.479903   -0.490624   -0.680521

But you can also automatically let nlat_half be calculated from the shape of the data. Note that you have to provide the name of the field though, FullGaussianField here to create a Field on a FullGaussianGrid.

field = FullGaussianField(data, input_as=Matrix)

32-element, 4-ring FullGaussianField{Float32, 1} on Array on SpeedyWeather.Architectures.CPU{KernelAbstractions.CPU}:
 -1.488153
 -1.4027804
 -0.7793485
 -0.46730664
  0.8085065
  0.579691
  1.0097895
  0.6438906
 -0.07064872
  1.064311
  ⋮
 -0.4906239
 -0.59580165
 -1.4902976
 -2.019049
  0.14522256
  1.479685
  1.3337358
 -0.5829725
 -0.6805206

To return to the original data array you can reshape the data of a full grid (which is representable as a matrix) as follows

data == Matrix(FullGaussianField(data, input_as=Matrix))

true

which in general is Array(field, as=Vector) for no reshaping (equivalent to field.data including possible conversion to Array) and Array(field, as=Matrix) with reshaping (full grids only).

Visualising Fields

As only the full fields can be reshaped into a matrix, the underlying data structure of any AbstractField is a vector for consistency. As shown in the examples above, one can therefore inspect the data as if it was a vector. But a field has through field.grid all the geometric information available to plot it on a map, SpeedyWeather also implements extensions for Makie's and UnicodePlots's' heatmap, also see Visualisation via Makie and Visualisation via UnicodePlots.

using CairoMakie    # triggers loading of Makie extension, or do using UnicodePlots instead!
grid = OctahedralGaussianGrid(24)
field = randn(grid)
heatmap(field)

Reduced fields are automatically interpolated to the corresponding full fields so that they can be visualised as a matrix.

Indexing Fields

All fields 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 two-dimensional grid covering the sphere. 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_londlatds for latitudes, longitudes per grid point not per ring). Now we could calculate Coriolis and add it on the grid as follows

grid = OctahedralClenshawGrid(5)    # define a grid
field = randn(grid)                 # random data on that grid
latd = get_latd(grid)               # get a vector of latitudes [˚N] per ring, north to south

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

for (j, ring) in enumerate(eachring(field))     # loop over every ring j
    f = coriolis[j]
    for ij in ring                              # loop over every longitude point on that ring
        field[ij] += f
    end
end

eachring(field) accesses field.grid.rings a precomputed 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 horizontal grid points with eachgridpoint and over every other dimensions with eachlayer, e.g. for 2D fields you can do

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

or use eachindex instead. For 3D fields eachindex loops over all elements, including the 3rd dimension but eachgridpoint would only loop over the horizontal. To loop with an index k over all additional dimensions (vertical, time, ...) do

field = zeros(grid, 2, 3)           # 4D, 2D defined by grid x 2 x 3
for k in eachlayer(field)           # loop over 2 x 3
    for ij in eachgridpoint(field)  # loop over 2D grid points
        field[ij, k]
    end
end

Interpolation between grids

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 one grid to another or to request an interpolated value on some coordinates. Using a OctahedralGaussianField (data on an OctahedralGaussianGrid) from above we can use the interpolate function to get it onto a FullGaussianGrid (or any other grid for purpose)

grid = OctahedralGaussianGrid(4)
field = randn(Float32, grid)
interpolate(FullGaussianGrid, field)

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 any other grid as an instance not a type (or provide nlat_half as the second argument), e.g.

output_grid = HEALPixGrid(6)
interpolate(output_grid, field)
interpolate(HEALPixGrid, 6, field)  # equivalent

To change the number format during interpolation you can preallocate the output field

output_field = Field(Float16, output_grid)
interpolate!(output_field, field)

Which will convert from Float64 to Float16 on the fly too. Note that we use interpolate! here not interpolate without ! as the former writes into output_field and therefore changes it in place.

Interpolation on coordinates

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

interpolate(30.0, 10.0, field)

-0.342801f0

we interpolated the data from field 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], field)

3-element Vector{Float32}:
 -0.342801
 -0.4338002
 -0.23844747

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

Performance for RingGrid interpolation

Every time an interpolation like interpolate(30.0, 10.0, field) 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 array
input field
interpolator

The output array 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. Such a destination can be a field but it does not have to be (see Interpolation on coordinates). The input field 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 of the field 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 field 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 field. So in the case you want to interpolate from field A to field B many times, you can just reuse the same interpolator. If you want to change the coordinates of the output but its total number of points remain constant 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 and want to interpolate between them

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

AnvilInterpolator{Float64} for 7-ring HEALPixGrid
├ nlat_half=4 (48 points, ~29.3˚, reduced)
└ architecture: SpeedyWeather.Architectures.CPU{KernelAbstractions.CPU}
└ 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 (it has only seen grids, no fields). We can now do

field_in = rand(grid_in)
field_out = zeros(grid_out)
interpolate!(field_out, field_in, interp)

which is identical to interpolate(field_out, field_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)

field_out = zeros(Float16, grid_out)
interpolate!(field_out, field_in, interp)

and we have converted data from a HEALPixField{Float64} (Float64 is always default if not specified) to a FullClenshawField{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

field_in = randn(OctahedralGaussianField{Float16}, 24)
field_out = zeros(FullClenshawField{Float16}, 24)
interp = RingGrids.interpolator(field_out, field_in, NF=Float32)
interpolate!(field_out, field_in, interp)

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
grid = HEALPixGrid(24)
interp = AnvilInterpolator(grid, npoints, NF=Float32)

AnvilInterpolator{Float32} for 47-ring HEALPixGrid
├ nlat_half=24 (1728 points, ~4.9˚, reduced)
└ architecture: SpeedyWeather.Architectures.CPU{KernelAbstractions.CPU}
└ onto: 10 points

with the first argument being the input grid and then the number of points to interpolate onto. The number format used for the interpolator is provided as NF. 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.

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

now we can update the locator inside our interpolator as follows

RingGrids.update_locator!(interp, londs, latds)

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

output_vec = zeros(10)
field_in = rand(grid)

we can use the interpolator as follows

interpolate!(output_vec, field_in, interp)

10-element Vector{Float64}:
 0.34364689471798743
 0.3312177153044105
 0.8082961890012526
 0.2932432348798577
 0.6400178433709259
 0.6485216040115607
 0.48586444010848395
 0.7810032573522634
 0.6626493973018864
 0.527938969252993

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

interpolate(londs, latds, field_in)

10-element Vector{Float64}:
 0.34364689380787067
 0.3312177157586497
 0.8082961951099354
 0.29324323561665655
 0.6400178434357897
 0.6485216065230774
 0.48586443603138346
 0.7810032532585667
 0.6626494032132237
 0.5279389693024032

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

Create a new grid of type Grid with resolution parameter nlat_half. architecture is the device type (CPU/GPU). Precomputes the ring indices rings.