Merge pull request #12 from jbisits/joey-consistent-symbols

Homogenise variables and methods, add density computations from `RasterStacks` and `RasterSeries`.

src/OceanRasterConversions.jl

@@ -2,7 +2,10 @@ module OceanRasterConversions
 
 using Rasters, GibbsSeaWater
 
-export convert_ocean_vars, depth_to_pressure, Sₚ_to_Sₐ, θ_to_Θ
+export
+    convert_ocean_vars,
+    depth_to_pressure, Sₚ_to_Sₐ, θ_to_Θ,
+    in_situ_density, potential_density
 
 include("oceanconversions.jl")
 

src/oceanconversions.jl

@@ -13,16 +13,16 @@ longitude, latitude, depth and time have a variable for pressure. A density vari
 computed which, by default, is _in-situ_ density. Potential density at a reference pressure
 can be computed instead by passing a the keyword argument `ref_pressure`.
 
-The name of the variables for potential temperature and salinity
-(either practical or absolute) must be passed in as a `NamedTuple` of the form
-`(sp = :salt_name, pt = :potential_temp_name)` where `:potential_temp_name` and `:salt_name`
-are the name of the potential temperature and salinity in the `Raster`.
+The name of the variables for potential temperature and practical salinity must be passed in
+as a `NamedTuple` of the form `(Sₚ = :salt_name, θ = :potential_temp_name)` where
+`:potential_temp_name` and `:salt_name` are the name of the potential temperature and
+practical salinity in the `Raster`.
 """
 function convert_ocean_vars(stack::RasterStack, var_names::NamedTuple;
                             ref_pressure = nothing)
 
-    Sₚ = read(stack[var_names.sp])
-    θ = read(stack[var_names.pt])
+    Sₚ = read(stack[var_names.Sₚ])
+    θ = read(stack[var_names.θ])
     rs_dims = get_dims(Sₚ)
     p = depth_to_pressure(Sₚ, rs_dims)
     find_nm = @. !ismissing(Sₚ) && !ismissing(θ)
@@ -73,17 +73,21 @@ function depth_to_pressure(raster::Raster, rs_dims::Tuple)
     return Raster(p, rs_dims)
 
 end
+"""
+    function depth_to_pressure(stack::RasterStack)
+Convert the depth dimension `Z` from a `RasterStack`.
+"""
 depth_to_pressure(stack::RasterStack) = depth_to_pressure(stack[keys(stack)[1]],
                                                           get_dims(stack[keys(stack)[1]]))
 """
     function Sₚ_to_Sₐ(raster::Raster, p::raster, rs_dims::Tuple, find_nm::Raster)
 Convert a `Raster` of practical salinity (`Sₚ`) to absolute salinity (`Sₐ`) using
-`gsw_sa_from_sp` from GibbsSeaWater.jl.
+`gsw_sa_from_sp` from GibbsSeaWater.jl. This conversion depends on pressure.
 """
 function Sₚ_to_Sₐ(Sₚ::Raster, p::Raster, rs_dims::Tuple, find_nm::Raster)
 
     lons, lats, z, time = rs_dims
-    Sₐ = similar(Array(Sₚ), Union{Float64, Missing})
+    Sₐ = similar(Array(Sₚ))
 
     if isnothing(time)
 
@@ -109,11 +113,18 @@ function Sₚ_to_Sₐ(Sₚ::Raster, p::Raster, rs_dims::Tuple, find_nm::Raster)
     return Raster(Sₐ, rs_dims)
 
 end
-Sₚ_to_Sₐ(stack::RasterStack, sp::Symbol) = Sₚ_to_Sₐ(stack[sp],
+"""
+    function Sₚ_to_Sₐ(stack::RasterStack, Sₚ::Symbol)
+    function Sₚ_to_Sₐ(series::RasterSeries, Sₚ::Symbol)
+Convert only the practical salinity, returning the absolute salinity, from a `RasterStack`
+or `RasterSeries`. The symbol for the practical salinity in the `RasterStack/Series` must be
+passed in.
+"""
+Sₚ_to_Sₐ(stack::RasterStack, Sₚ::Symbol) = Sₚ_to_Sₐ(stack[Sₚ],
                                                     depth_to_pressure(stack),
-                                                    get_dims(stack[sp]),
-                                                    .!ismissing.(stack[sp]))
-Sₚ_to_Sₐ(series::RasterSeries, sp::Symbol) = Sₚ_to_Sₐ.(series, sp)
+                                                    get_dims(stack[Sₚ]),
+                                                    .!ismissing.(stack[Sₚ]))
+Sₚ_to_Sₐ(series::RasterSeries, Sₚ::Symbol) = Sₚ_to_Sₐ.(series, Sₚ)
 
 """
     function θ_to_Θ(raster::Raster, Sₐ::raster, rs_dims::Tuple, find_nm::Raster)
@@ -123,7 +134,7 @@ Convert a `Raster` of potential temperature (`θ`) to conservative temperature (
 function θ_to_Θ(θ::Raster, Sₐ::Raster, rs_dims::Tuple, find_nm::Raster)
 
     lons, lats, z, time = rs_dims
-    Θ = similar(Array(θ), Union{Float64, Missing})
+    Θ = similar(Array(θ))
 
     if isnothing(time)
 
@@ -139,16 +150,26 @@ function θ_to_Θ(θ::Raster, Sₐ::Raster, rs_dims::Tuple, find_nm::Raster)
     return Raster(Θ, rs_dims)
 
 end
-θ_to_Θ(stack::RasterStack, pt::Symbol, sp::Symbol) = θ_to_Θ(stack[pt],
-                                                            Sₚ_to_Sₐ(stack, sp),
-                                                            get_dims(stack[pt]),
-                                                            .!ismissing.(stack[pt]) .&&
-                                                            .!ismissing.(stack[sp]))
-θ_to_Θ(series::RasterSeries, pt::Symbol, sp::Symbol) = θ_to_Θ.(series, pt, sp)
+"""
+    function θ_to_Θ(stack::RasterStack, var_names::NamedTuple)
+    function θ_to_Θ(series::RasterSeries, var_names::NamedTuple)
+Convert the potential temperature, returning only the conservative temperature, from a
+`RasterStack` or a `RasterSeries`. The `var_names` must be passed in as for
+`convert_ocean_vars` ---  that is, as a named tuple in the form
+`(Sₚ = :salt_name, θ = :potential_temp_name)` where `:potential_temp_name` and
+`:salt_name` are the name of the potential temperature and salinity in the `Raster`.
+"""
+θ_to_Θ(stack::RasterStack, var_names::NamedTuple) = θ_to_Θ(stack[var_names.θ],
+                                                           Sₚ_to_Sₐ(stack, var_names.Sₚ),
+                                                           get_dims(stack[var_names.θ]),
+                                                           .!ismissing.(stack[var_names.θ]) .&&
+                                                           .!ismissing.(stack[var_names.Sₚ]))
+θ_to_Θ(series::RasterSeries, var_names::NamedTuple) = θ_to_Θ.(series, Ref(var_names))
 
 """
     function in_situ_density(Sₐ::Raster, Θ::Raster, p::Raster, rs_dims::Tuple, find_nm::Raster)
-Compute in-situ density using `gsw_rho` from GibbsSeaWater.jl.
+Compute in-situ density using `gsw_rho` from GibbsSeaWater.jl. This conversion depends on
+absolute salinity (`Sₐ`), conservative temperature (`Θ`) and pressure (`p`).
 """
 function in_situ_density(Sₐ::Raster, Θ::Raster, p::Raster, rs_dims::Tuple, find_nm::Raster)
 
@@ -169,12 +190,27 @@ function in_situ_density(Sₐ::Raster, Θ::Raster, p::Raster, rs_dims::Tuple, fi
     return Raster(ρ, rs_dims)
 
 end
+"""
+    function in_situ_density(stack::RasterStack, var_names::NamedTuple)
+    function in_situ_density(series::RasterStack, var_names::NamedTuple)
+Compute and return the in-situ density `ρ` from a `RasterStack` or `RasterSeries`.
+This computation depends on absolute salinity `Sₐ`, conservative temperature `Θ`
+and pressure `p`. The variable names must be passed into the function as a `NamedTuple` in
+the form `(Sₐ = :salt_var, Θ = :temp_var, p = :pressure_var)`.
+"""
+in_situ_density(stack::RasterStack, var_names::NamedTuple) =
+    in_situ_density(stack[var_names.Sₐ], stack[var_names.Θ], stack[var_names.p],
+                    get_dims(stack[var_names.Sₐ]),
+                   .!ismissing.(stack[var_names.Sₐ]) .&& .!ismissing.(stack[var_names.Θ]))
+in_situ_density(series::RasterSeries, var_names::NamedTuple) = in_situ_density.(series, Ref(var_names))
 
 """
     function potential_density(Sₐ::Raster, Θ::Raster, p::Float64, rs_dims::Tuple, find_nm::Raster)
 Compute potential density at reference pressure `p`, `σₚ`, using `gsw_rho`  from GibbsSeaWater.jl.
+This conversion depends on absolute salinity (`Sₐ`), conservative temperature (`Θ`) and a
+user entered reference pressure (`p`).
 """
-function potential_density(Sₐ::Raster, Θ::Raster, p::Float64, rs_dims::Tuple, find_nm::Raster)
+function potential_density(Sₐ::Raster, Θ::Raster, p::Number, rs_dims::Tuple, find_nm::Raster)
 
     lons, lats, z, time = rs_dims
     σₚ = similar(Array(Θ))
@@ -193,7 +229,20 @@ function potential_density(Sₐ::Raster, Θ::Raster, p::Float64, rs_dims::Tuple,
     return Raster(σₚ, rs_dims)
 
 end
-
+"""
+    function potential_density(stack::RasterStack, var_names::NamedTuple)
+    function potential_density(series::RasterStack, var_names::NamedTuple)
+Compute and return the potential density `σₚ` at reference pressure `p` from a
+`RasterStack` or `RasterSeries`. This computation depends on absolute salinity `Sₐ`,
+conservative temperature `Θ` and a reference pressure `p`. The variable names must be
+passed into the function as a `NamedTuple` in the form
+`(Sₐ = :salt_var, Θ = :temp_var, p = ref_pressure)`. Note `p` in this case is a number.
+"""
+potential_density(stack::RasterStack, var_names::NamedTuple) =
+    potential_density(stack[var_names.Sₐ], stack[var_names.Θ], var_names.p,
+                      get_dims(stack[var_names.Sₐ]),
+                      .!ismissing.(stack[var_names.Sₐ]) .&& .!ismissing.(stack[var_names.Θ]))
+potential_density(series::RasterSeries, var_names::NamedTuple) = potential_density.(series, Ref(var_names))
 """
     function get_dims(raster::Raster)
 Get the dimensions of a `Raster`.

test/runtests.jl

@@ -10,6 +10,10 @@ include("test_oceanrasterconversions.jl")
     @test isequal(converted_Sₚ, Sₐ_)
     ## `θ_to_Θ`
     @test isequal(converted_θ, Θ)
+    ## `in_situ_density`
+    @test isequal(converted_ρ_stack, ρ)
+    ## `potential_density`
+    @test isequal(converted_σₚ_stack, σₚ)
     ## `convert_ocean_vars`
     # In situ density
     for (i, var) ∈ enumerate(test_vars_in_situ)
@@ -29,6 +33,10 @@ end
     @test isequal(converted_Sₚ_series, Sₐ_)
     ## `θ_to_Θ`
     @test isequal(converted_θ_series, Θ)
+     ## `in_situ_density`
+     @test isequal(converted_ρ_series, ρ)
+     ## `potential_density`
+     @test isequal(converted_σₚ_series, σₚ)
 
     ## `convert_ocean_vars`
     # In situ density
@@ -48,8 +56,8 @@ end
 
 @testset "Argument errors" begin
 
-    @test_throws ArgumentError convert_ocean_vars(rs_stack_NoX, (sp = :Sₚ, pt = :θ))
-    @test_throws ArgumentError convert_ocean_vars(rs_stack_NoY, (sp = :Sₚ, pt = :θ))
-    @test_throws ArgumentError convert_ocean_vars(rs_stack_NoZ, (sp = :Sₚ, pt = :θ))
+    @test_throws ArgumentError convert_ocean_vars(rs_stack_NoX, (Sₚ = :Sₚ, θ = :θ))
+    @test_throws ArgumentError convert_ocean_vars(rs_stack_NoY, (Sₚ = :Sₚ, θ = :θ))
+    @test_throws ArgumentError convert_ocean_vars(rs_stack_NoZ, (Sₚ = :Sₚ, θ = :θ))
 
 end

test/test_oceanrasterconversions.jl

@@ -30,17 +30,17 @@ Sₚ_noZ = rand(Sₚ_vals, X(lons), Y(lats), Ti(time))
 θ_noZ = rand(θ_vals, X(lons), Y(lats), Ti(time))
 test_vars_noZ = (Sₚ = Sₚ_noZ, θ = θ_noZ)
 rs_stack_NoZ = RasterStack(test_vars_noZ, (X(lons), Y(lats), Ti(time)))
-
+Sₚ
 ## Output to test
 converted_p = depth_to_pressure(rs_stack)
 converted_Sₚ = Sₚ_to_Sₐ(rs_stack, :Sₚ)
 converted_Sₚ_series = Rasters.combine(Sₚ_to_Sₐ(rs_series, :Sₚ), Ti)
-converted_θ = θ_to_Θ(rs_stack, :θ, :Sₚ)
-converted_θ_series = Rasters.combine(θ_to_Θ(rs_series, :θ, :Sₚ), Ti)
-rs_stack_res_in_situ = convert_ocean_vars(rs_stack, (sp = :Sₚ, pt = :θ))
-rs_stack_res_pd = convert_ocean_vars(rs_stack, (sp = :Sₚ, pt = :θ); ref_pressure)
-rs_series_res_in_situ = convert_ocean_vars(rs_series, (sp = :Sₚ, pt = :θ))
-rs_series_res_pd = convert_ocean_vars(rs_series, (sp = :Sₚ, pt = :θ); ref_pressure)
+converted_θ = θ_to_Θ(rs_stack, (Sₚ = :Sₚ, θ = :θ))
+converted_θ_series = Rasters.combine(θ_to_Θ(rs_series, (Sₚ = :Sₚ, θ = :θ)), Ti)
+rs_stack_res_in_situ = convert_ocean_vars(rs_stack, (Sₚ = :Sₚ, θ = :θ))
+rs_stack_res_pd = convert_ocean_vars(rs_stack, (Sₚ = :Sₚ, θ = :θ); ref_pressure)
+rs_series_res_in_situ = convert_ocean_vars(rs_series, (Sₚ = :Sₚ, θ = :θ))
+rs_series_res_pd = convert_ocean_vars(rs_series, (Sₚ = :Sₚ, θ = :θ); ref_pressure)
 
 test_vars_in_situ = keys(rs_stack_res_in_situ)
 test_vars_pd = keys(rs_stack_res_pd)
@@ -78,5 +78,12 @@ for t ∈ time
     end
 end
 
+test_stack = RasterStack((Sₐ = Sₐ, Θ = Θ, p = p), (X(lons), Y(lats), Z(z), Ti(time)))
+test_series = RasterSeries([test_stack[Ti(t)] for t ∈ time], Ti)
+converted_ρ_stack = in_situ_density(test_stack, (Sₐ = :Sₐ, Θ = :Θ, p = :p))
+converted_ρ_series = Rasters.combine(in_situ_density(test_series, (Sₐ = :Sₐ, Θ = :Θ, p = :p)), Ti)
+converted_σₚ_stack = potential_density(test_stack, (Sₐ = :Sₐ, Θ = :Θ, p = ref_pressure))
+converted_σₚ_series = Rasters.combine(potential_density(test_series, (Sₐ = :Sₐ, Θ = :Θ, p = ref_pressure)), Ti)
+
 vars_in_situ = (p, Sₐ, Θ, ρ)
 vars_pd = (p, Sₐ, Θ, σₚ)