Expose interpolation method to users

Author: ph-kev

NEWS.md

@@ -4,6 +4,26 @@ ClimaUtilities.jl Release Notes
 main
 ------
 
+#### Interpolation method for InterpolationsRegridder
+
+`InterpolationsRegridder` now accepts the new keyword argument
+`interpolation_method` for choosing which kind of interpolation to do. The
+default is `Interpolations.Linear()`.
+
+```julia
+import ClimaUtilities.Regridders
+import ClimaCore, Interpolations
+
+linear_reg = Regridders.InterpolationsRegridder(
+    target_space,
+    interpolation_method = Intp.Linear(),
+)
+constant_reg = Regridders.InterpolationsRegridder(
+    target_space,
+    interpolation_method = Intp.Constant(),
+)
+```
+
 v0.1.23
 ------
 

docs/src/regridders.md

@@ -65,6 +65,13 @@ regridded_u = Regridders.regrid(reg, target_date)
 along each direction) and returns a `ClimaCore` `Field` defined on the
 `target_space`.
 
+!!! note "Did you know?"
+    With versions of ClimaUtilities after v0.1.24, you can change the
+    interpolation method to perform any gridded interpolation supported by
+    Interpolations.jl. For example, to do constant interpolation, you can pass
+    `interpolation_method = Interpolations.Constant()` as a keyword argument
+    when constructing the regridder.
+
 Currently, `InterpolationsRegridder` only supports spherical shells and extruded
 spherical shells (but it could be easily extended to other domains).
 

ext/InterpolationsRegridderExt.jl

@@ -11,6 +11,7 @@ import ClimaUtilities.Regridders
 struct InterpolationsRegridder{
     SPACE <: ClimaCore.Spaces.AbstractSpace,
     FIELD <: ClimaCore.Fields.Field,
+    IM,
     BC,
     DT <: Tuple,
 } <: Regridders.AbstractRegridder
@@ -21,6 +22,9 @@ struct InterpolationsRegridder{
     """ClimaCore.Field of physical coordinates over which the data will be interpolated"""
     coordinates::FIELD
 
+    """Method of interpolation as accepted by Interpolations.jl"""
+    interpolation_method::IM
+
     """Tuple of extrapolation conditions as accepted by Interpolations.jl"""
     extrapolation_bc::BC
 
@@ -37,7 +41,9 @@ totuple(pt::ClimaCore.Geometry.XYZPoint) = pt.x, pt.y, pt.z
 
 """
     InterpolationsRegridder(target_space::ClimaCore.Spaces.AbstractSpace
-                            [; extrapolation_bc::Tuple])
+                            [; extrapolation_bc::Tuple,
+                               dim_increasing::Union{Nothing, Tuple},
+                               interpolation_method = Interpolations.Linear()])
 
 An online regridder that uses Interpolations.jl
 
@@ -67,6 +73,7 @@ function Regridders.InterpolationsRegridder(
     target_space::ClimaCore.Spaces.AbstractSpace;
     extrapolation_bc::Union{Nothing, Tuple} = nothing,
     dim_increasing::Union{Nothing, Tuple} = nothing,
+    interpolation_method = Intp.Linear(),
 )
     coordinates = ClimaCore.Fields.coordinate_field(target_space)
     # set default values for the extrapolation_bc and dim_increasing if they are not provided
@@ -89,6 +96,7 @@ function Regridders.InterpolationsRegridder(
     return InterpolationsRegridder(
         target_space,
         coordinates,
+        interpolation_method,
         extrapolation_bc,
         dim_increasing,
     )
@@ -119,7 +127,7 @@ function Regridders.regrid(regridder::InterpolationsRegridder, data, dimensions)
         Intp.interpolate(
             dimensions_FT,
             FT.(data_transformed),
-            Intp.Gridded(Intp.Linear()),
+            Intp.Gridded(regridder.interpolation_method),
         ),
         regridder.extrapolation_bc,
     )
@@ -148,7 +156,11 @@ function Regridders.regrid!(
         "Dimensions must be monotonically increasing to use regrid!. Sort the dimensions first, or use regrid.",
     )
     itp = Intp.extrapolate(
-        Intp.interpolate(dimensions, data, Intp.Gridded(Intp.Linear())),
+        Intp.interpolate(
+            dimensions,
+            data,
+            Intp.Gridded(regridder.interpolation_method),
+        ),
         regridder.extrapolation_bc,
     )
     gpuitp = Adapt.adapt(ClimaComms.array_type(regridder.target_space), itp)

test/regridders.jl

@@ -267,6 +267,40 @@ end
     end
 end
 
+@testset "Interpolation method" begin
+    # Test constant interpolation method
+    lon, lat, z =
+        collect(0.0:1:360), collect(-90.0:1:90), collect(0.0:1.0:100.0)
+    dimensions2D = (lon, lat)
+    dimensions3D = (lon, lat, z)
+    size2D = (361, 181)
+    size3D = (361, 181, 101)
+    ones_2Ddata = ones(size2D)
+    ones_3Ddata = ones(size3D)
+
+    for FT in (Float32, Float64)
+        spaces = make_spherical_space(FT; context)
+        horzspace = spaces.horizontal
+        hv_center_space = spaces.hybrid
+
+        reg_horz = Regridders.InterpolationsRegridder(
+            horzspace,
+            interpolation_method = Interpolations.Constant(),
+        )
+        reg_hv = Regridders.InterpolationsRegridder(
+            hv_center_space,
+            interpolation_method = Interpolations.Constant(),
+        )
+
+        regridded_2Ddata =
+            Regridders.regrid(reg_horz, ones_2Ddata, dimensions2D)
+        regridded_3Ddata = Regridders.regrid(reg_hv, ones_3Ddata, dimensions3D)
+
+        @test all(x -> x == one(x), parent(regridded_2Ddata))
+        @test all(x -> x == one(x), parent(regridded_3Ddata))
+    end
+end
+
 @testset "InterpolationsRegridderXYZPoint" begin
     helem = (10, 10)
     Nq = 4