From 6376e6154f31473879b3f63497051936a33bd21c Mon Sep 17 00:00:00 2001 From: Alexander Barth Date: Fri, 24 Mar 2023 11:39:57 +0100 Subject: [PATCH 1/3] Update README.md --- README.md | 1 - 1 file changed, 1 deletion(-) diff --git a/README.md b/README.md index adfdb09..edc7913 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,6 @@ [![Build Status](https://github.com/gher-uliege/PhysOcean.jl/workflows/CI/badge.svg)](https://github.com/gher-uliege/PhysOcean.jl/actions) [![Build Status Windows](https://ci.appveyor.com/api/projects/status/github/gher-uliege/PhysOcean.jl?branch=master&svg=true)](https://ci.appveyor.com/project/Alexander-Barth/physocean-jl) [![codecov.io](http://codecov.io/github/gher-uliege/PhysOcean.jl/coverage.svg?branch=master)](http://codecov.io/github/gher-uliege/PhysOcean.jl?branch=master) - [![documentation stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://gher-uliege.github.io/PhysOcean.jl/stable/) [![documentation latest](https://img.shields.io/badge/docs-latest-blue.svg)](https://gher-uliege.github.io/PhysOcean.jl/latest/) From 08305013b814bd1e738a8d0b998af914a74dec21 Mon Sep 17 00:00:00 2001 From: Zhiwei Date: Wed, 19 Apr 2023 15:37:38 +0800 Subject: [PATCH 2/3] modify the nansum and nanmean function --- src/PhysOcean.jl | 116 ++++++++++++++++++++++++----------------------- 1 file changed, 59 insertions(+), 57 deletions(-) diff --git a/src/PhysOcean.jl b/src/PhysOcean.jl index 6aef217..f1965eb 100644 --- a/src/PhysOcean.jl +++ b/src/PhysOcean.jl @@ -20,27 +20,29 @@ const OMEGA = 7.2921150E-5 # mean radius of earth (A.E. Gill, 1982. Atmosphere-Ocean Dynamics) const MEAN_RADIUS_EARTH = 6371000 # m -function nansum(x) - return sum(x[.!isnan.(x)]) +function nanmean(x::AbstractArray{T}) where {T<:Number} + sum = zero(T) + count = 0 + for i in eachindex(x) + if !isnan(x[i]) + sum += x[i] + count += 1 + end + end + return sum / count end -function nansum(x,dim) - m = isnan.(x) - x2 = copy(x) - x2[m] .= 0 - return sum(x2,dims = dim) +function nanmean(x::AbstractArray{T}, dims) where {T<:Number} + return mapslices(nanmean, x; dims=dims) end -function nanmean(x) - return mean(x[.!isnan.(x)]) -end -function nanmean(x,dim) - m = isnan.(x) - x2 = copy(x) - x2[m] .= 0 +function nansum(arr::AbstractArray{T}) where {T<:Number} + return (sum(x -> !isnan(x) * x, arr)) +end - return sum(x2,dims = dim) ./ sum(.!m,dims = dim) +function nansum(arr::AbstractArray{T}, dims) where {T<:Number} + return (sum(x -> !isnan(x) * x, arr, dims=dims)) end @@ -52,7 +54,7 @@ Return DateTime from matlab's and octave's datenum function datetime_matlab(datenum) # even if datenum is Int32, the computations must be done with Int64 to # prevent overflow - return DateTime(1970,1,1) + Dates.Millisecond(round(Int64,(datenum-Int64(719529)) * 24*60*60*1000)) + return DateTime(1970, 1, 1) + Dates.Millisecond(round(Int64, (datenum - Int64(719529)) * 24 * 60 * 60 * 1000)) end @@ -68,22 +70,22 @@ the relative humidity `r` (0 ≤ r ≤ 1, pressure ratio, not percentage), the wind speed `w` (m/s) and the air pressure (hPa). """ -function latentflux(Ts,Ta,r,w,pa) +function latentflux(Ts, Ta, r, w, pa) - Da = 1.5e-3; - rhoa = 1.3; # kg m-3 - Lh = 2.5e6; # J - epsilon = 0.622; + Da = 1.5e-3 + rhoa = 1.3 # kg m-3 + Lh = 2.5e6 # J + epsilon = 0.622 - esta = vaporpressure(Ta); - ests = vaporpressure(Ts); + esta = vaporpressure(Ta) + ests = vaporpressure(Ts) - qqa = r * esta * epsilon / pa; - qqs = ests * epsilon / pa; + qqa = r * esta * epsilon / pa + qqs = ests * epsilon / pa - Qe = Da * rhoa * abs(w) * (qqs-qqa) * Lh; + Qe = Da * rhoa * abs(w) * (qqs - qqa) * Lh return Qe end @@ -96,16 +98,16 @@ the sea surface temperature `Ts` (degree Celsius), the air temperature `Ta` (degree Celsius), the wate vapour pressure `e` (hPa) and the total cloud coverage `ttc` (0 ≤ tcc ≤ 1). """ -function longwaveflux(Ts,Ta,e,tcc) - epsilon = 0.98; - sigma = 5.67e-8; - lambda = 0.69; +function longwaveflux(Ts, Ta, e, tcc) + epsilon = 0.98 + sigma = 5.67e-8 + lambda = 0.69 # degree C to degree K Ts = Ts + TK Ta = Ta + TK - Qb = epsilon * sigma * Ts^4 * (0.39-0.05*sqrt(e))*(1-lambda*tcc^2)+4 * epsilon * sigma * Ts^3 *(Ts-Ta) + Qb = epsilon * sigma * Ts^4 * (0.39 - 0.05 * sqrt(e)) * (1 - lambda * tcc^2) + 4 * epsilon * sigma * Ts^3 * (Ts - Ta) return Qb end @@ -118,12 +120,12 @@ the wind speed `w` (m/s), the sea surface temperature `Ts` (degree Celsius), the air temperature `Ta` (degree Celsius). """ -function sensibleflux(Ts,Ta,w) - Sta = 1.45e-3; - ca = 1000; - rhoa = 1.3; +function sensibleflux(Ts, Ta, w) + Sta = 1.45e-3 + ca = 1000 + rhoa = 1.3 - Qc = Sta * ca * rhoa * w * (Ts-Ta); + Qc = Sta * ca * rhoa * w * (Ts - Ta) return Qc end @@ -133,8 +135,8 @@ end Compute the solar heat flux (W/m²) """ -function solarflux(Q,al) - Qs = Q*(1-al) +function solarflux(Q, al) + Qs = Q * (1 - al) return Qs end @@ -148,7 +150,7 @@ Monteith, J.L., and Unsworth, M.H. 2008. Principles of Environmental Physics. Th """ function vaporpressure(T) # Monteith and Unsworth (2008), https://en.wikipedia.org/wiki/Tetens_equation - e = 6.1078 * exp((17.27 * T)./(T + 237.3)); + e = 6.1078 * exp((17.27 * T) ./ (T + 237.3)) return e end @@ -159,9 +161,9 @@ end Return a Gaussian window with `N` points with a standard deviation of (N-1)/(2 α). """ -function gausswin(N, α = 2.5) - sigma = (N-1)/(2 * α) - return [exp(- n^2 / (2*sigma^2)) for n = -(N-1)/2:(N-1)/2] +function gausswin(N, α=2.5) + sigma = (N - 1) / (2 * α) + return [exp(-n^2 / (2 * sigma^2)) for n = -(N - 1)/2:(N-1)/2] end """ @@ -169,19 +171,19 @@ end Filter the vector `x` with a `N`-point Gaussian window. """ -function gaussfilter(x,N) - b = gausswin(N); - c = b/sum(b); - imax = size(x,1); +function gaussfilter(x, N) + b = gausswin(N) + c = b / sum(b) + imax = size(x, 1) xf = zeros(imax) - s=1; + s = 1 - for i=1:imax - xf[i] = sum(x[s:s+N-1] .* c); + for i = 1:imax + xf[i] = sum(x[s:s+N-1] .* c) - s = s+1; - if s>=size(x,1)-N - s=size(x,1)-N; + s = s + 1 + if s >= size(x, 1) - N + s = size(x, 1) - N end end @@ -194,7 +196,7 @@ end Provides coriolisfrequency et given latidudes in DEGREES from -90 Southpole to +90 Northpole """ function coriolisfrequency(latitude) - return 2*OMEGA*sind(latitude) + return 2 * OMEGA * sind(latitude) end """ @@ -203,9 +205,9 @@ end Provides gravity in m/s2 at ocean surface at given latidudes in DEGREES from -90 Southpole to +90 Northpole """ function earthgravity(latitude) - return 9.780327*(1.0026454-0.0026512* - cosd(2*latitude)+0.0000058*(cosd(2*latitude))^2 - ) + return 9.780327 * (1.0026454 - 0.0026512 * + cosd(2 * latitude) + 0.0000058 * (cosd(2 * latitude))^2 + ) end @@ -237,7 +239,7 @@ end # return (m,eofs,relvar,totvar) # end -function addprefix!(prefix,obsids) +function addprefix!(prefix, obsids) if prefix == "" return end From ee5b793db8d1096acf3e67e9e102bbe75dc0a256 Mon Sep 17 00:00:00 2001 From: Alexander Barth Date: Wed, 19 Apr 2023 10:45:00 +0200 Subject: [PATCH 3/3] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 0f3f994..392dfcd 100644 --- a/Project.toml +++ b/Project.toml @@ -4,7 +4,7 @@ keywords = ["physical-oceanography", "sea-water", "julia", "density", "fluxes", license = "LGPL" desc = "Utility functions for physical oceanography" author = ["Alexander Barth", "Jean-Marie Beckers", "Charles Troupin", "Aida Alvera Azcarate"] -version = "0.6.7" +version = "0.6.8" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"