From 1d88bfb4e10f4ad2f7b49fb4cb9bc1122c3ceec8 Mon Sep 17 00:00:00 2001 From: George Datseris Date: Wed, 5 Jan 2022 15:44:43 +0100 Subject: [PATCH] Update to DimensionalData.jl (#84) * just port dd to 0.20 * add `gnv` function * LookUp -> LookupArrays * climarray can be initialized and printed * fix missing format function * fix all statistical tests * maybe fix all tests...? * correct rebuild with format * even more fixes hopefully * replace boardcast typo * change all between usages to .. * correct lookup array specification for intervalset * use `gnv` in `Coord` dimensions * split tests into files * add some basic tests for Coordinates * better naming for space types Closes #79 * fix remaining small tests --- Project.toml | 2 +- docs/src/index.md | 2 +- docs/src/plotting.md | 2 +- docs/src/statistics.md | 10 +- src/core/aggregation.jl | 16 +- src/core/coredefs.jl | 86 ++--- src/core/interpolation.jl | 9 +- src/core/nc_io.jl | 26 +- src/core/prettyprint.jl | 7 +- src/physical_dimensions/spatial.jl | 36 +- src/physical_dimensions/spatial_equalarea.jl | 102 +++--- src/physical_dimensions/temporal.jl | 5 +- src/plotting/geomakie.jl | 4 +- src/plotting/python.jl | 4 +- test/advanced_tests.jl | 19 ++ test/general_tests.jl | 56 ++++ test/io_tests.jl | 48 +++ test/runtests.jl | 331 ++----------------- test/space_coord_tests.jl | 49 +++ test/space_tests.jl | 71 ++++ test/temporal_tests.jl | 138 ++++++++ 21 files changed, 558 insertions(+), 465 deletions(-) create mode 100644 test/advanced_tests.jl create mode 100644 test/general_tests.jl create mode 100644 test/io_tests.jl create mode 100644 test/space_coord_tests.jl create mode 100644 test/space_tests.jl create mode 100644 test/temporal_tests.jl diff --git a/Project.toml b/Project.toml index cd4df5aa..c81d94ff 100644 --- a/Project.toml +++ b/Project.toml @@ -15,7 +15,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" [compat] -DimensionalData = "0.18" +DimensionalData = "0.20" Interpolations = "0.13.2" NCDatasets = "0.10, 0.11" Requires = "1" diff --git a/docs/src/index.md b/docs/src/index.md index 3badd160..68bb3de1 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -21,7 +21,7 @@ Make sure your installed version coincides with the one in this docs (see bottom ## `ClimArray`: the core data structure This project treats "climate data" as an instance of [`ClimArray`](@ref). -At the moment `ClimArray` is a subtype of `DimensionalArray` from DimensionalData.jl. +At the moment `ClimArray` is a subtype of `DimArray` from DimensionalData.jl. A brief introduction to DimensionalData.jl is copied here from its docs, because basic handling of a `ClimArray` comes from DimensionalData.jl. DimensionalData.jl allows to dimensionally-index data by their values. diff --git a/docs/src/plotting.md b/docs/src/plotting.md index 94a16509..a07f12e8 100644 --- a/docs/src/plotting.md +++ b/docs/src/plotting.md @@ -5,4 +5,4 @@ ## Python based plotting with `cartopy` For now, you can use PyCall.jl, matplotlib, and the Python library cartopy. -In the file [`ClimateBase/plotting/python.jl`](https://github.com/JuliaClimate/ClimateBase.jl/tree/master/plotting/python.jl) we provide two functions that plot maps of `ClimArray` in arbitrary projections: `earthsurface` for `LonLatGrid` and `earthscatter` for `UnstructuredGrid`. You can incorporate these in your source code as a temporary solution. +In the file [`ClimateBase/plotting/python.jl`](https://github.com/JuliaClimate/ClimateBase.jl/tree/master/plotting/python.jl) we provide two functions that plot maps of `ClimArray` in arbitrary projections: `earthsurface` for `OrthogonalSpace` and `earthscatter` for `CoordinateSpace`. You can incorporate these in your source code as a temporary solution. diff --git a/docs/src/statistics.md b/docs/src/statistics.md index d09271fd..e9bdcb69 100644 --- a/docs/src/statistics.md +++ b/docs/src/statistics.md @@ -19,7 +19,7 @@ sametimespan ## Spatial -All functions in this section work for both types of space, see [Types of spatial coordinates](@ref). +All functions in this section work for both types of space, see [Types of spatial information](@ref). ```@docs zonalmean latmean @@ -32,11 +32,11 @@ lonlatfirst longitude_circshift ``` -### Types of spatial coordinates -At the moment the following type of spatial coordinates are supported: +### Types of spatial information +Spatial information (excluding height/pressure dimensions) in ClimateBase.jl exists in one of two forms: ```@docs -LonLatGrid -UnstructuredGrid +OrthogonalSpace +CoordinateSpace ``` ClimateBase.jl works with either type of spatial coordinate system. diff --git a/src/core/aggregation.jl b/src/core/aggregation.jl index bf3c9a00..a243fd1d 100644 --- a/src/core/aggregation.jl +++ b/src/core/aggregation.jl @@ -153,7 +153,7 @@ end ######################################################################### #= """ - dimwise(f, A, B) + broadcast_dims(f, A, B) Apply bivariate function `f` (e.g. multiplication `*`) across the (first) matching dimension of `A` and `B`, where `B` has only one dimension. For normal Julia arrays match is done on the first matching dimension size. @@ -168,18 +168,18 @@ for every `i` across the matching dimension. Useful when wanting to weight a field `A` based on its latitude or time. """ -function dimwise(f, A::AbDimArray, B::AbDimArray) - @assert length(size(B)) == 1 "3rd argument of dimwise must be a vector" +function broadcast_dims(f, A::AbDimArray, B::AbDimArray) + @assert length(size(B)) == 1 "3rd argument of broadcast_dims must be a vector" m = find_matching_dim(A, B) # Array(data) necessary because reshapes somewhow fail to play well with the code... - C = dimwise(f, Array(data(A)), Array(data(B)), m) + C = broadcast_dims(f, Array(data(A)), Array(data(B)), m) return DimensionalData.basetypeof(A)(C, dims(A)) end -dimwise(f, A::AbstractArray, B, m::Int = find_matching_dim(A, B)) = -dimwise!(copy(A), f, A, B, m) +broadcast_dims(f, A::AbstractArray, B, m::Int = find_matching_dim(A, B)) = +broadcast_dims!(copy(A), f, A, B, m) -function dimwise!(C, f, A::AbstractArray{T, N}, B::AbstractVector, m::Int) where {T, N} +function broadcast_dims!(C, f, A::AbstractArray{T, N}, B::AbstractVector, m::Int) where {T, N} shape = ntuple(i -> (i==m) ? length(B) : 1, N) b = reshape(B, shape) C .= f.(A, b) @@ -200,4 +200,4 @@ function find_matching_dim(A::AbstractArray, B::AbstractVector)::Int end =# -export dimwise +export broadcast_dims diff --git a/src/core/coredefs.jl b/src/core/coredefs.jl index 5e3efe7c..ca07ee19 100644 --- a/src/core/coredefs.jl +++ b/src/core/coredefs.jl @@ -2,8 +2,8 @@ # Basic imports and dimension definitions ########################################################################################## using DimensionalData -using DimensionalData: @dim, hasdim, Dimension, IndependentDim -using DimensionalData: basetypeof +using DimensionalData.Dimensions, DimensionalData.LookupArrays +using DimensionalData: basetypeof, broadcast_dims using Dates Time = DimensionalData.Ti @@ -19,18 +19,18 @@ dimindex(A, dim) = DimensionalData.dimnum(A, dim) """ gnv(object) → x -Short for "get numeric value", this function will return the pure numeric value -of the given object. Convenience function for quickly gettin the numeric data of +Short for "get numeric value", this function will return the pure numeric value(s) +of the given object. Convenience function for quickly getting the numeric data of dimensional arrays or dimensions. """ gnv(x) = x -gnv(x::AbDimArray) = DimensionalData.parent(x) -gnv(x::Dimension) = x.val +gnv(x::Union{AbDimArray, LookupArray}) = parent(x) +gnv(x::Dimension) = parent(parent(x)) -export At, Between, Near # Selectors from DimensionalArrays.jl +export At, (..), Between, Near # Selectors from DimensionalArrays.jl export hasdim, dims, dimindex export Time, Lon, Lat, dims, Coord, Hei, Pre, Ti -export UnstructuredGrid, LonLatGrid, spacestructure +export CoordinateSpace, OrthogonalSpace, spacestructure export DimensionalData # for accessing its functions export gnv @@ -69,55 +69,53 @@ const COMMONNAMES = Dict( ########################################################################################## # the following traits for the the way space is configured. currently the options are -# UnstructuredGrid, which is for equal area "grids" (or better points/coordinates) -# while the LonLatGrid is for standard Longitude x Latitude dimensions. -# Dispatch on `UnstructuredGrid` or other types is type-stable +# CoordinateSpace, which is for points/coordinates +# while the OrthogonalSpace is for standard Longitude x Latitude dimensions. +# Dispatch on `CoordinateSpace` or other types is type-stable # because it comes from the underlying dimension abstract type SpaceType end """ -Space coordinates are represented by two orthogonal dimensions `Lon, Lat`, +Space information is represented by two orthogonal dimensions `Lon, Lat`, one being longitude and the other being latitude. """ -struct LonLatGrid <: SpaceType end +struct OrthogonalSpace <: SpaceType end """ -Space coordinates are represented by a single dimension `Coord`, whose -elements are coordinate locations, i.e. 2-element `SVector(longitude, latitude)`. -Each coordinate represents an **equal area polygon** corresponding to the point in space. +Space information is represented by a single dimension `Coord`, whose +elements are coordinates, i.e. 2-element `SVector(longitude, latitude)`. +Each coordinate represents the center of an arbitrary polygon in space. The actual limits of each polygon are not included in the dimension for performance reasons. -This dimension allows indexing according to the underlying `Lon, Lat` representation, -e.g. you can do +In statistical functions such as [`spaceagg`](@ref), it is assumed that entry of the +coordinates covers **an equal amount of area**. If this is not the case, you can simply +provide an additional weights vector which would correspond to the area covered. + +This dimension also allows indexing by latitude ranges, e.g. you can do ```julia -A # some `ClimArray` with unstructured grid type. -A[Coord(Lon(Between(0, 30)), Lat(Between(-30, 30)))] +A # some `ClimArray` with a `Coord` dimension +A[Coord(Lat(-30..30)))] ``` -To use functions such as [`zonalmean`](@ref) or [`hemispheric_means`](@ref) with this grid, -you must first sort the `ClimArray` so that the latitudes -of its coordinates are sorted in ascending order. I.e. +Most functions of ClimateBase.jl implicitly assume that the coordinates are +sorted by latidude. You can achieve this with the following code: ```julia -A # some `ClimArray` with unstructured grid type. -coords = dims(A, Coord).val -si = sortperm(coords, by = reverse) +A # some `ClimArray` with a `Coord` dimension +coords = gnv(dims(A, Coord)) +si = sortperm(coords; by = reverse) A = A[Coord(si)] ``` **This is done automatically by [`ncread`](@ref).** - -!!! warn - `UnstructuredGrid` functionality is currently in an **experimental phase**! - Notice that non-equal area unstructured grids are not supported yet. """ -struct UnstructuredGrid <: SpaceType end +struct CoordinateSpace <: SpaceType end spacestructure(a::AbDimArray) = spacestructure(dims(a)) function spacestructure(dims) if hasdim(dims, Coord) - UnstructuredGrid() + CoordinateSpace() elseif hasdim(dims, Lon) || hasdim(dims, Lat) - LonLatGrid() + OrthogonalSpace() else error("Array does not have any spatial dimensions: `Lon`, `Lat`, or `Coord`.") end @@ -126,17 +124,18 @@ end ########################################################################################## # ClimArray definition and DimensionalData.jl extensions ########################################################################################## +# TODO: remove `refdims` field export ClimArray -struct ClimArray{T,N,D<:Tuple,R<:Tuple,A<:AbstractArray{T,N},Me} <: AbstractDimensionalArray{T,N,D,A} +struct ClimArray{T,N,D<:Tuple,R<:Tuple,A<:AbstractArray{T,N},Me} <: AbDimArray{T,N,D,A} data::A dims::D refdims::R name::Symbol attrib::Me end -ClimArray(A::DimensionalArray) = ClimArray(A.data, A.dims, A.refdims, A.name, A.metadata) +ClimArray(A::AbDimArray) = ClimArray(A.data, A.dims, A.refdims, A.name, A.metadata) ClimArray(A::ClimArray, dims::Tuple = A.dims; name = A.name, attrib = A.attrib) = -ClimArray(A.data, DimensionalData.formatdims(A.data, dims), A.refdims, Symbol(name), attrib) +ClimArray(A.data, format(dims, A.data), A.refdims, Symbol(name), attrib) """ ClimArray(A::Array, dims::Tuple; name = "", attrib = nothing) @@ -146,7 +145,7 @@ information, a name and an `attrib` field (typically a dictionary) that holds ge attributes. You can think of `ClimArray` as a in-memory representation of a CFVariable. -At the moment, a `ClimArray` is using `DimensionalArray` from DimensionalData.jl, and +At the moment, a `ClimArray` is using `DimArray` from DimensionalData.jl, and all basic handling of `ClimArray` is offered by `DimensionalData` (see below). `ClimArray` is created by passing in standard array data `A` and a @@ -167,23 +166,24 @@ A = ClimArray(data, dimensions) ``` """ ClimArray(A::AbstractArray, dims::Tuple; refdims=(), name="", attrib=nothing) = - ClimArray(A, DimensionalData.formatdims(A, dims), refdims, Symbol(name), attrib) + ClimArray(A, format(dims, A), refdims, Symbol(name), attrib) ClimArray(A::AbstractArray, dims::Tuple, name; refdims=(), attrib=nothing) = - ClimArray(A, DimensionalData.formatdims(A, dims), refdims, Symbol(name), attrib) + ClimArray(A, format(dims, A), refdims, Symbol(name), attrib) Base.parent(A::ClimArray) = A.data Base.@propagate_inbounds Base.setindex!(A::ClimArray, x, I::Vararg{DimensionalData.StandardIndices}) = setindex!(A.data, x, I...) DimensionalData.metadata(A::ClimArray) = A.attrib -DimensionalData.rebuild(A::ClimArray, data::Any, dims::Tuple=dims(A), refdims=DimensionalData.refdims(A), -name="", attrib=nothing) = ClimArray(data, dims, refdims, Symbol(name), attrib) DimensionalData.basetypeof(::ClimArray) = ClimArray +DimensionalData.rebuild( + ::ClimArray, data, dims, refdims, name, metadata +) = ClimArray(data, format(dims, data), refdims, name, metadata) DimensionalData.rebuild( A::ClimArray; - data=data(A), dims=dims(A), refdims=refdims(A), name=name(A), metadata=metadata(A) -) = ClimArray(data, dims, refdims, name, metadata) + data=gnv(A), dims=dims(A), refdims=refdims(A), name=name(A), metadata=metadata(A) +) = ClimArray(data, format(dims, data), refdims, name, metadata) # The following basic methods allow indexing with tuples, (Time(5), Lon(3)) diff --git a/src/core/interpolation.jl b/src/core/interpolation.jl index 904737ba..3ea350ac 100644 --- a/src/core/interpolation.jl +++ b/src/core/interpolation.jl @@ -86,24 +86,21 @@ function interpolate_height2pressure(A::ClimArray,pressure_levels::Vector; extra error("Height is not an ascending coordinate.") end - pressure = height2pressure.(dims(A,Hei()).val) + pressure = height2pressure.(gnv(dims(A,Hei()))) # construct output Array: pre = Pre(pressure_levels; metadata = Dict()) dims_no_height = otherdims(A, Hei()) - out_dims = (dims_no_height...,pre) + out_dims = (dims_no_height..., pre) dimension_lengths = length.(out_dims) - int_array = ClimArray(zeros(eltype(A), Tuple(dimension_lengths)), out_dims ; name = A.name, attrib = A.attrib) + int_array = ClimArray(zeros(eltype(A), Tuple(dimension_lengths)), out_dims; name = A.name, attrib = A.attrib) for i in otheridxs(A, Hei()) - itp = LinearInterpolation(pressure,A[i],extrapolation_bc=extrapolation_bc) int_array[i] = itp(pressure_levels) - end return int_array - end diff --git a/src/core/nc_io.jl b/src/core/nc_io.jl index 60278459..0a068c34 100644 --- a/src/core/nc_io.jl +++ b/src/core/nc_io.jl @@ -113,8 +113,8 @@ NCDatasets.jl and then convert to some kind of dimensional container. ## Keywords * `name` optionally rename loaded array. -* `grid = nothing` optionally specify whether the underlying grid is `grid = LonLatGrid()` - or `grid = UnstructuredGrid()`, see [Types of spatial coordinates](@ref). +* `grid = nothing` optionally specify whether the underlying grid is `grid = OrthogonalSpace()` + or `grid = CoordinateSpace()`, see [Types of spatial information](@ref). If `nothing`, we try to deduce automatically based on the names of dimensions and other keys of the `NCDataset`. * `lon, lat`. These two keywords are useful in unstructured grid data where the grid @@ -126,7 +126,7 @@ NCDatasets.jl and then convert to some kind of dimensional container. lon = Array(ds["clon"]); lat = Array(ds["clat"]); ``` - If `lon, lat` are given, `grid` is automatically assumed `UnstructuredGrid()`. + If `lon, lat` are given, `grid` is automatically assumed `CoordinateSpace()`. """ function ncread(path::Union{String, Vector{String}}, args...; kwargs...) NCDataset(path) do ds @@ -147,14 +147,14 @@ function ncread(ds::NCDatasets.AbstractDataset, var, selection = nothing; name = var2name(var), grid = nothing, lon = nothing, lat = nothing, ) if lon isa Vector && lat isa Vector - gridtype = UnstructuredGrid() + gridtype = CoordinateSpace() elseif isnothing(grid) gridtype = autodetect_grid(ds) else gridtype = grid end - if gridtype == UnstructuredGrid() + if gridtype == CoordinateSpace() return ncread_unstructured(ds, var, name, lon, lat, selection) else return ncread_lonlat(ds, var, name, selection) @@ -163,11 +163,11 @@ end function autodetect_grid(ds) if haskey(ds, "reduced_points") || haskey(ds.dim, "ncells") || haskey(ds, "clon") - return UnstructuredGrid() + return CoordinateSpace() elseif haskey(ds, "lat") && length(size(ds["lat"])) > 1 - return UnstructuredGrid() + return CoordinateSpace() else - return LonLatGrid() + return OrthogonalSpace() end end @@ -175,7 +175,7 @@ var2name(var) = string(var) var2name(var::Tuple) = join(var, "_") ######################################################################### -# Reading: LonLatGrid and main reading functionality +# Reading: OrthogonalSpace and main reading functionality ######################################################################### # Notice that this function properly loads even without any spatial coordinate function ncread_lonlat(ds::NCDatasets.AbstractDataset, var, name, selection) @@ -291,7 +291,7 @@ vector2range(r::AbstractRange) = r ######################################################################### -# Reading: UnstructuredGrid +# Reading: CoordinateSpace ######################################################################### export SVector @@ -342,7 +342,7 @@ function ncread_unstructured( # Make coordinate dimension lonlat = lonlat[sel[i]] - si = sortperm(lonlat, by = reverse) + si = sortperm(lonlat; by = reverse) coords = Coord(lonlat, (Lon, Lat)) insert!(actualdims, i, coords) @@ -474,7 +474,7 @@ function ncwrite(file::String, Xs; globalattr = Dict()) # into longitude and latitude. if any(X -> hasdim(X, Coord), Xs) error(""" - Outputing `UnstructuredGrid` coordinates to .nc files is not yet supported, + Outputing `CoordinateSpace` coordinates to .nc files is not yet supported, but it is an easy fix, see source of `ncwrite`. """) end @@ -511,7 +511,7 @@ function add_dims_to_ncfile!(ds::NCDatasets.AbstractDataset, dimensions::Tuple) eltype(v) == Date && (v = DateTime.(v)) l = length(v) NCDatasets.defDim(ds, d, l) # add dimension entry - attrib = dimensions[i].metadata + attrib = DimensionalData.metadata(dimensions[i]) if (isnothing(attrib) || attrib == DimensionalData.NoMetadata()) && haskey(DEFAULT_ATTRIBS, d) @warn "Dimension $d has no attributes, adding default attributes (mandatory)." attrib = DEFAULT_ATTRIBS[d] diff --git a/src/core/prettyprint.jl b/src/core/prettyprint.jl index b661fb52..87c4f5ce 100644 --- a/src/core/prettyprint.jl +++ b/src/core/prettyprint.jl @@ -5,10 +5,10 @@ function Base.show(io::IO, ::MIME"text/plain", A::ClimArray) summary(io, A) print(io, "and") printstyled(io, " data: "; color=:green) - dataA = data(A) + dataA = gnv(A) print(io, summary(dataA), "\n") x = 2length(dims(A)) + attriblength(A.attrib) + 5 - custom_show(io, data(A), x) + custom_show(io, dataA, x) end attriblength(d::AbstractDict) = length(d) @@ -23,8 +23,7 @@ function Base.summary(io::IO, A::ClimArray) print(io, ")") end - print(io, " with dimensions:") - DimensionalData._layout_dims(io, MIME"text/plain"(), dims(A)) + DimensionalData.Dimensions.print_dims(io, MIME"text/plain"(), dims(A)) print(io, "\n") if !isnothing(A.attrib) diff --git a/src/physical_dimensions/spatial.jl b/src/physical_dimensions/spatial.jl index b679987c..28e9753b 100644 --- a/src/physical_dimensions/spatial.jl +++ b/src/physical_dimensions/spatial.jl @@ -21,7 +21,7 @@ end Works for all types of space (`...` is necessary because `i` is a `Tuple`). """ spatialidxs(A::AbDimArray) = spatialidxs(spacestructure(A), A) -function spatialidxs(::LonLatGrid, A) +function spatialidxs(::OrthogonalSpace, A) lons = (Lon(i) for i in 1:size(A, Lon)) lats = (Lat(i) for i in 1:size(A, Lat)) return Iterators.product(lons, lats) @@ -95,12 +95,12 @@ export latmean, spacemean, zonalmean, spaceagg, uniquelats """ zonalmean(A::ClimArray [, W]) -Return the zonal mean of `A`. Works for both [`LonLatGrid`](@ref) as well as -[`UnstructuredGrid`](@ref). Optionally provide statistical weights `W`. +Return the zonal mean of `A`. Works for both [`OrthogonalSpace`](@ref) as well as +[`CoordinateSpace`](@ref). Optionally provide statistical weights `W`. These can be the same `size` as `A` or only having the same latitude structure as `A`. """ zonalmean(A::AbDimArray, W = nothing) = zonalmean(spacestructure(A), A, W) -zonalmean(::LonLatGrid, A::AbDimArray, W) = dropagg(mean, A, Lon, W) +zonalmean(::OrthogonalSpace, A::AbDimArray, W) = dropagg(mean, A, Lon, W) """ latmean(A::ClimArray) @@ -110,7 +110,7 @@ This function properly weights by the cosine of the latitude. function latmean(A::AbDimArray) lw = _latweights(A) if ndims(A) > 1 - return dropagg(sum, dimwise(*, A, lw), Lat) + return dropagg(sum, broadcast_dims(*, A, lw), Lat) else return sum(A .* lw) end @@ -143,7 +143,7 @@ longitude and latitude) of `A`, weighting every part of `A` by its spatial area. `ClimArray` with same spatial information as `A`, or having exactly same dimensions as `A`. """ spaceagg(f, A::AbDimArray, exw=nothing) = spaceagg(spacestructure(A), f, A, exw) -function spaceagg(::LonLatGrid, f, A::AbDimArray, w=nothing) +function spaceagg(::OrthogonalSpace, f, A::AbDimArray, w=nothing) wtype = spaceweightassert(A, w) cosweights = repeat(cosd.(dims(A, Lat).val)', size(A, Lon)) if dimindex(A, Lon) > dimindex(A, Lat) @@ -165,8 +165,6 @@ function spaceagg(::LonLatGrid, f, A::AbDimArray, w=nothing) if wtype != :dany r = map(i -> f(view(A, i), W), oidxs) else - # TODO: This multiplication .* here assumes that the lon-lat grid is the - # first dimension. r = map(i -> f(view(A, i), weights(view(w, i) .* cosweights)), oidxs) end return ClimArray(r, other, A.name) @@ -205,9 +203,9 @@ appropriately translating the latitudes of `south` so that both arrays have the latitudinal dimension (and thus can be compared and do opperations between them). """ hemispheric_functions(A) = hemispheric_functions(spacestructure(A), A) -function hemispheric_functions(::LonLatGrid, A) - nh = A[Lat(Between(0, 90))] - sh = A[Lat(Between(-90, 0))] +function hemispheric_functions(::OrthogonalSpace, A) + nh = A[Lat(0..90)] + sh = A[Lat(-90..0)] # TODO: this can be a function "reverse dim" di = dimindex(sh, Lat) newdims = [dims(sh)...] @@ -226,20 +224,20 @@ Use [`hemispheric_functions`](@ref) to just split `A` into two hemispheres. Optionally provide weights `W` that need to have the same structure as [`spaceagg`](@ref). """ hemispheric_means(A, args...) = hemispheric_means(spacestructure(A), A, args...) -function hemispheric_means(::LonLatGrid, A::AbDimArray) +function hemispheric_means(::OrthogonalSpace, A::AbDimArray) @assert hasdim(A, Lat) if hasdim(A, Lon) B = zonalmean(A) else B = A end - nh = latmean(B[Lat(Between(0, 90))]) - sh = latmean(B[Lat(Between(-90, 0))]) + nh = latmean(B[Lat(0..90)]) + sh = latmean(B[Lat(-90..0)]) return nh, sh end latitudes(A) = latitudes(spacestructure(A), A) -latitudes(::LonLatGrid, A) = dims(A, Lat).val +latitudes(::OrthogonalSpace, A) = dims(A, Lat).val ######################################################################### # Tropics/extratropics ######################################################################### @@ -252,11 +250,11 @@ having [-higher_lat:-lower_lat, lower_lat:higher_lat]. tropics_extratropics(A, args...; kwargs...) = tropics_extratropics(spacestructure(A), A, args...; kwargs...) -function tropics_extratropics(::LonLatGrid, A; lower_lat=30, higher_lat=90) - tropics = A[Lat(Between(-lower_lat, lower_lat))] +function tropics_extratropics(::OrthogonalSpace, A; lower_lat=30, higher_lat=90) + tropics = A[Lat(-lower_lat..lower_lat)] latdim = dims(A, Lat) - extra_idxs_sh = DimensionalData.sel2indices(latdim, Between(-higher_lat, -lower_lat)) - extra_idxs_nh = DimensionalData.sel2indices(latdim, Between(lower_lat, higher_lat)) + extra_idxs_sh = DimensionalData.selectindices(latdim, -higher_lat..(-lower_lat)) + extra_idxs_nh = DimensionalData.selectindices(latdim, lower_lat..higher_lat) extra_idxs = vcat(extra_idxs_sh, extra_idxs_nh) extratropics = A[Lat(extra_idxs)] return tropics, extratropics diff --git a/src/physical_dimensions/spatial_equalarea.jl b/src/physical_dimensions/spatial_equalarea.jl index 07bb4ca5..cf3073d3 100644 --- a/src/physical_dimensions/spatial_equalarea.jl +++ b/src/physical_dimensions/spatial_equalarea.jl @@ -1,3 +1,31 @@ +######################################################################### +# Basics +######################################################################### +""" + 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. +""" +uniquelats(A::AbDimArray) = uniquelats(gnv(dims(A, Coord))) +function uniquelats(c) + @assert issorted(c; by = x -> x[2]) + idxs = Vector{UnitRange{Int}}() + lats = eltype(eltype(c))[] + iprev = 1 + for i in 2:length(c) + if c[i][2] != c[i-1][2] + push!(lats, c[i-1][2]) + push!(idxs, iprev:(i-1)) + iprev = i + end + end + push!(lats, c[end][2]) + push!(idxs, iprev:length(c)) + return idxs, lats +end + + """ transform_to_coord(A::ClimArray) → B Transform given `A` to a new `B::ClimArray` so that the `Lon, Lat` dimensions in `A` @@ -6,13 +34,11 @@ are transformed to a `Coord` dimension in `B`. function transform_to_coord(A; name = A.name, attrib = A.attrib) hasdim(A, Coord) && return A @assert hasdim(A, Lon) && hasdim(A, Lat) - # First make coord dimension londim = dims(A, Lon) latdim = dims(A, Lat) lonlat = [SVector(l, lat) for lat in latdim for l in londim] coorddim = Coord(lonlat, (Lon, Lat)) - # Then, reshape A to have this dimension is = dimindex(A, (Lon, Lat)) sizes = [size(A)...] @@ -20,30 +46,27 @@ function transform_to_coord(A; name = A.name, attrib = A.attrib) deleteat!(sizes, is[2]) B = reshape(copy(A.data), sizes...) i = is[1] - # Then, make new dimensions - newdims = [dims(A)...] + newdims = Any[dims(A)...] deleteat!(newdims, is) insert!(newdims, i, coorddim) X = ClimArray(B, Tuple(newdims)) - # Don't forget to sort lonlat coordinates si = sortperm(lonlat, by = reverse) X = X[Coord(si)] return ClimArray(X; name = Symbol(name), attrib) end - ######################################################################### # Spatial aggregation functions ######################################################################### -spaceagg(::UnstructuredGrid, f, A, W = nothing) = dropagg(f, A, Coord, W) +spaceagg(::CoordinateSpace, f, A, W = nothing) = dropagg(f, A, Coord, W) -function spatialidxs(::UnstructuredGrid, A) +function spatialidxs(::CoordinateSpace, A) return ((Coord(i),) for i in 1:size(A, Coord)) end -function zonalmean(::UnstructuredGrid, A::AbDimArray, ::Nothing) +function zonalmean(::CoordinateSpace, A::AbDimArray, ::Nothing) idxs, lats = uniquelats(A) other = otherdims(A, Coord()) r = zeros(eltype(A), (length(lats), size.(Ref(A), other)...)) @@ -55,7 +78,7 @@ function zonalmean(::UnstructuredGrid, A::AbDimArray, ::Nothing) end return R end -function zonalmean(::UnstructuredGrid, A::AbDimArray{T, 1}, ::Nothing) where {T} +function zonalmean(::CoordinateSpace, A::AbDimArray{T, 1}, ::Nothing) where {T} idxs, lats = uniquelats(A) res = zeros(T, length(lats)) for (i, r) in enumerate(idxs) @@ -65,7 +88,7 @@ function zonalmean(::UnstructuredGrid, A::AbDimArray{T, 1}, ::Nothing) where {T} end # zonal mean with weights -function zonalmean(::UnstructuredGrid, A::ClimArray, W::AbstractArray) +function zonalmean(::CoordinateSpace, A::ClimArray, W::AbstractArray) @assert size(A) == size(W) idxs, lats = uniquelats(A) other = otherdims(A, Coord()) @@ -78,7 +101,7 @@ function zonalmean(::UnstructuredGrid, A::ClimArray, W::AbstractArray) end return R end -function zonalmean(::UnstructuredGrid, A::ClimArray{T, 1}, W::AbstractArray) where {T} +function zonalmean(::CoordinateSpace, A::ClimArray{T, 1}, W::AbstractArray) where {T} @assert size(A) == size(W) idxs, lats = uniquelats(A) res = zeros(T, length(lats)) @@ -89,39 +112,15 @@ function zonalmean(::UnstructuredGrid, A::ClimArray{T, 1}, W::AbstractArray) whe end -""" - 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. -""" -uniquelats(A::AbDimArray) = uniquelats(dims(A, Coord).val) -function uniquelats(c) - @assert issorted(c; by = x -> x[2]) - idxs = Vector{UnitRange{Int}}() - lats = eltype(eltype(c))[] - iprev = 1 - for i in 2:length(c) - if c[i][2] != c[i-1][2] - push!(lats, c[i-1][2]) - push!(idxs, iprev:(i-1)) - iprev = i - end - end - push!(lats, c[end][2]) - push!(idxs, iprev:length(c)) - return idxs, lats -end - ######################################################################### # Special latitude splitting ######################################################################### -function hemispheric_functions(::UnstructuredGrid, A) - c = dims(A, Coord).val +function hemispheric_functions(::CoordinateSpace, A) + c = gnv(dims(A, Coord)) nhi, shi = hemisphere_indices(c) nh = A[Coord(nhi)] sh = A[Coord(shi)] - oldc = dims(sh, Coord).val + oldc = gnv(dims(sh, Coord)) si = sortperm(oldc, by = x -> x[2], rev = true) newc = [SVector(x[1], abs(x[2])) for x in oldc[si]] di = dimindex(sh, Coord) @@ -132,14 +131,14 @@ function hemispheric_functions(::UnstructuredGrid, A) return nh, sh end -function hemispheric_means(::UnstructuredGrid, A::AbDimArray) +function hemispheric_means(::CoordinateSpace, A::AbDimArray) nhi, shi = hemisphere_indices(A) nh = dropagg(mean, A[Coord(nhi)], Coord) sh = dropagg(mean, A[Coord(shi)], Coord) return nh, sh end -function hemispheric_means(::UnstructuredGrid, A::AbDimArray, W) +function hemispheric_means(::CoordinateSpace, A::AbDimArray, W) @assert size(A) == size(W) nhi, shi = hemisphere_indices(A) nh = dropagg(mean, A[Coord(nhi)], Coord, W[Coord(nhi)]) @@ -150,20 +149,22 @@ end """ hemisphere_indices(coords) → nhi, shi Return the indices of coordinates belonging to the north and south hemispheres. +Coordinates with latitude exactly 0 are included in _both_ hemispheres. """ function hemisphere_indices(c) idxs, lats = uniquelats(c) i = findfirst(x -> x > 0, lats) + j = findlast(x -> x < 0, lats) shi = idxs[1][1]:idxs[i-1][end] - nhi = idxs[i][1]:idxs[end][end] + nhi = idxs[j+1][1]:idxs[end][end] return nhi, shi end -latitudes(::UnstructuredGrid, A) = unique!([x[2] for x in dims(A, Coord)]) +latitudes(::CoordinateSpace, A) = unique!([x[2] for x in dims(A, Coord)]) -function tropics_extratropics(::UnstructuredGrid, A; lower_lat=30) +function tropics_extratropics(::CoordinateSpace, A; lower_lat=30) # TODO: Support `higher_lat` keyword - c = dims(A, Coord).val + c = gnv(dims(A, Coord)) idxs, lats = uniquelats(c) i1 = findlast(x -> x < -lower_lat, lats) i2 = findfirst(x -> x > lower_lat, lats) @@ -195,7 +196,14 @@ function coord_latitudes_between(c, l1, l2) end # This modifies what happens on A[Coord(Lat(Between(x,y)))] -function DimensionalData.sel2indices(c::Coord, sel::Tuple{<:Lat{ <: Between}}) +function DimensionalData.selectindices(c::Coord, sel::Tuple{<:Lat{ <: Between}}) l1, l2 = sel[1].val.val - return coord_latitudes_between(c.val, l1, l2) # this is Vector{Int} + return coord_latitudes_between(gnv(c), l1, l2) # this is Vector{Int} +end + +# This modifies what happens on A[Coord(Lat(x..y))] +function DimensionalData.selectindices(c::Coord, + sel::Tuple{<:Lat{ <: DimensionalData.LookupArrays.IntervalSets.AbstractInterval}}) + l1 = sel[1].val.left; l2 = sel.val.right + return coord_latitudes_between(gnv(c), l1, l2) # this is Vector{Int} end \ No newline at end of file diff --git a/src/physical_dimensions/temporal.jl b/src/physical_dimensions/temporal.jl index 72cd81c2..44cb4e43 100644 --- a/src/physical_dimensions/temporal.jl +++ b/src/physical_dimensions/temporal.jl @@ -33,6 +33,7 @@ temporal_sampling(t::Dimension) = temporal_sampling(t.val) function temporal_sampling(t::AbstractVector{<:TimeType}) function issame(f) + f === hour && eltype(t) <: Date && return true x0 = f(t[1]) return all(i -> f(t[i]) == x0, 2:length(t)) end @@ -228,7 +229,7 @@ function sametimespan(Xs...) mint = Date(year(mint), 1, 1) maxt = Date(year(maxt), 12, 31) end - map(X -> X[Time(Between(mint, maxt))], Xs) + map(X -> X[Time(mint..maxt)], Xs) end sametimespan(Xs::Tuple) = sametimespan(Xs...) @@ -301,7 +302,7 @@ function timeagg_monthly(f, A::AbDimArray, w) elseif w isa Vector tw .* w elseif w isa AbDimArray - _w = dimwise(*, w, ClimArray(tw, (Time(t),))) + _w = broadcast_dims(*, w, ClimArray(tw, (Time(t),))) @view _w[Time(1:mys)] end other = otherdims(A, Time) diff --git a/src/plotting/geomakie.jl b/src/plotting/geomakie.jl index b9961377..a09b8b98 100644 --- a/src/plotting/geomakie.jl +++ b/src/plotting/geomakie.jl @@ -3,7 +3,7 @@ export climscatter, climscatter!, climsurface!, climsurface, climplot """ climplot(A::ClimArray; kwargs...) → fig, ax, el, cb Main plotting function that dispatches to [`climscatter!`](@ref) if `A` has -an [`UnstructuredGrid`](@ref) dimension, or to [`climsurface!`](@ref) for [`LonLatGrid`](@ref). +an [`CoordinateSpace`](@ref) dimension, or to [`climsurface!`](@ref) for [`OrthogonalSpace`](@ref). Return the figure, axis, plotted element, and colorbar. @@ -13,7 +13,7 @@ Plotting from ClimateBase.jl also works with `Observable`s that enclose a `ClimA You can update the values of the observable with another `ClimArray` with the same spatial dimension, and the plot will be updated. See documentation online for examples. """ -function climplot(A, args...; scatter = spacestructure(A) == UnstructuredGrid(), +function climplot(A, args...; scatter = spacestructure(A) == CoordinateSpace(), source = "+proj=longlat +datum=WGS84", dest = "+proj=eqearth", colorbar = true, name = string(DimensionalData.name(A)), kwargs...) diff --git a/src/plotting/python.jl b/src/plotting/python.jl index 3542546f..2f2659f2 100644 --- a/src/plotting/python.jl +++ b/src/plotting/python.jl @@ -10,7 +10,7 @@ DEFPROJ = ccrs.Mollweide() """ earthscatter(A::ClimArray, projection = ccrs.Mollweide(); kwargs...) → fig, ax, cb -Plot a `ClimArray` with space type `UnstructuredGrid` as a scatter plot with the color +Plot a `ClimArray` with space type `CoordinateSpace` as a scatter plot with the color of the points being the value of `A` at these points. This requires that `A` has only one dimension, the coordinates. @@ -62,7 +62,7 @@ end """ earthsurface(A::ClimArray, projection = ccrs.Mollweide(); kwargs...) → fig, ax, cb -Plot a `ClimArray` with space type `LonLatGrid` as a surface plot with the color +Plot a `ClimArray` with space type `OrthogonalSpace` as a surface plot with the color of the surface being the value of `A`. This requires that `A` has exactly two dimensions, `Lon, Lat`. diff --git a/test/advanced_tests.jl b/test/advanced_tests.jl new file mode 100644 index 00000000..9759986f --- /dev/null +++ b/test/advanced_tests.jl @@ -0,0 +1,19 @@ +@testset "Vertical interpolation" begin + D = ClimArray([1.0:1.0:11.0 2.0:1.0:12.0 3.0:1.0:13.0], (Hei(0.0:2000.0:20000.0), Ti(1:3))) + pressure_levels = [950.0,850.0,650.0,350.0,250.0,150.0] .* 100.0 + D_pre = interpolate_height2pressure(D, pressure_levels,extrapolation_bc=NaN) + D_back = interpolate_pressure2height(D_pre, Vector(dims(D,Hei).val),extrapolation_bc=Line()) + my_dim = Dim{:My_Dim} + E = ClimArray([1.0:1.0:10.0 2.0:1.0:11.0 3.0:2.0:21.0], (my_dim(1:10), Ti(1:3))) + pressure = ClimArray([1000.0:-100.0:100.0 1000.0:-100.0:100.0 1000.0:-100.0:100.0] *100.0, (my_dim(1:10), Ti(1:3))) + E_pre = interpolation2pressure(E, pressure, pressure_levels; vertical_coord=my_dim, extrapolation_bc=NaN ) + E_pre2 = interpolation2pressure(reverse(E,dims=my_dim), reverse(pressure,dims=my_dim), pressure_levels; vertical_coord=my_dim, extrapolation_bc=NaN, descending = false ) + + @test hasdim(D_pre,Pre()) + @test gnv(dims(D_back,Hei)) == gnv(dims(D,Hei)) && gnv(dims(D_back,Ti)) == gnv(dims(D,Ti)) + @test E_pre.data == E_pre2.data + @test E_pre[Pre(1)] < E_pre[Pre(2)] < E_pre[Pre(3)] < E_pre[Pre(4)] < E_pre[Pre(5)] + @test D_pre[Pre(1)] < D_pre[Pre(2)] < D_pre[Pre(3)] < D_pre[Pre(4)] < D_pre[Pre(5)] + @test D_back[Hei(1)] < D_back[Hei(2)] < D_back[Hei(3)] < D_back[Hei(4)] < D_back[Hei(5)] + @test D[Hei(1)] < D_pre[Pre(3)] < D[Hei(11)] +end diff --git a/test/general_tests.jl b/test/general_tests.jl new file mode 100644 index 00000000..76bcf07e --- /dev/null +++ b/test/general_tests.jl @@ -0,0 +1,56 @@ +@testset "pretty printing" begin + P = sprint(show, MIME"text/plain"(), A) + @test contains(P, "ClimArray") + @test contains(P, "Lon") + @test contains(P, "data") + @test contains(P, "44") +end + +@testset "General Statistics" begin +@testset "Dropping dimensions" begin + dt = timemean(A) + @test dt isa ClimArray + @test !hasdim(dt, Time) + @test hasdim(dt, Lon) + @test hasdim(dt, Lat) + dl = zonalmean(A) + @test dl isa ClimArray + @test !hasdim(dl, Lon) + @test hasdim(dl, Time) + @test hasdim(dl, Lat) + da = latmean(A) + @test da isa ClimArray + @test !hasdim(da, Lat) + @test hasdim(da, Lon) + @test hasdim(da, Time) + ds = spacemean(A) + @test ds isa ClimArray + @test !hasdim(ds, Lon) + @test !hasdim(ds, Lat) + @test hasdim(ds, Time) + ds = dropagg(mean, A, Time) + @test dt isa ClimArray + @test !hasdim(dt, Time) + @test hasdim(dt, Lon) + @test hasdim(dt, Lat) +end + +@testset "dropagg with weighting" begin + # spatiotemporal weights: + W = zero(A) + W[Time(5)] .= 1.0 + res = dropagg(mean, A, Time, W) + @test all(res .≈ A[Time(5)]) + @test dims(A, Lon).val == dims(res, Lon).val + res = dropagg(std, A, Time, W) + @test all(x -> isapprox(x, 0; atol = 1e-8), res) + + # just time weight + w = zeros(length(t)) + w[5] = 1.0 + res = dropagg(mean, A, Time, w) + @test all(res .≈ A[Time(5)]) + @test dims(A, Lon).val == dims(res, Lon).val +end + +end # general statistics \ No newline at end of file diff --git a/test/io_tests.jl b/test/io_tests.jl new file mode 100644 index 00000000..83243c82 --- /dev/null +++ b/test/io_tests.jl @@ -0,0 +1,48 @@ +@testset "NetCDF IO" begin + +@testset "NetCDF standard tests" begin + globat = Dict("history" => "test") + ncwrite("test.nc", (A, B); globalattr = globat) + Aloaded = ncread("test.nc", "insolation") + Bloaded = ncread("test.nc", "x2") + + @test A.data == Aloaded.data + @test DimensionalData.metadata(dims(Aloaded, Lon))["units"] == "degrees_east" + @test B.data == Bloaded.data + @test string(Bloaded.name) == "x2" + @test DimensionalData.metadata(dims(Bloaded, Time))["standard_name"] == "time" + + Bpartly = ncread("test.nc", "x2", ([1,4,5], :, 1:3)) + @test Bpartly.data[:, :, :] == Bloaded.data[[1,4,5], :, 1:3] + + ds = NCDataset("test.nc") + @test ds.attrib["history"] == "test" + close(ds) + rm("test.nc") +end + +@testset "Missings handling" begin + m = rand(Float64, length.((lons, lats))) + midx = CartesianIndices(m)[1:2, :] + L = length(midx) + m[midx] .= -99.9 + M = ClimArray(m, (Lon(lons), Lat(lats)); attrib = Dict("_FillValue" => -99.9), name = "M") + ncwrite("missing_test.nc", M) + M2 = ncread("missing_test.nc", "M") + for i in midx; @test ismissing(M2[i]); end + # now test `missing_weights` + B, W = missing_weights(M2) + @test all(iszero, W.data[midx]) + @test all(isequal(-99.9), B.data[midx]) + mmean = mean(skipmissing(M2.data)) + bmean = mean(B, StatsBase.weights(W)) + @test bmean == mmean + # test `missing_weights` application to zonal mean + C = B[3:end, :] + z1 = zonalmean(C) + z2 = zonalmean(B, W) + @test all(z1.data .≈ z2.data) + rm("missing_test.nc") +end + +end # NetCDF tests \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index ff42cd38..97d1f0f9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,11 +1,11 @@ using ClimateBase, Test, Dates using Statistics -using Interpolations +using ClimateBase.Interpolations import StatsBase Time = ClimateBase.Time cd(@__DIR__) -# Create the artificial dimensional array A that will be used in tests +# Create the artificial dimensional arrays function monthly_insolation(t::TimeType, args...) d = ClimateBase.monthspan(t) mean(insolation(τ, args...) for τ in d) @@ -34,316 +34,25 @@ end A = ClimArray(A, d; name = "insolation") B = ClimArray(B, d; attrib = Dict("a" => 2)) # on purpose without name -# %% General Tests -@testset "Dropping dimensions" begin - dt = timemean(A) - @test dt isa ClimArray - @test !hasdim(dt, Time) - @test hasdim(dt, Lon) - @test hasdim(dt, Lat) - dl = zonalmean(A) - @test dl isa ClimArray - @test !hasdim(dl, Lon) - @test hasdim(dl, Time) - @test hasdim(dl, Lat) - da = latmean(A) - @test da isa ClimArray - @test !hasdim(da, Lat) - @test hasdim(da, Lon) - @test hasdim(da, Time) - ds = spacemean(A) - @test ds isa ClimArray - @test !hasdim(ds, Lon) - @test !hasdim(ds, Lat) - @test hasdim(ds, Time) - ds = dropagg(mean, A, Time) - @test dt isa ClimArray - @test !hasdim(dt, Time) - @test hasdim(dt, Lon) - @test hasdim(dt, Lat) -end - -@testset "dropagg with weighting" begin - # spatiotemporal weights: - W = zero(A) - W[Time(5)] .= 1.0 - res = dropagg(mean, A, Time, W) - @test all(res .≈ A[Time(5)]) - @test dims(A, Lon).val == dims(res, Lon).val - res = dropagg(std, A, Time, W) - @test all(x -> isapprox(x, 0; atol = 1e-8), res) - - # just time weight - w = zeros(length(t)) - w[5] = 1.0 - res = dropagg(mean, A, Time, w) - @test all(res .≈ A[Time(5)]) - @test dims(A, Lon).val == dims(res, Lon).val -end - - -# %% Time Tests -@testset "Temporal weighting" begin - # spatiotemporal weights: - W = zero(A) - W[Time(5)] .= 1.0 - res = timeagg(mean, A, W) - @test all(res .≈ A[Time(5)]) - @test dims(A, Lon).val == dims(res, Lon).val - res = timeagg(std, A, W) - @test all(x -> isapprox(x, 0; atol = 1e-8), res) +# Create similar arrays with Coord dimension +reduced_points = [round(Int, length(lons)*cosd(θ)) for θ in lats] +coords = ClimateBase.reduced_grid_to_points(lats, reduced_points) +sort!(coords; by = reverse) # critical! - # just time weight - w = zeros(length(t)) - w[5] = 1.0 - res = timeagg(mean, A, w) - @test all(res .≈ A[Time(5)]) - @test dims(A, Lon).val == dims(res, Lon).val - - # TODO: more tests needed here, e.g. for timeagg(mean, t, a, w) +coord_dim = Coord(coords, (Lon, Lat)) +C = zeros(length(coords), length(t)) +for (i, (lon, θ)) in enumerate(coords) + C[i, :] .= monthly_insolation.(t, θ) .* cosd(lon) end -@testset "Advanced temporal manipulation" begin - tyearly = Date(2000, 3, 1):Year(1):Date(2020, 3, 31) - thourly = DateTime(2000, 3, 1):Hour(1):DateTime(2001, 4, 15) - mdates = unique!([(year(d), month(d)) for d in tdaily]) - ydates = unique!([year(d) for d in tdaily]) - tranges = temporalrange(tdaily, Dates.month) - yranges = temporalrange(tdaily, Dates.year) - @testset "time sampling" begin - @test length(tranges) == length(mdates) - @test length(yranges) == length(ydates) - - @test temporal_sampling(t) == :monthly - @test temporal_sampling(tdaily) == :daily - @test temporal_sampling(tyearly) == :yearly - - # Test hourly stuff - @test temporal_sampling(thourly) == :hourly - @test temporal_sampling(collect(thourly)) == :hourly - hmys = maxyearspan(thourly, :hourly) - @test hmys < length(thourly) - @test year(thourly[1]) == year(thourly[hmys])-1 - - # Arbitrary random range - tother = collect(thourly) - tother[4] = DateTime(2000, 7, 2, 2, 2) - @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]))) - W = [0, 1, 0] - @test vec(timemean(X).data) == vec(mean(X.data; dims = 2)) - @test timemean(X, W).data == X.data[:, 2] - - C = ClimArray(zeros(length(lats), length(tdaily)), (Lat(lats), Time(tdaily))) - - # First version: just count number of days - for j in 1:length(tdaily) - for i in 1:length(lats) - C[i, j] = daysinmonth(tdaily[j]) - end - end - Cm = monthlyagg(C) - - @test length(unique(Array(Cm))) == 4 # there are four unique number of days - # test that each value when rounded to an integer is an integer (for first slice only - # all remaining slices are the same) - for e in Cm[:, 1] - @test round(Int, e) == e - end - @test step(dims(Cm, Time).val) == Month(1) - @test temporal_sampling(dims(Cm, Time).val) == :monthly - - for j in 1:length(tdaily) - for i in 1:length(lats) - C[i, j] = daysinyear(tdaily[j]) - end - end - Cy = yearlyagg(C) - @test length(unique(Array(Cy))) == 2 - for e in Cy - @test round(Int, e) == e - end - @test step(dims(Cy, Time).val) == Year(1) - @test temporal_sampling(dims(Cy, Time).val) == :yearly - - # Second version: test actual physics - for j in 1:length(tdaily) - for i in 1:length(lats) - C[i, j] = insolation(tdaily[j], lats[i]) - end - end - Bz = zonalmean(B) - Cm = monthlyagg(C) - @test size(Cm, Time) == size(B, Time) + 1 - # For times besides the first and last month, they sould be the same: - for k in 5:length(t)-5 - @test all(Cm[Time(k)] .≈ Bz[Time(k)]) - end - - Asea = seasonalyagg(A) - tsea = dims(Asea, Time).val - @test Base.step(tsea) == Month(3) - end - - @testset "interannual variability" begin - C = spacemean(B) - for y in unique(year.(t))[1:end-1] - C[Time(At(Date(y, 3, 15)))] = 0 - end - dates, vals = seasonality(C) - @test length(dates) == 12 - @test all(v -> length(v) == 20, vals) - @test month(dates[3]) == 3 - @test all(iszero, vals[3]) - @test all(!iszero, vals[2]) - end -end +C = ClimArray(C, (coord_dim, Time(t)); name = "has_coords", attrib = Dict("a" => 2)) -# %% Space tests -@testset "Default physical weights" begin - x, y = hemispheric_means(B) - @test timemean(x) != dropagg(mean, x, Time) - d1 = timemean(x) - timemean(y) - @test abs(d1) < 0.01 - d2 = dropagg(mean, x, Time) - dropagg(mean, y, Time) - @test d2 < -0.1 - x = B[:, :, 1] - @test zonalmean(latmean(x)) ≈ latmean(zonalmean(x)) ≈ spacemean(x) - @test spacemean(x) - mean(x) > 50 - @test all(y -> abs(y) < 1e-8, spacemean(A)) -end - -@testset "Spatial weighting" begin - W = zero(A) - W[Lon(5)] .= 1.0 - res = spaceagg(mean, A, W) - @test all(res .≈ latmean(A)[Lon(5)]) - @test dims(A, Time) == dims(res, Time) - res = spaceagg(std, A, W) - @test all(x -> x>0, res) - - W = zeros(length(lons), length(lats)) - W[5, :] .= 1.0 - W = ClimArray(W, dims(A, (Lon, Lat))) - res = spaceagg(mean, A, W) - @test all(res .≈ latmean(A)[Lon(5)]) - @test dims(A, Time) == dims(res, Time) -end - -@testset "Insolation/Hemispheres" begin - @test !any(x -> x < 0, B) - x = timemean(B) - x, y = hemispheric_functions(zonalmean(x)) - # @test x isa ClimArray - x[Lat(1)] == y[Lat(1)] - for j in 2:size(x, Lat) - @test x[j] < x[j-1] - @test y[j] < y[j-1] - end -end - -@testset "tropics/extratropics" begin - tropics, extratropics = tropics_extratropics(B) - tlats = dims(tropics, Lat).val - elats = dims(extratropics, Lat).val - @test all(l -> -30 ≤ l ≤ 30, tlats) - @test all(l -> abs(l) ≥ 30, elats) - @test hasdim(tropics, Ti) - - tmean = spacemean(timemean(tropics)) - emean = spacemean(timemean(extratropics)) - amean = spacemean(timemean(B)) - @test amean ≈ (emean + tmean)/2 rtol = 5 - - # Test for keyword options for latitude bounds - tropics, extratropics = tropics_extratropics(B, lower_lat=20, higher_lat=78) - tlats = dims(tropics, Lat).val - elats = dims(extratropics, Lat).val - @test all(l -> -20 ≤ l ≤ 20, tlats) - @test all(l -> abs(l) ≥ 20, elats) - @test maximum(abs.(elats)) <= 78 - @test hasdim(tropics, Ti) # probably redundant, already tested - -end - -# %% IO tests -@testset "NetCDF file IO" begin - globat = Dict("history" => "test") - ncwrite("test.nc", (A, B); globalattr = globat) - Aloaded = ncread("test.nc", "insolation") - Bloaded = ncread("test.nc", "x2") - - @test A.data == Aloaded.data - @test dims(Aloaded, Lon).metadata["units"] == "degrees_east" - @test B.data == Bloaded.data - @test string(Bloaded.name) == "x2" - @test dims(Bloaded, Time).metadata["standard_name"] == "time" - - Bpartly = ncread("test.nc", "x2", ([1,4,5], :, 1:3)) - @test Bpartly.data[:, :, :] == Bloaded.data[[1,4,5], :, 1:3] - - ds = NCDataset("test.nc") - @test ds.attrib["history"] == "test" - close(ds) - rm("test.nc") -end - -@testset "Vertical interpolation" begin - D = ClimArray([1.0:1.0:11.0 2.0:1.0:12.0 3.0:1.0:13.0], (Hei(0.0:2000.0:20000.0), Ti(1:3))) - pressure_levels = [950.0,850.0,650.0,350.0,250.0,150.0] .* 100.0 - D_pre = interpolate_height2pressure(D, pressure_levels,extrapolation_bc=NaN) - D_back = interpolate_pressure2height(D_pre, Vector(dims(D,Hei).val),extrapolation_bc=Line()) - my_dim = Dim{:My_Dim} - E = ClimArray([1.0:1.0:10.0 2.0:1.0:11.0 3.0:2.0:21.0], (my_dim(1:10), Ti(1:3))) - pressure = ClimArray([1000.0:-100.0:100.0 1000.0:-100.0:100.0 1000.0:-100.0:100.0] *100.0, (my_dim(1:10), Ti(1:3))) - E_pre = interpolation2pressure(E, pressure, pressure_levels; vertical_coord=my_dim, extrapolation_bc=NaN ) - E_pre2 = interpolation2pressure(reverse(E,dims=my_dim), reverse(pressure,dims=my_dim), pressure_levels; vertical_coord=my_dim, extrapolation_bc=NaN, descending = false ) - - @test hasdim(D_pre,Pre()) - @test dims(D_back,Hei).val == dims(D,Hei).val && dims(D_back,Ti).val == dims(D,Ti).val - @test E_pre.data == E_pre2.data - @test E_pre[Pre(1)] < E_pre[Pre(2)] < E_pre[Pre(3)] < E_pre[Pre(4)] < E_pre[Pre(5)] - @test D_pre[Pre(1)] < D_pre[Pre(2)] < D_pre[Pre(3)] < D_pre[Pre(4)] < D_pre[Pre(5)] - @test D_back[Hei(1)] < D_back[Hei(2)] < D_back[Hei(3)] < D_back[Hei(4)] < D_back[Hei(5)] - @test D[Hei(1)] < D_pre[Pre(3)] < D[Hei(11)] -end - -@testset "Missings handling" begin - m = rand(Float64, length.((lons, lats))) - midx = CartesianIndices(m)[1:2, :] - L = length(midx) - m[midx] .= -99.9 - M = ClimArray(m, (Lon(lons), Lat(lats)); attrib = Dict("_FillValue" => -99.9), name = "M") - ncwrite("missing_test.nc", M) - M2 = ncread("missing_test.nc", "M") - for i in midx; @test ismissing(M2[i]); end - # now test `missing_weights` - B, W = missing_weights(M2) - @test all(iszero, W.data[midx]) - @test all(isequal(-99.9), B.data[midx]) - mmean = mean(skipmissing(M2.data)) - bmean = mean(B, StatsBase.weights(W)) - @test bmean == mmean - # test `missing_weights` application to zonal mean - C = B[3:end, :] - z1 = zonalmean(C) - z2 = zonalmean(B, W) - @test all(z1.data .≈ z2.data) - rm("missing_test.nc") -end +# %% +@testset "ClimateBase.jl tests" begin +include("general_tests.jl") +include("temporal_tests.jl") +include("space_tests.jl") +include("space_coord_tests.jl") +include("io_tests.jl") +include("advanced_tests.jl") +end \ No newline at end of file diff --git a/test/space_coord_tests.jl b/test/space_coord_tests.jl new file mode 100644 index 00000000..4abf61d6 --- /dev/null +++ b/test/space_coord_tests.jl @@ -0,0 +1,49 @@ +@testset "Spatial Coordinates Tests" begin + @testset "Indexing and Basics" begin + # Needs https://github.com/rafaqz/DimensionalData.jl/issues/357 + # subsel = C[Lat(0..8)] + # @test hasdim(subsel, Coord) + # @test size(subsel, Coord) < size(C, Coord) + # @test all(0 ≤ c[2] ≤ 8 for c in gnv(dims(subsel, Coord))) + + @test length(collect(spatialidxs(C))) == size(C, Coord) + + Atran = transform_to_coord(A) + @test hasdim(Atran, Coord) + @test lons == sort!(unique!([c[1] for c in gnv(dims(Atran, Coord))])) + + end + + @testset "Zonal mean" begin + zC = zonalmean(C) + @test hasdim(zC, Lat) + @test gnv(dims(zC, Lat)) == lats + @test maximum(zC) < 1 # because we aggregate over cosine Lon, the result is small + # TODO: test with weights + end + + @testset "Hemispheric stuff" begin + # TODO: Would be useful to test this hemispheric stuff with coordinates + # that include 0 latitude and that do not cover symmetric latitdues + # even though it is super uncommon + nhi, shi = ClimateBase.hemisphere_indices(gnv(dims(C, Coord))) + @test length(nhi) == length(shi) + @test sort!(vcat(nhi, shi)) == 1:size(C, Coord) + + Cnh, Csh = hemispheric_functions(C) # symmetric, non-zero latitudes + @test hasdim(Cnh, Coord) && hasdim(Csh, Coord) + @test size(Cnh, Coord) == size(Csh, Coord) + + # hemispheric_functions makes latitudes the same by convention: + @test gnv(dims(Cnh, Coord)) == gnv(dims(Csh, Coord)) + @test minimum(uniquelats(Cnh)[2]) ≥ 0 + @test minimum(uniquelats(Csh)[2]) ≥ 0 + + Cnh, Csh = hemispheric_means(C) + @test !hasdim(Cnh, Coord) && !hasdim(Csh, Coord) + @test hasdim(Cnh, Time) + # TODO: test actual numeric values + end + + +end \ No newline at end of file diff --git a/test/space_tests.jl b/test/space_tests.jl new file mode 100644 index 00000000..f855ef6b --- /dev/null +++ b/test/space_tests.jl @@ -0,0 +1,71 @@ +@testset "Space Tests" begin + +@testset "Default physical weights" begin + x, y = hemispheric_means(B) + @test timemean(x) != dropagg(mean, x, Time) + d1 = timemean(x) - timemean(y) + @test abs(d1) < 0.01 + d2 = dropagg(mean, x, Time) - dropagg(mean, y, Time) + @test d2 < -0.1 + x = B[:, :, 1] + @test zonalmean(latmean(x)) ≈ latmean(zonalmean(x)) ≈ spacemean(x) + @test spacemean(x) - mean(x) > 50 + @test all(y -> abs(y) < 1e-8, spacemean(A)) +end + +@testset "Spatial weighting" begin + W = zero(A) + W[Lon(5)] .= 1.0 + res = spaceagg(mean, A, W) + @test all(res .≈ latmean(A)[Lon(5)]) + @test dims(A, Time) == dims(res, Time) + res = spaceagg(std, A, W) + @test all(x -> x>0, res) + + W = zeros(length(lons), length(lats)) + W[5, :] .= 1.0 + W = ClimArray(W, dims(A, (Lon, Lat))) + res = spaceagg(mean, A, W) + @test all(res .≈ latmean(A)[Lon(5)]) + @test dims(A, Time) == dims(res, Time) +end + +@testset "Insolation/Hemispheres" begin + @test !any(x -> x < 0, B) + x = timemean(B) + x, y = hemispheric_functions(zonalmean(x)) + # @test x isa ClimArray + x[Lat(1)] == y[Lat(1)] + for j in 2:size(x, Lat) + @test x[j] < x[j-1] + @test y[j] < y[j-1] + end +end + +@testset "tropics/extratropics" begin + tropics, extratropics = tropics_extratropics(B) + tlats = dims(tropics, Lat).val + elats = dims(extratropics, Lat).val + @test all(l -> -30 ≤ l ≤ 30, tlats) + @test all(l -> abs(l) ≥ 30, elats) + @test hasdim(tropics, Ti) + + tmean = spacemean(timemean(tropics)) + emean = spacemean(timemean(extratropics)) + amean = spacemean(timemean(B)) + @test amean ≈ (emean + tmean)/2 rtol = 5 + + # Test for keyword options for latitude bounds + tropics, extratropics = tropics_extratropics(B, lower_lat=20, higher_lat=78) + tlats = dims(tropics, Lat).val + elats = dims(extratropics, Lat).val + @test all(l -> -20 ≤ l ≤ 20, tlats) + @test all(l -> abs(l) ≥ 20, elats) + @test maximum(abs.(elats)) <= 78 + @test hasdim(tropics, Ti) # probably redundant, already tested + +end + + + +end # Space Tests \ No newline at end of file diff --git a/test/temporal_tests.jl b/test/temporal_tests.jl new file mode 100644 index 00000000..29fe04ec --- /dev/null +++ b/test/temporal_tests.jl @@ -0,0 +1,138 @@ +# %% Time Tests +@testset "Temporal Tests" begin + +@testset "Temporal weighting" begin + # spatiotemporal weights: + W = zero(A) + W[Time(5)] .= 1.0 + res = timeagg(mean, A, W) + @test all(res .≈ A[Time(5)]) + @test dims(A, Lon).val == dims(res, Lon).val + res = timeagg(std, A, W) + @test all(x -> isapprox(x, 0; atol = 1e-8), res) + + # just time weight + w = zeros(length(t)) + w[5] = 1.0 + res = timeagg(mean, A, w) + @test all(res .≈ A[Time(5)]) + @test dims(A, Lon).val == dims(res, Lon).val + + # TODO: more tests needed here, e.g. for timeagg(mean, t, a, w) +end + +@testset "Advanced temporal manipulation" begin + tyearly = Date(2000, 3, 1):Year(1):Date(2020, 3, 31) + thourly = DateTime(2000, 3, 1):Hour(1):DateTime(2001, 4, 15) + mdates = unique!([(year(d), month(d)) for d in tdaily]) + ydates = unique!([year(d) for d in tdaily]) + tranges = temporalrange(tdaily, Dates.month) + yranges = temporalrange(tdaily, Dates.year) + @testset "time sampling" begin + @test length(tranges) == length(mdates) + @test length(yranges) == length(ydates) + + @test temporal_sampling(t) == :monthly + @test temporal_sampling(tdaily) == :daily + @test temporal_sampling(tyearly) == :yearly + + # Test hourly stuff + @test temporal_sampling(thourly) == :hourly + @test temporal_sampling(collect(thourly)) == :hourly + hmys = maxyearspan(thourly, :hourly) + @test hmys < length(thourly) + @test year(thourly[1]) == year(thourly[hmys])-1 + + # Arbitrary random range + tother = collect(thourly) + tother[4] = DateTime(2000, 7, 2, 2, 2) + @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]))) + W = [0, 1, 0] + @test vec(timemean(X).data) == vec(mean(X.data; dims = 2)) + @test timemean(X, W).data == X.data[:, 2] + + C = ClimArray(zeros(length(lats), length(tdaily)), (Lat(lats), Time(tdaily))) + + # First version: just count number of days + for j in 1:length(tdaily) + for i in 1:length(lats) + C[i, j] = daysinmonth(tdaily[j]) + end + end + Cm = monthlyagg(C) + + @test length(unique(Array(Cm))) == 4 # there are four unique number of days + # test that each value when rounded to an integer is an integer (for first slice only + # all remaining slices are the same) + for e in Cm[:, 1] + @test round(Int, e) == e + end + @test step(dims(Cm, Time).val) == Month(1) + @test temporal_sampling(dims(Cm, Time).val) == :monthly + + for j in 1:length(tdaily) + for i in 1:length(lats) + C[i, j] = daysinyear(tdaily[j]) + end + end + Cy = yearlyagg(C) + @test length(unique(Array(Cy))) == 2 + for e in Cy + @test round(Int, e) == e + end + @test step(dims(Cy, Time).val) == Year(1) + @test temporal_sampling(dims(Cy, Time).val) == :yearly + + # Second version: test actual physics + for j in 1:length(tdaily) + for i in 1:length(lats) + C[i, j] = insolation(tdaily[j], lats[i]) + end + end + Bz = zonalmean(B) + Cm = monthlyagg(C) + @test size(Cm, Time) == size(B, Time) + 1 + # For times besides the first and last month, they sould be the same: + for k in 5:length(t)-5 + @test all(Cm[Time(k)] .≈ Bz[Time(k)]) + end + + Asea = seasonalyagg(A) + tsea = dims(Asea, Time).val + @test Base.step(tsea) == Month(3) + end + + @testset "interannual variability" begin + C = spacemean(B) + for y in unique(year.(t))[1:end-1] + C[Time(At(Date(y, 3, 15)))] = 0 + end + dates, vals = seasonality(C) + @test length(dates) == 12 + @test all(v -> length(v) == 20, vals) + @test month(dates[3]) == 3 + @test all(iszero, vals[3]) + @test all(!iszero, vals[2]) + end +end + + + +end # tempooral tests \ No newline at end of file