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
.
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.
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.879732 0.5415 1.12268 -1.50739
+ 0.285329 -0.624654 -0.761027 0.207444
+ 1.16048 0.836968 0.146837 -0.545373
+ 0.95938 -1.50522 1.19983 1.2712
+ 0.484405 -0.551759 0.714948 -1.1101
+ -0.810506 0.0394689 -0.299618 0.7849
+ 1.10882 -0.495427 -2.0568 -0.645783
+ -0.0262196 0.431795 -1.45129 0.544136
grid = FullGaussianGrid(map)
32-element, 4-ring FullGaussianGrid{Float32}:
+ 0.87973154
+ 0.28532878
+ 1.1604829
+ 0.95938015
+ 0.4844051
+ -0.81050587
+ 1.1088202
+ -0.02621964
+ 0.54150033
+ -0.6246537
+ ⋮
+ -1.4512941
+ -1.5073909
+ 0.20744395
+ -0.54537296
+ 1.271205
+ -1.1100991
+ 0.78489953
+ -0.6457826
+ 0.54413635
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.87973154
+ 0.28532878
+ 1.1604829
+ 0.95938015
+ 0.4844051
+ -0.81050587
+ 1.1088202
+ -0.02621964
+ 0.54150033
+ -0.6246537
+ ⋮
+ -1.4512941
+ -1.5073909
+ 0.20744395
+ -0.54537296
+ 1.271205
+ -1.1100991
+ 0.78489953
+ -0.6457826
+ 0.54413635
Which can be reshaped to reobtain map
from above. Alternatively you can Matrix(grid)
to do this in one step
map == Matrix(FullGaussianGrid(map))
true
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.12177
+ -1.136
+ -0.1753
+ 0.01573
+ 2.08
+ 0.1267
+ -0.1635
+ -0.1525
+ -1.063
+ 1.19
+ ⋮
+ -0.7705
+ -0.1923
+ 0.11285
+ 0.26
+ -0.2717
+ -2.05
+ -0.2732
+ -0.495
+ -0.648
and any element type T
can be used for OctahedralGaussianGrid{T}
and similar for other grid types.
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)
+plot(grid)
48-ring OctahedralGaussianGrid{Float64}
+ ┌────────────────────────────────────────────────────────────┐ 3
+ 90 │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ ┌──┐
+ │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
+ │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
+ │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
+ │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
+ │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
+ │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
+ ˚N │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
+ │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
+ │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
+ │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
+ │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
+ │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
+ │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
+ -90 │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ └──┘
+ └────────────────────────────────────────────────────────────┘ -3
+ 0 ˚E 360
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.
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}:
+ -0.5153889
+ 1.2055534
+ -2.0790656
+ -0.06593549
+ -0.2374435
+ -0.6407559
+ 0.5670635
+ -1.055078
+ 0.50161904
+ -0.9276229
+ ⋮
+ 1.2072601
+ 0.10708638
+ 2.1995041
+ 0.21880381
+ 0.547214
+ 0.3342631
+ -1.0475396
+ 1.5327549
+ 0.3171192
interpolate(FullGaussianGrid,grid)
128-element, 8-ring FullGaussianGrid{Float64}:
+ -0.5153889060020447
+ 0.3843986988067627
+ -1.072500467300415
+ -0.1945665031671524
+ -0.6407558917999268
+ 0.1615281105041504
+ -0.2767294645309448
+ -0.5703123807907104
+ -0.34673571586608887
+ 0.14624671638011932
+ ⋮
+ 0.7905355095863342
+ -0.12465041875839233
+ 0.9322167038917542
+ 1.1532952412962914
+ 0.7139788568019902
+ 0.5472139716148376
+ -0.011187583208084106
+ 0.2426077127456665
+ 0.6210281550884247
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.42008688951798034
+ 0.555417409849856
+ -0.7389725999233231
+ -0.7990407240123742
+ -0.15325064895476756
+ -0.2767481319297421
+ -0.5053629873185568
+ 0.1792811788754473
+ -0.41938988787787085
+ -0.2577479312795486
+ ⋮
+ 0.9256958547717504
+ 1.1882408858903037
+ 0.3272752016088992
+ 0.5134322488509792
+ 0.3927225021445712
+ -0.25803173308952543
+ 0.3062354900178891
+ 0.9081788858153896
+ 0.4646344104532309
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.73378944f0
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.73378944
+ 0.6409761
+ 0.343434
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
.
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)
SpeedyWeather.RingGrids.AnvilInterpolator{Float64, HEALPixGrid{Float64}}(SpeedyWeather.RingGrids.GridGeometry{HEALPixGrid{Float64}}(4, 7, 48, [90.0, 66.44353569089876, 41.8103148957786, 19.471220634490685, 0.0, -19.47122063449071, -41.81031489577862, -66.44353569089876, -90.00000000000001], [45.0, 135.0, 225.0, 315.0, 22.5, 67.5, 112.5, 157.5, 202.5, 247.49999999999997 … 112.5, 157.5, 202.5, 247.49999999999997, 292.5, 337.5, 45.0, 135.0, 225.0, 315.0], UnitRange{Int64}[1:4, 5:12, 13:20, 21:28, 29:36, 37:44, 45:48], [4, 8, 8, 8, 8, 8, 4], [45.0, 22.5, 0.0, 22.5, 0.0, 22.5, 45.0]), SpeedyWeather.RingGrids.AnvilLocator{Float64}(264, [0, 0, 0, 0, 0, 0, 0, 0, 0, 0 … 7, 7, 7, 7, 7, 7, 7, 7, 7, 7], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0 … 46, 46, 47, 47, 47, 47, 47, 48, 48, 48], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0 … 47, 47, 48, 48, 48, 48, 48, 45, 45, 45], [4, 4, 4, 1, 1, 1, 1, 1, 1, 2 … -1, -1, -1, -1, -1, -1, -1, -1, -1, -1], [1, 1, 1, 2, 2, 2, 2, 2, 2, 3 … -1, -1, -1, -1, -1, -1, -1, -1, -1, -1], [0.6367678868600251, 0.6367678868600251, 0.6367678868600251, 0.6367678868600251, 0.6367678868600251, 0.6367678868600251, 0.6367678868600251, 0.6367678868600251, 0.6367678868600251, 0.6367678868600251 … 0.3632321131399741, 0.3632321131399741, 0.3632321131399741, 0.3632321131399741, 0.3632321131399741, 0.3632321131399741, 0.3632321131399741, 0.3632321131399741, 0.3632321131399741, 0.3632321131399741], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.8333333333333333, 0.9999999999999998, 0.16666666666666652, 0.33333333333333304, 0.5, 0.6666666666666665, 0.8333333333333326, 0.0, 0.16666666666666652, 0.3333333333333326], [0.5, 0.6666666666666667, 0.8333333333333333, 0.0, 0.16666666666666657, 0.33333333333333315, 0.5, 0.6666666666666666, 0.8333333333333331, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]))
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.4183279674982502
+ 0.3306960264587438
+ 0.24306408541923763
+ 0.15543214437973124
+ 0.15721474249354195
+ 0.15899734060735266
+ 0.1607799387211634
+ 0.1625625368349741
+ 0.1643451349487848
+ 0.1661277330625955
+ ⋮
+ 0.4332268708124246
+ 0.40365038093864175
+ 0.374073891064859
+ 0.34449740119107614
+ 0.3149209113172934
+ 0.2853444214435107
+ 0.25576793156972777
+ 0.2634122069207885
+ 0.2710564822718492
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.4182
+ 0.3308
+ 0.243
+ 0.1554
+ 0.1572
+ 0.159
+ 0.1608
+ 0.1626
+ 0.1643
+ 0.1661
+ ⋮
+ 0.433
+ 0.4036
+ 0.374
+ 0.3445
+ 0.315
+ 0.2854
+ 0.2559
+ 0.2634
+ 0.271
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.4338
+ -0.471
+ -0.508
+ -0.545
+ -0.582
+ -0.6484
+ -0.6157
+ -0.583
+ -0.5503
+ -0.4473
+ ⋮
+ -0.09143
+ -0.2184
+ -0.2346
+ -0.251
+ -0.2673
+ -0.1527
+ -0.1847
+ -0.2166
+ -0.2485
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)
SpeedyWeather.RingGrids.AnvilInterpolator{Float32, HEALPixGrid}(SpeedyWeather.RingGrids.GridGeometry{HEALPixGrid}(24, 47, 1728, [90.0, 86.10076357950555, 82.19700324028634, 78.28414760510762, 74.35752898700072, 70.41233167174659, 66.44353569089876, 62.445854167002665, 58.41366190347208, 54.34091230386124 … -54.340912303861266, -58.41366190347208, -62.445854167002665, -66.44353569089876, -70.41233167174661, -74.35752898700072, -78.28414760510763, -82.19700324028634, -86.10076357950557, -90.00000000000001], [45.0, 135.0, 225.0, 315.0, 22.5, 67.5, 112.5, 157.5, 202.5, 247.49999999999997 … 112.5, 157.5, 202.5, 247.49999999999997, 292.5, 337.5, 45.0, 135.0, 225.0, 315.0], UnitRange{Int64}[1:4, 5:12, 13:24, 25:40, 41:60, 61:84, 85:112, 113:144, 145:180, 181:220 … 1509:1548, 1549:1584, 1585:1616, 1617:1644, 1645:1668, 1669:1688, 1689:1704, 1705:1716, 1717:1724, 1725:1728], [4, 8, 12, 16, 20, 24, 28, 32, 36, 40 … 40, 36, 32, 28, 24, 20, 16, 12, 8, 4], [45.0, 22.5, 14.999999999999998, 11.25, 9.0, 7.499999999999999, 6.428571428571429, 5.625, 5.0, 4.5 … 4.5, 5.0, 5.625, 6.428571428571429, 7.499999999999999, 9.0, 11.25, 14.999999999999998, 22.5, 45.0]), SpeedyWeather.RingGrids.AnvilLocator{Float32}(10, [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, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], Float32[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], Float32[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], Float32[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]))
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.3652828935637102
+ 0.5565613842240793
+ 0.6862341194002235
+ 0.5322939678805791
+ 0.16473156242942733
+ 0.5775924903748594
+ 0.5113146641130806
+ 0.6264092866151043
+ 0.32267942942789263
+ 0.3210232072918106
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.3652828952713613
+ 0.5565613838092105
+ 0.6862341205332358
+ 0.5322939694019978
+ 0.16473156616212448
+ 0.577592495239229
+ 0.5113146569167053
+ 0.6264092866208574
+ 0.3226794245424488
+ 0.321023207851591
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
.
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.
abstract type AbstractFullGrid{T} <: AbstractGrid{T} end
An AbstractFullGrid
is a horizontal grid with a constant number of longitude points across latitude rings. Different latitudes can be used, Gaussian latitudes, equi-angle latitdes, or others.
sourceabstract type AbstractGrid{T} <: AbstractVector{T} end
The abstract supertype for all spatial grids on the sphere supported by SpeedyWeather.jl. Every new grid has to be of the form
abstract type AbstractGridClass{T} <: AbstractGrid{T} end
+struct MyNewGrid{T} <: AbstractGridClass{T}
+ data::Vector{T} # all grid points unravelled into a vector
+ nlat_half::Int # resolution: latitude rings on one hemisphere (Equator incl)
+end
MyNewGrid
should belong to a grid class like AbstractFullGrid
, AbstractOctahedralGrid
or AbstractHEALPixGrid
(that already exist but you may introduce a new class of grids) that share certain features such as the number of longitude points per latitude ring and indexing, but may have different latitudes or offset rotations. Each new grid Grid
(or grid class) then has to implement the following methods (as an example, see octahedral.jl)
Fundamental grid properties getnpoints # total number of grid points nlatodd # does the grid have an odd number of latitude rings? getnlat # total number of latitude rings getnlat_half # number of latitude rings on one hemisphere incl Equator
Indexing getnlonmax # maximum number of longitudes points (at the Equator) getnlonperring # number of longitudes on ring j eachindexinring # a unit range that indexes all longitude points on a ring
Coordinates getcolat # vector of colatitudes (radians) getcolatlon # vectors of colatitudes, longitudes (both radians)
Spectral truncation truncationorder # linear, quadratic, cubic = 1,2,3 for grid gettruncation # spectral truncation given a grid resolution get_resolution # grid resolution given a spectral truncation
Quadrature weights and solid angles getquadratureweights # = sinθ Δθ for grid points on ring j for meridional integration getsolidangle # = sinθ Δθ Δϕ, solid angle of grid points on ring j
sourceabstract type AbstractHEALPixGrid{T} <: AbstractGrid{T} end
An AbstractHEALPixGrid
is a horizontal grid similar to the standard HEALPixGrid, but different latitudes can be used, the default HEALPix latitudes or others.
sourceabstract type AbstractInterpolator{NF,G} end
Supertype for Interpolators. Every Interpolator <: AbstractInterpolator is expected to have two fields,
- geometry, which describes the grid G to interpolate from
- locator, which locates the indices on G and their weights to interpolate onto a new grid.
NF is the number format used to calculate the interpolation, which can be different from the input data and/or the interpolated data on the new grid.
sourceAbstractLocator{NF}
Supertype of every Locator, which locates the indices on a grid to be used to perform an interpolation. E.g. AnvilLocator uses a 4-point stencil for every new coordinate to interpolate onto. Higher order stencils can be implemented by defining OtherLocator <: AbstractLocactor.
sourceabstract type AbstractOctaHEALPixGrid{T} <: AbstractGrid{T} end
An AbstractOctaHEALPixGrid
is a horizontal grid similar to the standard OctahedralGrid, but the number of points in the ring closest to the Poles starts from 4 instead of 20, and the longitude of the first point in each ring is shifted as in HEALPixGrid. Also, different latitudes can be used.
sourceabstract type AbstractOctahedralGrid{T} <: AbstractGrid{T} end
An AbstractOctahedralGrid
is a horizontal grid with 16+4i longitude points on the latitude ring i starting with i=1 around the pole. Different latitudes can be used, Gaussian latitudes, equi-angle latitdes, or others.
sourceAnvilLocator{NF<:AbstractFloat} <: AbtractLocator
Contains arrays that locates grid points of a given field to be uses in an interpolation and their weights. This Locator is a 4-point average in an anvil-shaped grid-point arrangement between two latitude rings.
sourceL = AnvilLocator( ::Type{NF}, # number format used for the interpolation
+ npoints::Integer # number of points to interpolate onto
+ ) where {NF<:AbstractFloat}
Zero generator function for the 4-point average AnvilLocator. Use update_locator! to update the grid indices used for interpolation and their weights. The number format NF is the format used for the calculations within the interpolation, the input data and/or output data formats may differ.
sourceG = FullClenshawGrid{T}
A FullClenshawGrid is a regular latitude-longitude grid with an odd number of nlat
equi-spaced latitudes, the central latitude ring is on the Equator. The same nlon
longitudes for every latitude ring. The grid points are closer in zonal direction around the poles. The values of all grid points are stored in a vector field data
that unravels the data 0 to 360˚, then ring by ring, which are sorted north to south.
sourceG = FullGaussianGrid{T}
A full Gaussian grid is a regular latitude-longitude grid that uses nlat
Gaussian latitudes, and the same nlon
longitudes for every latitude ring. The grid points are closer in zonal direction around the poles. The values of all grid points are stored in a vector field v
that unravels the data 0 to 360˚, then ring by ring, which are sorted north to south.
sourceG = FullHEALPixGrid{T}
A full HEALPix grid is a regular latitude-longitude grid that uses nlat
latitudes from the HEALPix grid, and the same nlon
longitudes for every latitude ring. The grid points are closer in zonal direction around the poles. The values of all grid points are stored in a vector field v
that unravels the data 0 to 360˚, then ring by ring, which are sorted north to south.
sourceG = FullOctaHEALPixGrid{T}
A full OctaHEALPix grid is a regular latitude-longitude grid that uses nlat
OctaHEALPix latitudes, and the same nlon
longitudes for every latitude ring. The grid points are closer in zonal direction around the poles. The values of all grid points are stored in a vector field v
that unravels the data 0 to 360˚, then ring by ring, which are sorted north to south.
sourceGridGeometry{G<:AbstractGrid}
contains general precomputed arrays describing the grid of G.
sourceG = GridGeometry( Grid::Type{<:AbstractGrid},
+ nlat_half::Integer)
Precomputed arrays describing the geometry of the Grid with resolution nlat_half. Contains latitudes and longitudes of grid points, their ring index j and their unravelled indices ij.
sourceH = HEALPixGrid{T}
A HEALPix grid with 12 faces, each nside
xnside
grid points, each covering the same area. The number of latitude rings on one hemisphere (incl Equator) nlat_half
is used as resolution parameter. The values of all grid points are stored in a vector field v
that unravels the data 0 to 360˚, then ring by ring, which are sorted north to south.
sourceH = OctaHEALPixGrid{T}
A OctaHEALPix grid with 4 base faces, each nlat_half
xnlat_half
grid points, each covering the same area. The values of all grid points are stored in a vector field data
that unravels the data 0 to 360˚, then ring by ring, which are sorted north to south.
sourceG = OctahedralClenshawGrid{T}
An Octahedral Clenshaw grid that uses nlat
equi-spaced latitudes. Like FullClenshawGrid, the central latitude ring is on the Equator. Like OctahedralGaussianGrid, the number of longitude points per latitude ring decreases towards the poles. Starting with 20 equi-spaced longitude points (starting at 0˚E) on the rings around the poles, each latitude ring towards the equator has consecuitively 4 more points, one for each face of the octahedron. E.g. 20,24,28,32,...nlon-4,nlon,nlon,nlon-4,...,32,28,24,20. The maximum number of longitue points is nlon
. The values of all grid points are stored in a vector field v
that unravels the data 0 to 360˚, then ring by ring, which are sorted north to south.
sourceG = OctahedralGaussianGrid{T}
An Octahedral Gaussian grid that uses nlat
Gaussian latitudes, but a decreasing number of longitude points per latitude ring towards the poles. Starting with 20 equi-spaced longitude points (starting at 0˚E) on the rings around the poles, each latitude ring towards the equator has consecuitively 4 more points, one for each face of the octahedron. E.g. 20,24,28,32,...nlon-4,nlon,nlon,nlon-4,...,32,28,24,20. The maximum number of longitue points is nlon
. The values of all grid points are stored in a vector field v
that unravels the data 0 to 360˚, then ring by ring, which are sorted north to south.
sourceMatrix!(M::AbstractMatrix,
+ G::OctaHEALPixGrid;
+ quadrant_rotation=(0,1,2,3),
+ matrix_quadrant=((2,2),(1,2),(1,1),(2,1)),
+ )
Sorts the gridpoints in G
into the matrix M
without interpolation. Every quadrant of the grid G
is rotated as specified in quadrant_rotation
, 0 is no rotation, 1 is 90˚ clockwise, 2 is 180˚ etc. Grid quadrants are counted eastward starting from 0˚E. The grid quadrants are moved into the matrix quadrant (i,j) as specified. Defaults are equivalent to centered at 0˚E and a rotation such that the North Pole is at M's midpoint.
sourceMatrix!(M::AbstractMatrix,
+ G::OctahedralClenshawGrid;
+ quadrant_rotation=(0,1,2,3),
+ matrix_quadrant=((2,2),(1,2),(1,1),(2,1)),
+ )
Sorts the gridpoints in G
into the matrix M
without interpolation. Every quadrant of the grid G
is rotated as specified in quadrant_rotation
, 0 is no rotation, 1 is 90˚ clockwise, 2 is 180˚ etc. Grid quadrants are counted eastward starting from 0˚E. The grid quadrants are moved into the matrix quadrant (i,j) as specified. Defaults are equivalent to centered at 0˚E and a rotation such that the North Pole is at M's midpoint.
sourceMatrix!(MGs::Tuple{AbstractMatrix{T},OctaHEALPixGrid}...;kwargs...)
Like Matrix!(::AbstractMatrix,::OctaHEALPixGrid)
but for simultaneous processing of tuples ((M1,G1),(M2,G2),...)
with matrices Mi
and grids Gi
. All matrices and grids have to be of the same size respectively.
sourceMatrix!(MGs::Tuple{AbstractMatrix{T},OctahedralClenshawGrid}...;kwargs...)
Like Matrix!(::AbstractMatrix,::OctahedralClenshawGrid)
but for simultaneous processing of tuples ((M1,G1),(M2,G2),...)
with matrices Mi
and grids Gi
. All matrices and grids have to be of the same size respectively.
sourceanvil_average(a, b, c, d, Δab, Δcd, Δy) -> Any
+
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. See schematic:
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.
sourceaverage_on_poles(
+ A::SpeedyWeather.RingGrids.AbstractGrid{NF<:AbstractFloat},
+ rings::Vector{<:UnitRange{<:Integer}}
+) -> Tuple{Any, Any}
+
Computes the average at the North and South pole from a given grid A
and it's precomputed ring indices rings
. The North pole average is an equally weighted average of all grid points on the northern-most ring. Similar for the South pole.
sourceaverage_on_poles(
+ A::SpeedyWeather.RingGrids.AbstractGrid{NF<:Integer},
+ rings::Vector{<:UnitRange{<:Integer}}
+) -> Tuple{Any, Any}
+
Method for A::Abstract{T<:Integer}
which rounds the averaged values to return the same number format NF
.
sourcei = each_index_in_ring(grid,j)
UnitRange i
to access data on grid grid
on ring j
.
sourceijs = eachgridpoint(grid)
UnitRange ijs
to access each grid point on grid grid
.
sourceeachring(grid::SpeedyWeather.RingGrids.AbstractGrid) -> Any
+
Vector{UnitRange} rings
to loop over every ring of grid grid
and then each grid point per ring. To be used like
rings = eachring(grid)
+for ring in rings
+ for ij in ring
+ grid[ij]
sourceeachring(
+ grid1::SpeedyWeather.RingGrids.AbstractGrid,
+ grids::Grid<:SpeedyWeather.RingGrids.AbstractGrid...
+) -> Any
+
Same as eachring(grid)
but performs a bounds check to assess that all grids in grids
are of same size.
sourcetrue/false = extrema_in(v::Vector,a::Real,b::Real)
For every element vᵢ in v does a<=vi<=b hold?
sourceget_nlons(
+ Grid::Type{<:SpeedyWeather.RingGrids.AbstractGrid},
+ nlat_half::Integer;
+ both_hemispheres
+) -> Any
+
Returns a vector nlons
for the number of longitude points per latitude ring, north to south. Provide grid Grid
and its resolution parameter nlat_half
. For both_hemisphere==false only the northern hemisphere (incl Equator) is returned.
sourcetrue/false = isdecreasing(v::Vector)
Check whether elements of a vector v
are strictly decreasing.
sourcetrue/false = isincreasing(v::Vector)
Check whether elements of a vector v
are strictly increasing.
sourcei_new,j_new = rotate_matrix_indices_180(i,j,s)
Rotate indices i,j
of a square matrix of size s x s by 180˚.
sourcei_new,j_new = rotate_matrix_indices_270(i,j,s)
Rotate indices i,j
of a square matrix of size s x s anti-clockwise by 270˚.
sourcei_new,j_new = rotate_matrix_indices_90(i,j,s)
Rotate indices i,j
of a square matrix of size s x s anti-clockwise by 90˚.
sourcewhichring(
+ ij::Integer,
+ rings::Vector{UnitRange{Int64}}
+) -> Int64
+
Obtain ring index j from gridpoint ij and Vector{UnitRange} describing rind indices as obtained from eachring(::Grid)
source