export ClimatologyInterpolator

struct ClimatologyInterpolator{I<:OneDimensionalInterpolator,ğ�’¯}
    x::Vector{ğ�’¯}
    mask::BitMatrix
    ğ�’»r::Matrix{I}
    ğ�’»T::Matrix{I}
    A::Matrix{ğ�’¯}
    f::Matrix{ğ�’¯}
    lat::Vector{ğ�’¯}
    n::Int64
    m::Int64
    L::Int64
end

function Base.show(io::IO, â„�::ClimatologyInterpolator{I,ğ�’¯}) where {I,ğ�’¯}
    @unpack n, m = â„�
    print(io, "ClimatologyInterpolator{$I,$ğ�’¯}, $n x $m")
end

Base.size(â„�::ClimatologyInterpolator) = (â„�.n, â„�.m, â„�.L)

function ClimatologyInterpolator(ğ�’�::AbstractVector{Climatology{ğ�’¯}},
                                 x::AbstractVector{<:Real},
                                 interpolator::Type=CubicSplineInterpolator,
                                 boundaries::Type=StrictBoundaries) where {ğ�’¯}
    @assert issorted(x) "x vector must be sorted in ascending order"
    @assert length(ğ�’�) > 1 "must have at least two Climatologies"
    n, m = size(ğ�’�[1])
    mask = ğ�’�[1].mask
    A = ğ�’�[1].A
    f = ğ�’�[1].f
    lat = ğ�’�[1].lat
    for ğ�’¸ âˆˆ ğ�’�
        #demand identical size
        @assert size(ğ�’¸) == (n,m) "Climatologies must all be the same size"
        #demand identical cell areas
        @assert all(ğ�’¸.A .â‰ˆ A) "cell area grids must be identical"
        #demand identical cell land fractions
        @assert all(ğ�’¸.f .â‰ˆ f) "cell land fraction grids must be identical"
        #demand identical latitudes
        @assert all(ğ�’¸.lat .â‰ˆ lat) "cell latitudes must be the same"
        #union of all masked cells
        mask .*= ğ�’¸.mask
    end
    #construct interpolators
    x = collect(ğ�’¯, x)
    @multiassign ğ�’»r, ğ�’»T = Matrix{interpolator}(undef, n, m)
    for i âˆˆ 1:n, j âˆˆ 1:m
        if mask[i,j]
            ğ�’»r[i,j] = interpolator(x, map(ğ�’¸->ğ�’¸.r[i,j], ğ�’�), boundaries())
            ğ�’»T[i,j] = interpolator(x, map(ğ�’¸->ğ�’¸.T[i,j], ğ�’�), boundaries())
        end
    end
    #construct unified interpolator
    ClimatologyInterpolator(x, mask, ğ�’»r, ğ�’»T, A, f, lat, n, m, length(ğ�’�))
end

function (â„�::ClimatologyInterpolator{I,ğ�’¯})(x) where {I,ğ�’¯}
    #pull out fields of struct
    @unpack mask, ğ�’»r, ğ�’»T, A, f, lat, n, m = â„�
    #assume NaN until unmasked
    @multiassign r, T = fill(convert(ğ�’¯, NaN), (n, m))
    #interpolate runoff and temperature at all points
    @inbounds for i âˆˆ 1:n, j âˆˆ 1:m
        if mask[i,j]
            #interpolate runoff
            r[i,j] = ğ�’»r[i,j](x)
            #interpolate temperature
            T[i,j] = ğ�’»T[i,j](x)
        end
    end
    #construct a new Climatology
    Climatology(mask, r, T, A, f, lat, n, m)
end
