Testbroadcast

Project.toml

@@ -19,6 +19,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
+Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 TimeZones = "f269a46b-ccf7-5d73-abea-4c690281aa53"
 
 [compat]
@@ -35,6 +36,7 @@ RecursiveArrayTools = "2.20.0"
 StaticArrays = "1.2.13"
 StatsBase = "0.33.12"
 Suppressor = "0.2"
+Tables = "1.6.0"
 TimeZones = "1.6.1"
 julia = "1.6"
 

inst/explore_broadcast_missing.jl

@@ -0,0 +1,14 @@
+x = [1.0, 2.0, missing]
+y = [3.0, 3.0, 3.0]
+f(x,y) = x+y
+function g(x,y)
+    z = f.(x,y)
+    if z isa Union
+        return typeof(z)
+        convert(typeof(z).b.a, z)
+    else
+        z
+    end
+end
+z = g(x,y)
+@code_warntype g(x,y)

src/Bigleaf.jl

@@ -4,6 +4,7 @@ module Bigleaf
 using Optim
 using FillArrays
 using DataFrames
+using Tables: columntable
 using Dates, TimeZones
 using Pipe
 using AstroLib
@@ -37,6 +38,9 @@ export Monin_Obukhov_length, Monin_Obukhov_length!, stability_parameter,
     roughness_parameters, Reynolds_Number
 export aerodynamic_conductance!, add_Ga, add_Ga!, compute_Ram, compute_Ram!, roughness_z0h
 
+# shorthand type: either DataFrame or column table (Tables.columntable)
+DFTable = Union{NamedTuple,AbstractDataFrame}
+
 include("util.jl")    
 include("bigleaf_constants.jl")
 include("unit_conversions.jl")

src/bigleaf_constants.jl

@@ -39,6 +39,7 @@ Default values can be overridden by the named arguments:
 - frac2percent Conversion between fraction and percent
 """
 function bigleaf_constants(;
+  # make sure to provide Float64 values (rather than Int) so that return is type-stable
   ## Physical constants
   cp         = 1004.834,        # specific heat of air for constant pressure (J K-1 kg-1)
   Rgas       = 8.31451,         # universal gas constant (J mol-1 K-1)
@@ -49,7 +50,7 @@ function bigleaf_constants(;
   eps        = 0.622,           # ratio of the molecular weight of water vapor to dry air (=Mw/Md)
   g          = 9.81,            # gravitational acceleration (m s-2)
   solar_constant = 1366.1,      # solar constant, i.e. solar radation at earth distance from the sun (W m-2)
-  pressure0  = 101325,          # reference atmospheric pressure at sea level (Pa)
+  pressure0  = 101325.0,        # reference atmospheric pressure at sea level (Pa)
   Tair0      = 273.15,          # reference air temperature (K)
   k          = 0.41,            # von Karman constant
   Cmol       = 0.012011,        # molar mass of carbon (kg mol-1)
@@ -61,17 +62,17 @@ function bigleaf_constants(;
   ## Conversion constants
   Kelvin       = 273.15,         # conversion degree Celsius to Kelvin
   DwDc         = 1.6,            # Ratio of the molecular diffusivities for water vapor and CO2
-  days2seconds = 86400,          # seconds per day
-  kPa2Pa       = 1000,           # conversion kilopascal (kPa) to pascal (Pa)
+  days2seconds = 86400.0,        # seconds per day
+  kPa2Pa       = 1000.0,         # conversion kilopascal (kPa) to pascal (Pa)
   Pa2kPa       = 0.001,          # conversion pascal (Pa) to kilopascal (kPa)
-  umol2mol     = 1e-06,          # conversion micromole (umol) to mole (mol)
-  mol2umol     = 1e06,           # conversion mole (mol) to micromole (umol)
-  kg2g         = 1000,           # conversion kilogram (kg) to gram (g)
+  umol2mol     = 1.0e-6,         # conversion micromole (umol) to mole (mol)
+  mol2umol     = 1.0e6,          # conversion mole (mol) to micromole (umol)
+  kg2g         = 1000.0,         # conversion kilogram (kg) to gram (g)
   g2kg         = 0.001,          # conversion gram (g) to kilogram (kg)
-  kJ2J         = 1000,           # conversion kilojoule (kJ) to joule (J)
+  kJ2J         = 1000.0,         # conversion kilojoule (kJ) to joule (J)
   J2kJ         = 0.001,          # conversion joule (J) to kilojoule (kJ)
   se_median    = 1.253,          # conversion standard error (SE) of the mean to SE of the median (http://influentialpoints.com/Training/standard_error_of_median.htm)
-  frac2percent = 100             # conversion between fraction and percent
+  frac2percent = 100.0           # conversion between fraction and percent
 )
   Dict(
 :cp => cp, :Rgas => Rgas, :Rv => Rv, :Rd => Rd, :Md => Md, :Mw => Mw, :eps => eps, :g => g,

src/boundary_layer_conductance.jl

@@ -41,10 +41,13 @@ true
 ``` 
 """
 function compute_Gb!(df::AbstractDataFrame, approach::Val{:Thom_1972}; kwargs...)
-  compute_Gb_!(df, approach, :ustar; kwargs...) # inputcols
+  ft(ustar) = Gb_Thom(ustar; kwargs...)
+  transform!(df, :ustar => ByRow(ft) => AsTable)
 end
-function compute_Gb!(df::AbstractDataFrame, approach::Val{:constant_kB1}; kwargs...)
-  compute_Gb_!(df, approach, :ustar; kwargs...) # inputcols
+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...)
+  transform!(df, :ustar => ft => AsTable)
 end
 function compute_Gb!(df::AbstractDataFrame, approach::Val{:Choudhury_1988}; kwargs...)
   Gb_Choudhury!(df; kwargs...)
@@ -53,14 +56,6 @@ function compute_Gb!(df::AbstractDataFrame, approach::Val{:Su_2001}; kwargs...)
   Gb_Su!(df; kwargs...)
 end
 
-function compute_Gb_!(df::AbstractDataFrame, approach, inputcols; 
-  Sc::AbstractVector=SA[], kwargs...
-  )
-  fGb(args...) = compute_Gb(approach, args...; kwargs...)
-  transform!(df, :ustar => ByRow(fGb) => SA[:Rb_h, :Gb_h, :kB_h, :Gb_CO2])
-end
-compute_Gb(::Val{:Thom_1972}, args...; kwargs...) = Gb_Thom(args...; kwargs...)
-compute_Gb(::Val{:constant_kB1}, args...; kB_h, kwargs...) = Gb_constant_kB1(args..., kB_h; kwargs...)
 
 """
     add_Gb(Gb_h::Union{Missing,Number}, Sc::Vararg{Pair,N}; constants)
@@ -270,8 +265,8 @@ function Gb_Choudhury!(df::AbstractDataFrame; leafwidth, LAI, wind_zh=nothing,
     wind_zh = wind_profile(df, zh, d, z0m; zh, zr, stab_formulation, constants)
   end
   # Broadcasting does not work over keyword arguments, need to pass as positional
-  fwind(ustar, wind_zh; kwargs...) = Gb_Choudhury(ustar;wind_zh, kwargs...)
-  ft(ustar) = fwind.(ustar, wind_zh; leafwidth, LAI, constants)
+  fwind(ustar, wind_zh, leafwidth, LAI; kwargs...) = Gb_Choudhury(ustar;wind_zh, leafwidth, LAI, kwargs...)
+  ft(ustar) = fwind.(ustar, wind_zh, leafwidth, LAI; constants)
   transform!(df, :ustar => ft => AsTable)
 end
 
@@ -365,7 +360,7 @@ compute_Gb!(df,Val(:Su_2001), zh=25,zr=40,Dl=0.1,LAI=1.5)
 true
 ``` 
 """
-function Gb_Su(Tair,pressure,ustar; zh, wind_zh, Dl, fc, N=2, Cd=0.2, hs=0.01,
+function Gb_Su(Tair,pressure,ustar; wind_zh, Dl, fc, N=2, Cd=0.2, hs=0.01,
   constants=bigleaf_constants()
   )
   wind_zh = max(0.01,wind_zh) ## avoid zero windspeed
@@ -390,13 +385,13 @@ function Gb_Su!(df::AbstractDataFrame; wind_zh=nothing, Dl, fc=nothing,
   end
   if isnothing(fc)
     isnothing(LAI) && error("one of 'fc' or 'LAI' must be provided")
-    fc = (1-exp(-LAI/2)) 
+    fc = length(LAI) == 1 ? (1-exp(-LAI/2)) : @. (1-exp(-LAI/2))
   end
   isnothing(Dl) && error(
     "need to provide keyword argument Dl with :Su_2001 method")
   inputcols = SA[:Tair,:pressure,:ustar]
   # Broadcasting does not work over keyword arguments, need to pass as positional
-  fwind(wind_zh, args...; kwargs...) = Gb_Su(args...; wind_zh, kwargs...)
-  ft(args...) = fwind.(wind_zh, args...; zh, Dl, fc, N=2, Cd=0.2, hs=0.01, constants)
+  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)
   transform!(df, inputcols => ft => AsTable)
 end

src/stability_correction.jl

@@ -56,9 +56,10 @@ 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)
 end
-function Monin_Obukhov_length(df;constants=bigleaf_constants())
-  Monin_Obukhov_length!(copy(df, copycols = false); 
-    constants).MOL
+function Monin_Obukhov_length(df::DFTable; kwargs...)
+    tmp = Monin_Obukhov_length.(df.Tair, df.pressure, df.ustar, df.H; kwargs...)
+    # eltype_agg = Union{eltype(df.Tair), eltype(df.pressure), eltype(df.ustar), eltype(df.H)}
+    # convert(Vector{eltype_agg}, tmp)
 end
 
 """
@@ -110,7 +111,7 @@ function stability_parameter!(df::AbstractDataFrame; zr,d, MOL=nothing,
   df[!,:zeta] .= stability_parameter(df; zr,d,MOL,constants)
   df
 end
-function stability_parameter(df::AbstractDataFrame; zr,d, MOL=nothing,
+function stability_parameter(df::DFTable; zr,d, MOL=nothing,
   constants=bigleaf_constants())
   if isnothing(MOL) MOL = Monin_Obukhov_length(df; constants); end
   stability_parameter.(zr,d,MOL)
@@ -219,14 +220,22 @@ function get_stability_coefs_unstable(::Val{:Dyer_1970}, zeta)
 end
 
 #TODO z or zr here?
-function stability_correction(Tair,pressure,ustar,H, z,d; 
+function stability_correction(Tair::Union{Missing,Number},pressure,ustar,H, z,d; 
   stab_formulation=Val(:Dyer_1970), constants = bigleaf_constants())
   stab_formulation isa Val{:no_stability_correction} && return(
     (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)
 end
+function stability_correction(Tair,pressure,ustar,H, z,d; kwargs...) 
+  Tables.columns(
+    stability_correction.(Tair,pressure,ustar,H, z,d; kwargs...))
+end
+function stability_correction(df::DFTable, z,d; kwargs...)   
+  tmp1 = stability_correction.(df.Tair, df.pressure, df.ustar, df.H, z, d; kwargs...)
+  tmp2 = Tables.columns(tmp1)
+end
 
 function stability_correction!(df; zeta=df.zeta, 
   stab_formulation=Val(:Dyer_1970), constants = bigleaf_constants())
@@ -242,6 +251,6 @@ end
 
 function stability_correction!(df, z, d; 
   stab_formulation=Val(:Dyer_1970), constants = bigleaf_constants())
-  ft(args...) = stability_correction(args..., z, d; stab_formulation, constants)
-  transform!(df, SA[:Tair,:pressure,:ustar,:H] => ByRow(ft) => AsTable)
+  ft(args...) = stability_correction.(args..., z, d; stab_formulation, constants)
+  transform!(df, SA[:Tair,:pressure,:ustar,:H] => ft => AsTable)
 end

src/surface_roughness.jl

@@ -142,7 +142,8 @@ function roughness_parameters(::Val{:canopy_height}, zh; frac_d=0.7, frac_z0m=0.
   (;d, z0m, z0m_se)
 end
 
-function roughness_parameters(::Val{:canopy_height_LAI}, zh, LAI; cd=0.2, hs=0.01)
+function roughness_parameters(::Val{:canopy_height_LAI}, zh, LAI; 
+  cd=0.2, hs=0.01)
   X = cd * LAI
   d = 1.1 * zh * log(1 + X^(1/4))
   z0m = ifelse(0 <= X <= 0.2, hs + 0.3 * X^(1/2), 0.3 * zh * (1 - d/zh))
@@ -152,13 +153,12 @@ end
 
 # docu: supply psi_m = 0 for no stability correction, default method
 # if psi_m is given df only needs wind and ustar
-function roughness_parameters(::Val{:wind_profile}, df, zh, zr;
+function roughness_parameters(::Val{:wind_profile}, df::DFTable, zh, zr;
   d = 0.7*zh, psi_m = nothing, stab_formulation=Val(:Dyer_1970), 
   constants=bigleaf_constants()
   )
   if isnothing(psi_m)
-    psi_m = stability_correction!(copy(df, copycols=false), zr, d; 
-      stab_formulation, constants).psi_m
+    psi_m = stability_correction(df, zr, d; stab_formulation, constants).psi_m
   end
   z0m_all = allowmissing(@. (zr - d) * exp(-constants[:k]*df.wind / df.ustar - psi_m))
   #z0m_all[(z0m_all .> zh)] .= missing # problems with missings
@@ -263,7 +263,7 @@ function wind_profile(stab_formulation, z, ustar, Tair,pressure,H, d, z0m;
   wind_profile(z, ustar, d, z0m, psi_m)
 end
 
-function wind_profile(df::AbstractDataFrame, z, d, z0m = nothing; 
+function wind_profile(df::DFTable, z, d, z0m = nothing; 
   zh = nothing, zr = nothing, 
   stab_formulation = Val(:Dyer_1970), constants = bigleaf_constants()
   )
@@ -275,13 +275,12 @@ function wind_profile(df::AbstractDataFrame, z, d, z0m = nothing;
     # uses psi_m at zr
     z0m = roughness_parameters(Val(:wind_profile), df, zh, zr).z0m
   end
-  psi_m = stability_correction!(
-    copy(df, copycols=false), z, d; stab_formulation, constants).psi_m
+  psi_m = stability_correction(df, z, d; stab_formulation, constants).psi_m
   #wind_profile(df, z, d, z0m, psi_m; constants)
   wind_profile(df, z, d, z0m, psi_m; constants)
 end
 
-function wind_profile(df::AbstractDataFrame, z, d, z0m, psi_m::AbstractVector;
+function wind_profile(df::DFTable, z, d, z0m, psi_m::AbstractVector;
   constants = bigleaf_constants())
   # when psi_m is given, df is not used any more, but keep for consistency
   wind_profile.(z, df.ustar, d, z0m, psi_m; constants)

test/Project.toml

@@ -7,5 +7,6 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
+Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TimeZones = "f269a46b-ccf7-5d73-abea-4c690281aa53"