Remapping

interpolate_array(field, xpts, ypts)
interpolate_array(field, xpts, ypts, zpts)

Interpolate a field to a regular array using pointwise interpolation.

This is primarily used for plotting and diagnostics.

Examples

longpts = range(Geometry.LongPoint(-180.0), Geometry.LongPoint(180.0), length = 21)
latpts = range(Geometry.LatPoint(-80.0), Geometry.LatPoint(80.0), length = 21)
zpts = range(Geometry.ZPoint(0.0), Geometry.ZPoint(1000.0), length = 21)

interpolate_array(field, longpts, latpts, zpts)

interpolate(remapper::Remapper, fields) interpolate!(dest, remapper::Remapper, fields)

Interpolate the given field(s) as prescribed by remapper.

The optimal number of fields passed is the buffer_length of the remapper. If more fields are passed, the remapper will batch work with size up to its buffer_length.

This call mutates the internal (private) state of the remapper.

Horizontally, interpolation is performed with the barycentric formula in [14], equation (3.2). Vertical interpolation is linear except in the boundary elements where it is 0th order.

interpolate! writes the output to the given destiniation. dest is expected to be defined on the root process and to be nothing for the other processes.

Note: interpolate allocates new arrays and has some internal type-instability, interpolate! is non-allocating and type-stable.

When using interpolate!, the destination has to be the same array type as the device in use (e.g., CuArray for CUDA runs).

Example

Given field1,field2, two Field defined on a cubed sphere.

longpts = range(-180.0, 180.0, 21)
latpts = range(-80.0, 80.0, 21)
zpts = range(0.0, 1000.0, 21)

hcoords = [Geometry.LatLongPoint(lat, long) for long in longpts, lat in latpts]
zcoords = [Geometry.ZPoint(z) for z in zpts]

space = axes(field1)

remapper = Remapper(space, hcoords, zcoords)

int1 = interpolate(remapper, field1)
int2 = interpolate(remapper, field2)

# Or
int12 = interpolate(remapper, [field1, field2])
# With int1 = int12[1, :, :, :]

   interpolate(field::ClimaCore.Fields;
               hresolution = 180,
               zresolution = nothing,
               target_hcoords = default_target_hcoords(space; hresolution),
               target_zcoords = default_target_zcoords(space; zresolution)
               )

Interpolate the given fields on the Cartesian product of target_hcoords with target_zcoords. For Spaces without horizontal/vertical component, the relevant target_*coords are ignored.

When zresolution is set to nothing, do not perform any interpolation in the vertical direction. When zresolution, perform linear interpolation with uniformly spaced levels.

Coordinates have to be ClimaCore.Geometry.Points.

Note: do not use this method when performance is important. Instead, define a Remapper and call interpolate(remapper, fields). Different Fields defined on the same Space can share a Remapper, so that interpolation can be optimized.

Example

Given field, a Field defined on a cubed sphere.

By default, a target uniform grid is chosen (with resolution hresolution and zresolution), so remapping is simply

julia> interpolate(field)

This will return an array of interpolated values.

Resolution can be specified

julia> interpolate(field; hresolution = 100, zresolution = 50)

Coordinates can be also specified directly:

julia> longpts = range(-180.0, 180.0, 21)
julia> latpts = range(-80.0, 80.0, 21)
julia> zpts = range(0.0, 1000.0, 21)

julia> hcoords = [Geometry.LatLongPoint(lat, long) for long in longpts, lat in latpts]
julia> zcoords = [Geometry.ZPoint(z) for z in zpts]

julia> interpolate(field, target_hcoords, target_zcoords)

If you need the array of coordinates, you can call default_target_hcoords (or default_target_zcoords) passing axes(field). This will return an array of Geometry.Points. The functions Geometry.Components and Geometry.Component can be used to extract the components as numeric values. For example,

julia> Geometry.components.(Geometry.components.([
           Geometry.LatLongPoint(x, y) for x in range(-180.0, 180.0, length = 180),
           y in range(-90.0, 90.0, length = 180)
       ]))
180×180 Matrix{StaticArraysCore.SVector{2, Float64}}:
 [-180.0, -90.0]    [-180.0, -88.9944]    …  [-180.0, 88.9944]    [-180.0, 90.0]
  ⋮                                        ⋱
 [180.0, -90.0]     [180.0, -88.9944]        [180.0, 88.9944]     [180.0, 90.0]

To extract only long or lat, one can broadcast getindex

julia> lats = getindex.(Geometry.components.([Geometry.LatLongPoint(x, y)
                                              for x in range(-180.0, 180.0, length = 180),
                                                  y in range(-90.0, 90.0, length = 180)
                                             ]),
                        1)

This can be used directly for plotting.