New sametimespan function and some bugfixes. (#76)

* fix ambiguity bug in zonal mean

* impelemnt tropics_extratropics for Coord data

* correct sametimespan function

* aDD TESTS for sametimespan

* bump version

* fix rtests

Project.toml

@@ -1,7 +1,7 @@
 name = "ClimateBase"
 uuid = "35604d93-0fb8-4872-9436-495b01d137e2"
 authors = ["Datseris <[email protected]>", "Philippe Roy <[email protected]>"]
-version = "0.13.6"
+version = "0.13.7"
 
 [deps]
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"

docs/src/index.md

@@ -162,6 +162,7 @@ temporal_sampling
 realtime_days
 realtime_milliseconds
 seasonality
+sametimespan
 ```
 
 ## Spatial

src/ClimateBase.jl

@@ -17,4 +17,6 @@ include("climate/albedo.jl")
 include("tsa/continuation.jl")
 include("tsa/decomposition.jl")
 
+include("exports.jl")
+
 end # module

src/exports.jl

@@ -0,0 +1,12 @@
+# This file contains the exports for ClimateBase module.
+# Having all exports at a single file is overall cleaner
+# and leads to easier-to-identify public API.
+
+# TODO: Actually move all exports here.
+
+# Temporal
+export monthday_indices, maxyearspan, daymonth, realtime_days, realtime_milliseconds,
+temporal_sampling, timemean, timeagg, monthlyagg, yearlyagg, temporalrange, seasonalyagg, 
+season, DAYS_IN_ORBIT, HOURS_IN_ORBIT, seasonality, sametimespan
+
+#

src/physical_dimensions/spatial_equalarea.jl

@@ -27,7 +27,7 @@ function zonalmean(::UnstructuredGrid, A::AbDimArray{T, 1}, ::Nothing) where {T}
 end
 
 # zonal mean with weights
-function zonalmean(::UnstructuredGrid, A::ClimArray, W)
+function zonalmean(::UnstructuredGrid, A::ClimArray, W::AbstractArray)
 	@assert size(A) == size(W)
     idxs, lats = uniquelats(A)
     other = otherdims(A, Coord())
@@ -40,7 +40,7 @@ function zonalmean(::UnstructuredGrid, A::ClimArray, W)
     end
     return R
 end
-function zonalmean(::UnstructuredGrid, A::ClimArray{T, 1}, W) where {T}
+function zonalmean(::UnstructuredGrid, A::ClimArray{T, 1}, W::AbstractArray) where {T}
 	@assert size(A) == size(W)
     idxs, lats = uniquelats(A)
     res = zeros(T, length(lats))
@@ -52,7 +52,7 @@ end
 
 
 """
-    uniquelats(A::AbDimArray) → idxs, lats
+    uniquelats(A::ClimArray) → idxs, lats
     uniquelats(c::Vector{<:AbstractVector}) → idxs, lats
 Find the unique latitudes of `A`. Return the indices (vector of ranges) that each latitude
 in `lats` covers, as well as the latitudes themselves.
@@ -121,13 +121,33 @@ Return the indices of coordinates belonging to the north and south hemispheres.
 function hemisphere_indices(c)
     idxs, lats = uniquelats(c)
     i = findfirst(x -> x > 0, lats)
-    shi = idxs[1][1]:idxs[i][end]
-    nhi = idxs[i+1][1]:idxs[end][end]
+    shi = idxs[1][1]:idxs[i-1][end]
+    nhi = idxs[i][1]:idxs[end][end]
     return nhi, shi
 end
 
 latitudes(::UnstructuredGrid, A) = unique!([x[2] for x in dims(A, Coord)])
 
+function tropics_extratropics(::UnstructuredGrid, A; lower_lat=30)
+    # TODO: Support `higher_lat` keyword
+    c = dims(A, Coord).val
+    idxs, lats = uniquelats(c)
+    i1 = findlast(x -> x < -lower_lat, lats)
+    i2 = findfirst(x -> x > lower_lat, lats)
+    # tropics indices (accounting for hemispheric only data as well)
+    t1 = isnothing(i1) ? 0 : i1
+    t2 = isnothing(i2) ? length(idxs)+1 : i2
+    i_tropics = idxs[t1+1][1]:idxs[t2-1][end]
+    # extratropics indices (accounting for hemispheric only data as well)
+    i_sh_extra = isnothing(i1) ? (1:0) : idxs[1][1]:idxs[i1][end]
+    i_nh_extra = isnothing(i2) ? (1:0) : idxs[i2][1]:idxs[end][end]
+    i_extra = vcat(i_sh_extra, i_nh_extra)
+
+    tropics = A[Coord(i_tropics)]
+    extratropics = A[Coord(i_extra)]
+    return tropics, extratropics
+end
+
 #########################################################################
 # Extention of convenience indexing of `Coord`
 #########################################################################
@@ -145,4 +165,4 @@ end
 function DimensionalData.sel2indices(c::Coord, sel::Tuple{<:Lat{ <: Between}})
     l1, l2 = sel[1].val.val
     return coord_latitudes_between(c.val, l1, l2) # this is Vector{Int}
-end
+end

src/physical_dimensions/temporal.jl

@@ -2,12 +2,6 @@
 Handling of time in data as a physical quantity, and time-related data processing
 =#
 using Statistics, StatsBase
-export monthday_indices, maxyearspan, daymonth, realtime_days, realtime_milliseconds
-export temporal_sampling
-export timemean, timeagg
-export monthlyagg, yearlyagg, temporalrange, seasonalyagg, season
-export DAYS_IN_ORBIT, HOURS_IN_ORBIT
-export seasonality
 #########################################################################
 # Datetime related
 #########################################################################
@@ -62,12 +56,12 @@ function temporal_sampling(t::AbstractVector{<:TimeType})
         return :other
     end
 end
-temporal_sampling(t::AbstractVector) = :other
-temporal_sampling(t::StepRange{<:Any,Month}) = :monthly
-temporal_sampling(t::StepRange{<:Any,Year}) = :yearly
-temporal_sampling(t::StepRange{<:Any,Day}) = :daily
-temporal_sampling(t::StepRange{<:Any,Hour}) = :hourly
-temporal_sampling(t::StepRange{<:Any,<:Any}) = :other
+temporal_sampling(::AbstractVector) = :other
+temporal_sampling(::StepRange{<:Any,Month}) = :monthly
+temporal_sampling(::StepRange{<:Any,Year}) = :yearly
+temporal_sampling(::StepRange{<:Any,Day}) = :daily
+temporal_sampling(::StepRange{<:Any,Hour}) = :hourly
+temporal_sampling(::StepRange{<:Any,<:Any}) = :other
 
 "return the appropriate subtype of `Dates.Period` or `nothing`."
 function tsamp2period(tsamp)
@@ -214,6 +208,30 @@ end
 realtime_milliseconds(A) = realtime_milliseconds(dims(A, Ti).val)
 
 
+"""
+    sametimespan(Xs...) → Ys
+Given several `ClimArray`s, return the same `ClimArray`s but now accessed in the `Time`
+dimension so that they all have span the same time interval.
+
+`sametimespan` also has more intelligent handling of monthly or yearly sampled data.
+"""
+function sametimespan(Xs...)
+    mint = maximum(minimum(dims(X, Time).val) for X in Xs)
+    maxt = minimum(maximum(dims(X, Time).val) for X in Xs)
+    # Make an intelligent decision for monthly/yearly sampled data
+    tsamps = temporal_sampling.(Xs)
+    if all(isequal(:monthly), tsamps)
+        mint = Date(year(mint), month(mint), 1)
+        d = daysinmonth(maxt)
+        maxt = Date(year(maxt), month(maxt), d)
+    elseif all(isequal(:yearly), tsamps)
+        mint = Date(year(mint), 1, 1)
+        maxt = Date(year(maxt), 12, 31)
+    end
+    map(X -> X[Time(Between(mint, maxt))], Xs)
+end
+sametimespan(Xs::Tuple) = sametimespan(Xs...)
+
 #########################################################################
 # temporal statistics
 #########################################################################

test/runtests.jl

@@ -131,6 +131,18 @@ end
         @test temporal_sampling(tother) == :other
     end
 
+    @testset "sametimespan" begin
+        tm1 = Date(2000, 3, 15):Month(1):Date(2010, 3, 15)
+        tm2 = Date(2001, 1, 1):Month(1):Date(2011, 1, 1)
+        A1 = ClimArray(rand(length(tm1)), (Time(tm1),))
+        A2 = ClimArray(rand(length(tm2)), (Time(tm2),))
+        B1, B2 = sametimespan(A1, A2)
+        @test size(B1) == size(B2)
+        @test dims(B1, Time)[1] == Date(2001, 1, 15)
+        @test dims(B1, Time)[end] == Date(2010, 3, 15)
+        @test dims(B2, Time)[end] == Date(2010, 3, 1)
+    end
+
     @testset "monthlyagg and co." begin
         # test yearly temporal weights (i.e., no weighting)
         X = ClimArray(rand(3,3), (Lon(1:3), Time(tyearly[1:3])))