use anonymous functions with DataFrames to avoid world-age problems

docs/src/walkthrough.md

@@ -405,7 +405,7 @@ wind_heights = 22:2:60.0
 d = 0.7 * thal.zh
 z0m = roughness_parameters(Val(:wind_profile), thas; zh=thal.zh, zr=thal.zr).z0m
 wp = map(wind_heights) do z
-  wind_profile(z, thas,d, z0m; zh=thal.zh, zr=thal.zr)
+  wind_profile(z, thas,d, z0m)
 end;
 nothing # hide
 ```

inst/explore_broadcast_missing.jl

@@ -1,5 +1,5 @@
-x = [1.0, 2.0, missing]
-y = [3.0, 3.0, 3.0]
+x = [1.0, 2.0, missing, missing]
+y = [3.0, 3.0, 3.0, 3.0]
 f(x,y) = x+y
 function g(x,y)
     z = f.(x,y)
@@ -11,4 +11,27 @@ function g(x,y)
     end
 end
 z = g(x,y)
-@code_warntype g(x,y)
+@code_warntype g(x,y)
+
+
+typeof(x), typeof(x[1:2]), typeof(x[3:4]) # preserves eltype
+typeof(x.*1.0), typeof(x[1:2].*1.0), typeof(x[3:4].*1.0) # does not
+typeof(x*1.0), typeof(x[1:2]*1.0), typeof(x[3:4]*1.0) # does not
+
+
+x = [1.0, 2.0, missing, missing]
+y = 1.0
+typeof(x.*y), typeof(x[1:2].*y), typeof(x[3:4].*y) 
+# (Vector{Union{Missing, Float64}}, Vector{Float64}, Vector{Missing})
+res = x[1:2].*y
+res_eltype = Union{eltype(x),typeof(y)}
+convert(Vector{res_eltype}, res)::Vector{res_eltype}
+
+f1(x,y) = x[1:2].*y
+@code_warntype f1(x,y)
+
+function f2(x,y)
+    res_eltype = Union{eltype(x),typeof(y)}
+    convert(Vector{res_eltype}, x[1:2].*y)::Vector{res_eltype}
+end
+@inferred f2(x,y)

src/aerodynamic_conductance.jl

@@ -102,7 +102,7 @@ function aerodynamic_conductance!(df; Gb_model = Val(:Thom_1972), Ram_model = Va
     if isnothing(z0m) 
       z0m = roughness_parameters(Val(:wind_profile), df; zh, zr, psi_m = df.psi_m).z0m
     end
-    wind_zh = wind_profile(zh, df, d, z0m; zh, zr, stab_formulation, constants)
+    wind_zh = wind_profile(zh, df, d, z0m; stab_formulation, constants)
   end
   #
   # calculate canopy boundary layer conductance (Gb)
@@ -117,14 +117,14 @@ function aerodynamic_conductance!(df; Gb_model = Val(:Thom_1972), Ram_model = Va
   # calculate aerodynamic risistance for momentum (Ra_m)
   Ram_model isa Val{:wind_profile} && compute_Ram!(df, Ram_model; zr, d, z0m, constants) 
   Ram_model isa Val{:wind_zr} && compute_Ram!(df, Ram_model) 
-  function ft(Ra_m, Gb_h, Gb_CO2) 
+  fr = (Ra_m, Gb_h, Gb_CO2) -> begin
     Ga_m = 1/Ra_m
     Ra_h = Ra_m + 1/Gb_h
     Ga_h = 1/Ra_h
     Ga_CO2 =  1/(Ra_m + 1/Gb_CO2)
     (;Ga_m, Ga_h, Ra_h, Ga_CO2)
   end
-  transform!(df, [:Ra_m, :Gb_h, :Gb_CO2] => ByRow(ft) => AsTable)
+  transform!(df, [:Ra_m, :Gb_h, :Gb_CO2] => ByRow(fr) => AsTable)
 end
 
 
@@ -230,16 +230,16 @@ function compute_Ram!(df, method::Val{:wind_profile};
   # put keyword arguments to positional arguments for proper broadcast
   ftpos(ustar, zr, d, z0m, psi_h) = 
     compute_Ram(method, ustar; zr, d, z0m, psi_h, kwargs...)
-  ft(ustar) = ftpos.(ustar, zr, d, z0m, psi_h)
+  ft = (ustar) -> ftpos.(ustar, zr, d, z0m, psi_h)
   transform!(df, :ustar => ft => :Ra_m)
 end
 
 function compute_Ram(::Val{:wind_zr}, ustar::Union{Missing,Number}, wind)
   Ra_m = wind / ustar^2
 end
 function compute_Ram!(df, method::Val{:wind_zr}; kwargs...)
-  ft(ustar, wind) = compute_Ram(method, ustar, wind; kwargs...)
-  transform!(df, SA[:ustar, :wind] => ByRow(ft) => :Ra_m)
+  fr = (ustar, wind) -> compute_Ram(method, ustar, wind; kwargs...)
+  transform!(df, SA[:ustar, :wind] => ByRow(fr) => :Ra_m)
 end
 
 
@@ -286,7 +286,7 @@ function add_Ga!(df::AbstractDataFrame, Sc::Vararg{Pair,N};
   Gb_h = df.Gb_h, Ga_m = df.Ga_m, kwargs...) where N
   N == 0 && return(df)
   Scn, Scv = get_names_and_values("Ga_", Sc...)
-  ft() = add_Ga_.(Gb_h, Ga_m, Ref(Scn), Ref(Scv); kwargs...)
+  ft = () -> add_Ga_.(Gb_h, Ga_m, Ref(Scn), Ref(Scv); kwargs...)
   transform!(df, [] => ft => AsTable)
 end
 function add_Ga(Gb_h::Union{Missing,Number}, Ga_m::Union{Missing,Number}, 

src/boundary_layer_conductance.jl

@@ -36,12 +36,12 @@ true
 ``` 
 """
 function compute_Gb!(df::AbstractDataFrame, approach::Val{:Thom_1972}; kwargs...)
-  ft(ustar) = Gb_Thom(ustar; kwargs...)
-  transform!(df, :ustar => ByRow(ft) => :Gb_h)
+  fr = (ustar) -> Gb_Thom(ustar; kwargs...)
+  transform!(df, :ustar => ByRow(fr) => :Gb_h)
 end
 function compute_Gb!(df::AbstractDataFrame, approach::Val{:constant_kB1}; kB_h, kwargs...)
   # do not use ByRow because kb_H can be a vector
-  ft(ustar) = Gb_constant_kB1.(ustar, kB_h; kwargs...)
+  ft = (ustar) -> Gb_constant_kB1.(ustar, kB_h; kwargs...)
   transform!(df, :ustar => ft => :Gb_h)
 end
 function compute_Gb!(df::AbstractDataFrame, approach::Val{:Choudhury_1988}; kwargs...)
@@ -76,7 +76,7 @@ function compute_Gb_quantities(Gb_h, ustar; constants=bigleaf_constants())
   (;Rb_h, Gb_h, kB_h, Gb_CO2)
 end
 function compute_Gb_quantities!(df::AbstractDataFrame; constants=bigleaf_constants())
-  ft(args...) = compute_Gb_quantities.(args...; constants)
+  ft = (args...) -> compute_Gb_quantities.(args...; constants)
   transform!(df, SA[:Gb_h, :ustar] => ft => AsTable)
 end
 
@@ -121,7 +121,7 @@ true
 function add_Gb!(df::AbstractDataFrame, Sc::Vararg{Pair,N}; Gb_h = df.Gb_h, kwargs...) where N
   N == 0 && return(df)
   Scn, Scv = get_names_and_values("Gb_", Sc...)
-  ft() = add_Gb_.(Gb_h, Ref(Scn), Ref(Scv); kwargs...)
+  ft = () -> add_Gb_.(Gb_h, Ref(Scn), Ref(Scv); kwargs...)
   transform!(df, [] => ft => AsTable)
 end
 function add_Gb(Gb_h::Union{Missing,Number}, Sc::Vararg{Pair,N}; kwargs...) where N
@@ -266,7 +266,7 @@ function Gb_Choudhury!(
   # end
   # Broadcasting does not work over keyword arguments, need to pass as positional
   fwind(wind_zh, leafwidth, LAI; kwargs...) = Gb_Choudhury(;wind_zh, leafwidth, LAI, kwargs...)
-  ft() = fwind.(wind_zh, leafwidth, LAI; constants)
+  ft = () -> fwind.(wind_zh, leafwidth, LAI; constants)
   transform!(df, [] => ft => :Gb_h)
 end
 
@@ -333,7 +333,9 @@ to Su et al. 2001:
 using DataFrames
 df = DataFrame(Tair=25,pressure=100,wind=[3,4,5],ustar=[0.5,0.6,0.65],H=[200,230,250]) 
 zh = 25; zr = 40
-wind_zh = wind_profile(zh, df, 0.7*zh; zh, zr)
+z0m = roughness_parameters(
+  Val(:wind_profile), df.ustar, df.wind, df.Tair, df.pressure, df.H; zh, zr).z0m 
+wind_zh = wind_profile(zh, df, 0.7*zh, z0m)
 compute_Gb!(df,Val(:Su_2001); wind_zh, Dl=0.01, LAI=5)
 # the same meteorological conditions, but larger leaves
 compute_Gb!(df,Val(:Su_2001); wind_zh, Dl=0.1,LAI=5)
@@ -370,6 +372,6 @@ function Gb_Su!(df::AbstractDataFrame; wind_zh, Dl, fc=nothing,
   inputcols = SA[:Tair,:pressure,:ustar]
   # Broadcasting does not work over keyword arguments, need to pass as positional
   fwind(wind_zh, Dl, fc, N, Cd, hs, args...; kwargs...) = Gb_Su(args...; wind_zh, Dl, fc, N, Cd, hs, kwargs...)
-  ft(args...) = fwind.(wind_zh, Dl, fc, N, Cd, hs, args...; constants)
+  ft = (args...) -> fwind.(wind_zh, Dl, fc, N, Cd, hs, args...; constants)
   transform!(df, inputcols => ft => :Gb_h)
 end

src/stability_correction.jl

@@ -53,8 +53,8 @@ function Monin_Obukhov_length(Tair, pressure, ustar, H; constants=bigleaf_consta
   MOL  = (-rho*constants[:cp]*ustar^3*TairK) / (constants[:k]*constants[:g]*H)
 end
 function Monin_Obukhov_length!(df;constants=bigleaf_constants())
-  ft(args...) = Monin_Obukhov_length(args...; constants)
-  transform!(df, SA[:Tair, :pressure, :ustar, :H] => ByRow(ft) => :MOL)
+  fr = (args...) -> Monin_Obukhov_length(args...; constants)
+  transform!(df, SA[:Tair, :pressure, :ustar, :H] => ByRow(fr) => :MOL)
 end
 # function Monin_Obukhov_length(df::DFTable; kwargs...)
 #     tmp = Monin_Obukhov_length.(df.Tair, df.pressure, df.ustar, df.H; kwargs...)
@@ -113,11 +113,11 @@ end
 function stability_parameter!(df::AbstractDataFrame; z,d, MOL=nothing,
   constants=bigleaf_constants())
   if isnothing(MOL)
-    ft(args...) = stability_parameter.(z,d,args...; constants)
+    ft = (args...) -> stability_parameter.(z,d,args...; constants)
     transform!(df, SA[:Tair, :pressure, :ustar, :H] => ft => :zeta)
   else
-    ft2() = stability_parameter.(z,d,MOL)
-    transform!(df, [] => ft2 => :zeta)
+    ft = () -> stability_parameter.(z,d,MOL)
+    transform!(df, [] => ft => :zeta)
   end
 end
 # function stability_parameter(df::DFTable; z,d, MOL=nothing,
@@ -234,7 +234,8 @@ function stability_correction(z,d, Tair::Union{Missing,Number},pressure,ustar,H;
     (psi_h = zero(z), psi_m = zero(z)))
   MOL = Monin_Obukhov_length(Tair,pressure,ustar,H; constants)
   zeta  = stability_parameter(z,d,MOL)
-  stability_correction(zeta; stab_formulation)
+  psis = stability_correction(zeta; stab_formulation)
+  psis
 end
 # function stability_correction(Tair,pressure,ustar,H, z,d; kwargs...) 
 #   Tables.columns(
@@ -253,20 +254,18 @@ function stability_correction!(df; zeta=nothing, z=nothing, d=nothing,
     return(df)
   end
   if !isnothing(zeta)
-    ft() = stability_correction.(zeta; stab_formulation)
+    ft = () -> stability_correction.(zeta; stab_formulation)
     transform!(df, [] => ft => AsTable)
   else
-    function ft_met(args...)
-      stability_correction.(z, d, args...; stab_formulation, constants)
-    end
-    transform!(df, SA[:Tair,:pressure,:ustar,:H] => ft_met => AsTable)
+    ft = (args...) -> stability_correction.(z, d, args...; stab_formulation, constants)
+    transform!(df, SA[:Tair,:pressure,:ustar,:H] => ft => AsTable)
   end
 end
 
 function stability_correction(df::DFTable; z, d, 
   stab_formulation=Val(:Dyer_1970), constants = bigleaf_constants()
   )
-  # cannot dispatch on keyword argument, hence need if-clause
+  # do not provide zeta, because can simply invoke stability_correction.(zeta)
   if stab_formulation isa Val{:no_stability_correction}
     z = zero(first(skipmissing(df.ustar)))
     rows = map(Tables.rows(df)) do row

src/surface_roughness.jl

@@ -177,7 +177,10 @@ function roughness_parameters(method::Val{:wind_profile},
   #   zr, d, Tair, pressure, ustar, H; stab_formulation, constants)).psi_m
   psi_m = getindex.(stability_correction.(
     zr, d, Tair, pressure, ustar, H; stab_formulation, constants), :psi_m)
-  roughness_parameters(method, ustar, wind, psi_m; zh, zr, d, constants, kwargs...)
+  # https://discourse.julialang.org/t/eltype-changed-during-broadcast-with-missing/71312?u=bgctw
+  res_eltype = Union{eltype(ustar),eltype(wind),eltype(Tair),eltype(pressure),eltype(H)}
+  psi_m = convert(Vector{res_eltype}, psi_m)::Vector{res_eltype}    
+  rp = roughness_parameters(method, ustar, wind, psi_m; zh, zr, d, constants, kwargs...)
 end
 
 function roughness_parameters(method::Val{:wind_profile}, df::DFTable; 
@@ -194,15 +197,11 @@ function roughness_parameters(method::Val{:wind_profile}, df::DFTable;
   end
 end
 
-
-
-
-
+
 """                                                                                                                       
     wind_profile(z::Number, ustar, d, z0m, psi_m = zero(z); constants)
-    wind_profile(df::AbstractDataFrame, z, d, z0m, psi_m::AbstractVector; constants)
-    wind_profile(df::AbstractDataFrame, z, d, z0m = nothing; 
-      zh = nothing, zr = nothing, 
+    wind_profile(z, df::AbstractDataFrame, d, z0m, psi_m::AbstractVector; constants)
+    wind_profile(z, df::DFTable, d, z0m; psi_m = nothing,
       stab_formulation = Val(:Dyer_1970), constants)
 
 Wind speed at a given height above the canopy estimated from single-level
@@ -219,13 +218,9 @@ For DataFrame variant with supplying stability_parameter
   - `ustar`     : Friction velocity (m s-1)
 - `psi_m`     : value of the stability function for heat, see [`stability_correction`](@ref)
   Pass `psi_m = 0.0` to neglect stability correction.
-For DataFrame varinat where psi_m and z0m are to be estimated
-- `zh`        : canopy height (m)
-- `zr`        : Instrument (reference) height (m)
-- `df`:       : DataFrame with columns 
-  - `ustar`     : Friction velocity (m s-1)
+For DataFrame varinat where psi_m is to be estimated
+- additional columns in `df`:
   - `Tair`, `pressure`, `H` : see [`stability_correction`](@ref)
-  - `wind` : see [`roughness_parameters(Val(:wind_profile), ...)`]
 
 # Details
 The underlying assumption is the existence of a logarithmic wind profile
@@ -238,10 +233,8 @@ In this case, the wind speed at a given height z is given by:
 The roughness parameters zero-plane displacement height (d) and roughness length (z0m)
 can be approximated from [`roughness_parameters`](@ref). 
 ``\\psi_m`` is the stability correction. Set it to zero (not recommended) to neglect
-statbility correction.
-If `z0m` is not provided, i.e. defaults to nothing, then 
-`z0m` is first estimated from the wind profile equation and then used in the equation
-above for the calculation of `u(z)` (see e.g. Newman & Klein 2014).        
+statbility correction. By default it is estimated from wind profile using
+[`stability_correction`](@ref)
 
 # Note
 Note that this equation is only valid for z >= d + z0m, and it is not 
@@ -266,9 +259,12 @@ using DataFrames
 heights = 18:2:40  # heights above ground for which to calculate wind speed
 df = DataFrame(Tair=25,pressure=100,wind=[3,4,5],ustar=[0.5,0.6,0.65],H=[200,230,250]) 
 zr=40;zh=25;d=16
-z0m = roughness_parameters(Val(:wind_profile), df; zh, zr).z0m
+# z0m and MOL are independent of height, compute before
+MOL = Monin_Obukhov_length.(df.Tair, df.pressure, df.ustar, df.H)
+z0m = roughness_parameters(
+  Val(:wind_profile), df.ustar, df.wind, df.Tair, df.pressure, df.H; zh, zr).z0m 
 ws = map(heights) do z
-  wind_profile(z,df,d,z0m; zr, zh)
+  wind_profile(z,df,d,z0m; MOL)
 end
 using Plots # plot wind profiles for the three rows in df
 plot(first.(ws), heights, ylab = "height (m)", xlab = "wind speed (m/s)", legend=:topleft)
@@ -278,31 +274,31 @@ nothing
 # output
 ``` 
 """  
-function wind_profile(z::Number, ustar, d, z0m, psi_m; constants=bigleaf_constants())
+function wind_profile(
+  z::Number, ustar::Union{Missing,Number}, d, z0m, psi_m; constants=bigleaf_constants()
+  )
   wind_heights = max(0,(ustar / constants[:k]) * (log(max(0,(z - d)) / z0m) - psi_m))
 end
 
-function wind_profile(z, ustar::Union{Missing,Number}, d, z0m, Tair,pressure,H,
+function wind_profile(z::Number, ustar::Union{Missing,Number}, d, z0m, Tair,pressure,H,
   stab_formulation=Val(:Dyer_1970), constants=bigleaf_constants())
-  psi_m = stability_correction(z,d, Tair,pressure,ustar,H; stab_formulation, constants).psi_m
+  psi_m = stability_correction(
+    z,d, Tair,pressure,ustar,H; stab_formulation, constants).psi_m
   wind_profile(z, ustar, d, z0m, psi_m)
 end
 
-function wind_profile(z, df::DFTable, d, z0m = nothing; 
-  zh = nothing, zr = nothing, psi_m = nothing,
+function wind_profile(z::Number, df::DFTable, d, z0m; psi_m = nothing, MOL = nothing,
   stab_formulation = Val(:Dyer_1970), constants = bigleaf_constants()
   )
-  if isnothing(z0m)
-    (isnothing(zh) || isnothing(zr)) && error(
-      "wind_profile: If 'z0m' is not given, must specify both 'zh' and 'zr' optional " *
-      "arguments to estimate 'z0m' from 'df.wind' and 'df.ustar' measured at height zr.")
-    # uses psi_m at zr
-    z0m = roughness_parameters(Val(:wind_profile), df; zh, zr).z0m
-  end
   if isnothing(psi_m)
-    psi_m = stability_correction(df; z, d, stab_formulation, constants).psi_m
+    if isnothing(MOL) 
+      psi_m = stability_correction(df; z, d, stab_formulation, constants).psi_m
+    else
+      zeta = stability_parameter.(z,d,MOL)
+      psi_m = getindex.(stability_correction.(zeta), :psi_m)
+    end
   end
   #wind_profile(df, z, d, z0m, psi_m; constants)
-  wind_profile.(z, df.ustar, d, z0m, psi_m; constants)
+  windz = wind_profile.(z, df.ustar, d, z0m, psi_m; constants)
 end
 

test/boundary_layer_conductance.jl

@@ -81,7 +81,25 @@ end
     #
     # DataFrame variant
     df = copy(tha48)
-    wind_zh = wind_profile(zh, df, 0.7*zh; zh, zr)
+    # sum(skipmissing(...)) not inferrable in Julia 1.6
+    # @descend_code_warntype roughness_parameters(
+    #     Val(:wind_profile), df.ustar, df.wind, df.Tair, df.pressure, df.H; 
+    #     zh, zr)
+    # @code_warntype roughness_parameters(
+    #     Val(:wind_profile), df.ustar, df.wind, df.Tair, df.pressure, df.H; 
+    #     zh, zr)
+    # z0m = (@inferred roughness_parameters(
+    #     Val(:wind_profile), df.ustar, df.wind, df.Tair, df.pressure, df.H; 
+    #     zh, zr)).z0m
+    z0m = roughness_parameters(
+        Val(:wind_profile), df.ustar, df.wind, df.Tair, df.pressure, df.H; 
+        zh, zr).z0m
+    # function f1(zh, ustar, z0m, Tair, pressure, H) 
+    #     wind_zh = wind_profile.(zh, ustar, 0.7*zh, z0m, Tair, pressure, H)
+    #     wind_zh * 2
+    # end
+    # @code_warntype f(zh, df.ustar, z0m,df.Tair, df.pressure, df.H)    
+    wind_zh = wind_profile.(zh, df.ustar, 0.7*zh, z0m, df.Tair, df.pressure, df.H)
     @inferred compute_Gb!(df, Val(:Choudhury_1988); leafwidth, LAI, wind_zh)
     @test last(propertynames(df)) == :Gb_h
     @test df.Gb_h[1] ≈ Gb_Choud rtol=1e-6

test/surface_roughness.jl

@@ -85,18 +85,15 @@ end
     @test windzc[1] == u30c
     #plot(windz)
     #plot!(windz2)
+    #
+    # with providing psi_m
     psi_m = stability_correction(df; z, d).psi_m
     windzc2 = @inferred wind_profile(z, columntable(dfd), d, z0m; psi_m)    
     @test windzc2 == windzc
     #
-    # estimate z0m
-    # need to give zh and zr in addition to many variables in df
-    @test_throws Exception wind_profile(z, df, d)    
-    windzc3 = @inferred wind_profile(z, columntable(dfd), d; zh=thal.zh, zr=thal.zr, 
-        stab_formulation = Val(:Dyer_1970))    
-    # may have used slightly different estimated z0m
-    #windzc3 - windzc
-    @test all(isapprox.(windzc3, windzc, atol=0.1))
-    @test windzc3[1] ≈ 2.764203 rtol=1e-3 # from R
+    # with providing L
+    MOL = Monin_Obukhov_length.(df.Tair, df.pressure, df.ustar, df.H)
+    windzc3 = @inferred wind_profile(z, columntable(dfd), d, z0m; MOL)    
+    @test windzc3 == windzc
 end