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