Add XYZPoint regridding interpolation

This commit adds support for `ClimaCore.Geometry.XYZPoint` to
`InterpolationsRegridder` which was needed in order to support
interpolation of observational data used for simulation forcing to
single column models in ClimaAtmos. Note that because single column
models are currently represented as boxes we needed support for XYZPoint
and not ZPoint objects. As such data that is read into these objects
must be expanded to include ``dummy" x and y dimensions. More generally,
this addition supports interpolations over boxes.

To be conservative by default we throw an error if extrapolating outside
of the Z range of the data. For ERA5 however this might need to be
changed to `Linear` extrapolation as typically simulations will have
points closer to the ground than the lowest point in ERA5.

Co-authored-by: Julian Schmitt <[email protected]>
Co-authored-by: Gabriele Bozzola <[email protected]>

Author: Julians42

NEWS.md

@@ -3,6 +3,16 @@ ClimaUtilities.jl Release Notes
 
 main
 ------
+
+v0.1.23
+------
+
+#### Support for Box Interpolations Regridder. PR[#151](https://github.com/CliMA/ClimaUtilities.jl/pull/151)
+
+`Regridders.IterpolationsRegridder` now supports regridding on 
+`ClimaCore.Geometry.XYZPoint` objects which allows for interpolation 
+onto boxes and single column simulations. 
+
 v0.1.22
 ------
 

Project.toml

@@ -1,7 +1,7 @@
 name = "ClimaUtilities"
 uuid = "b3f4f4ca-9299-4f7f-bd9b-81e1242a7513"
 authors = ["Gabriele Bozzola <[email protected]>", "Julia Sloan <[email protected]>"]
-version = "0.1.22"
+version = "0.1.23"
 
 [deps]
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

docs/src/regridders.md

@@ -12,7 +12,7 @@ Currently, `Regridders` comes with two implementations:
    directly with files.
 2. `InterpolationsRegridder` uses
    [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl) to
-   perform non-conservative linear interpolation onto lat-long(-z) grids.
+   perform non-conservative linear interpolation onto lat-long(-z) and x-y-z grids.
    `InterpolationsRegridder` works directly with data.
 
 !!! note

ext/InterpolationsRegridderExt.jl

@@ -33,6 +33,7 @@ end
 # should be first, so files will have longitude as first dimension.
 totuple(pt::ClimaCore.Geometry.LatLongZPoint) = pt.long, pt.lat, pt.z
 totuple(pt::ClimaCore.Geometry.LatLongPoint) = pt.long, pt.lat
+totuple(pt::ClimaCore.Geometry.XYZPoint) = pt.x, pt.y, pt.z
 
 """
     InterpolationsRegridder(target_space::ClimaCore.Spaces.AbstractSpace
@@ -77,8 +78,12 @@ function Regridders.InterpolationsRegridder(
         isnothing(extrapolation_bc) &&
             (extrapolation_bc = (Intp.Periodic(), Intp.Flat(), Intp.Throw()))
         isnothing(dim_increasing) && (dim_increasing = (true, true, true))
+    elseif eltype(coordinates) <: ClimaCore.Geometry.XYZPoint
+        isnothing(extrapolation_bc) &&
+            (extrapolation_bc = (Intp.Flat(), Intp.Flat(), Intp.Throw()))
+        isnothing(dim_increasing) && (dim_increasing = (true, true, true))
     else
-        error("Only lat-long, lat-long-z spaces are supported")
+        error("Only lat-long, lat-long-z, and x-y-z spaces are supported")
     end
 
     return InterpolationsRegridder(

test/TestTools.jl

@@ -1,5 +1,7 @@
 import ClimaCore
 import ClimaComms
+import ClimaCore: CommonSpaces, Grids
+
 @static pkgversion(ClimaComms) >= v"0.6" && ClimaComms.@import_required_backends
 
 function make_spherical_space(FT; context = ClimaComms.context())
@@ -92,3 +94,51 @@ function make_regional_space(FT; context = ClimaComms.context())
         vertical = vert_center_space,
     )
 end
+
+function make_box_space(FT; context = ClimaComms.context())
+    helem = (10, 10)
+    xrange = (0.0, 1.0)
+    yrange = (0.0, 1.0)
+    zrange = (0.0, 1.0)
+    Nq = 4
+    zelem = 10
+
+    vertical_domain = ClimaCore.Domains.IntervalDomain(
+        ClimaCore.Geometry.ZPoint{FT}(zrange[1]),
+        ClimaCore.Geometry.ZPoint{FT}(zrange[2]),
+        boundary_names = (:bottom, :top),
+    )
+
+    x_domain = ClimaCore.Domains.IntervalDomain(
+        ClimaCore.Geometry.XPoint{FT}(xrange[1]),
+        ClimaCore.Geometry.XPoint{FT}(xrange[2]),
+        boundary_names = (:west, :east),
+    )
+
+    y_domain = ClimaCore.Domains.IntervalDomain(
+        ClimaCore.Geometry.YPoint{FT}(yrange[1]),
+        ClimaCore.Geometry.YPoint{FT}(yrange[2]),
+        boundary_names = (:south, :north),
+    )
+
+    horzdomain = ClimaCore.Domains.RectangleDomain(x_domain, y_domain)
+    horzmesh = ClimaCore.Meshes.RectilinearMesh(horzdomain, helem...)
+    horztopology = ClimaCore.Topologies.Topology2D(context, horzmesh)
+    quad = ClimaCore.Spaces.Quadratures.GLL{Nq}()
+    horzspace = ClimaCore.Spaces.SpectralElementSpace2D(horztopology, quad)
+
+    vertmesh = ClimaCore.Meshes.IntervalMesh(vertical_domain, nelems = zelem)
+
+    vert_center_space = ClimaCore.Spaces.CenterFiniteDifferenceSpace(
+        ClimaComms.device(context),
+        vertmesh,
+    )
+
+    hybrid = ClimaCore.Spaces.ExtrudedFiniteDifferenceSpace(
+        horzspace,
+        vert_center_space,
+    )
+
+    return hybrid
+
+end

test/regridders.jl

@@ -215,6 +215,50 @@ end
     end
 end
 
+@testset "InterpolationsRegridderXYZPoint" begin
+    helem = (10, 10)
+    Nq = 4
+    zelem = 10
+    FT = Float64 # test for Float64
+
+    x, y, z = collect(0.0:1:5), collect(0.0:1:6), collect(0.0:1:7)
+
+    dimensions3D = (x, y, z)
+    size3D = (6, 7, 8)
+
+    data_x3D = zeros(size3D)
+    data_y3D = zeros(size3D)
+    data_z3D = zeros(size3D)
+
+    for i in 1:length(x)
+        for j in 1:length(y)
+            for k in 1:length(z)
+                data_x3D[i, j, k] = x[i]
+                data_y3D[i, j, k] = y[j]
+                data_z3D[i, j, k] = z[k]
+            end
+        end
+    end
+    # create the box space with helper function
+    spaces = make_box_space(Float64; context)
+    # create interpolation object with InterpolationsRegridder for XYZPoint object
+    reg_box = Regridders.InterpolationsRegridder(spaces)
+
+    regridded_x = Regridders.regrid(reg_box, data_x3D, dimensions3D)
+
+    regridded_y = Regridders.regrid(reg_box, data_y3D, dimensions3D)
+
+    regridded_z = Regridders.regrid(reg_box, data_z3D, dimensions3D)
+
+    err_x = reg_box.coordinates.x .- regridded_x
+    err_y = reg_box.coordinates.y .- regridded_y
+    err_z = reg_box.coordinates.z .- regridded_z
+
+    @test maximum(err_x) < 1e-5
+    @test maximum(err_y) < 1e-5
+    @test maximum(err_z) < 1e-5
+end
+
 @testset "TempestRegridder" begin
     for FT in (Float32, Float64)
         data_path = joinpath(@__DIR__, "test_data", "era5_1979_1.0x1.0_lai.nc")