diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml new file mode 100644 index 0000000..cd83f87 --- /dev/null +++ b/.github/workflows/CI.yml @@ -0,0 +1,41 @@ +# https://github.com/Alexander-Barth/NCDatasets.jl/blob/master/.github/workflows/ci.yml +name: CI +on: + - push + - pull_request +jobs: + test: + name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + version: + # - '1.3' + - '1' + # - 'nightly' + os: + - ubuntu-latest + - macOS-latest + # - windows-latest + arch: + - x64 + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: julia-actions/cache@v1 + - uses: julia-actions/julia-buildpkg@latest + - run: | + git config --global user.name Tester + git config --global user.email te@st.er + - uses: julia-actions/julia-runtest@latest + continue-on-error: ${{ matrix.version == 'nightly' }} + + - uses: julia-actions/julia-processcoverage@v1 + - uses: codecov/codecov-action@v4 + with: + token: ${{ secrets.CODECOV_TOKEN }} + file: lcov.info diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml new file mode 100644 index 0000000..eda5134 --- /dev/null +++ b/.github/workflows/CompatHelper.yml @@ -0,0 +1,24 @@ +name: CompatHelper + +on: + schedule: + - cron: '15 23 */4 * *' + +jobs: + CompatHelper: + runs-on: ${{ matrix.os }} + strategy: + matrix: + julia-version: [1.2.0] + julia-arch: [x86] + os: [ubuntu-latest] + steps: + - uses: julia-actions/setup-julia@latest + with: + version: ${{ matrix.julia-version }} + - name: Pkg.add("CompatHelper") + run: julia -e 'using Pkg; Pkg.add("CompatHelper")' + - name: CompatHelper.main() + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + run: julia -e 'using CompatHelper; CompatHelper.main()' diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml new file mode 100644 index 0000000..28beba9 --- /dev/null +++ b/.github/workflows/Documenter.yml @@ -0,0 +1,20 @@ +name: Documenter +on: + push: + branches: [master, main] + tags: [v*] + pull_request: + branches: [master, main] +jobs: + Documenter: + name: Documentation + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v3 + - uses: julia-actions/julia-buildpkg@latest + - uses: julia-actions/julia-docdeploy@latest + env: + RASTERDATASOURCES_PATH: ".." + GKSwstype: "100" + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml new file mode 100644 index 0000000..edccfdc --- /dev/null +++ b/.github/workflows/TagBot.yml @@ -0,0 +1,18 @@ +# https://github.com/JuliaRegistries/TagBot +name: TagBot +on: + issue_comment: + types: + - created + workflow_dispatch: +jobs: + TagBot: + if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot' + runs-on: ubuntu-latest + steps: + - uses: JuliaRegistries/TagBot@v1 + with: + token: ${{ secrets.GITHUB_TOKEN }} + # Edit the following line to reflect the actual name of the GitHub Secret containing your private key + ssh: ${{ secrets.DOCUMENTER_KEY }} + # ssh: ${{ secrets.NAME_OF_MY_SSH_PRIVATE_KEY_SECRET }} diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..5a730b8 --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +Manifest.toml +*debug* +.vscode +docs/build diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..0887a9b --- /dev/null +++ b/Project.toml @@ -0,0 +1,13 @@ +name = "ClimateStats" +uuid = "f2027f19-6765-4257-993d-235c5fe36d77" +authors = ["Dongdong Kong "] +version = "0.1.0" + +[deps] +Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" +DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +Ipaper = "e58298cb-69f7-4186-aecd-5834b6793426" +NaNStatistics = "b946abbf-3ea7-4610-9019-9858bfdeaf2d" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 0000000..d2b75cf --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,6 @@ +[deps] +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 0000000..7393101 --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,34 @@ +using Documenter, ClimateStats +CI = get(ENV, "CI", nothing) == "true" + +# Logging.disable_logging(Logging.Warn) + +# Make the docs, without running the tests again +# We need to explicitly add all the extensions here +makedocs( + modules=[ + ClimateStats + # Base.get_extension(Ipaper, :IpaperSlopeExt) + ], + format=Documenter.HTML( + prettyurls=CI, + ), + pages=[ + "Introduction" => "index.md" + "R Base" => "RBase.md" + "Statistics" => "Statistics.md" + "Slope" => "Slope.md" + "Climate" => "Climate.md" + "Parallel" => "Parallel.md" + ], + sitename="ClimateStats.jl", + warnonly=true, + clean=false, +) + +# Enable logging to console again +# Logging.disable_logging(Logging.BelowMinLevel) + +deploydocs( + repo="github.com/jl-pkgs/ClimateStats.jl.git", +) diff --git a/docs/src/Climate.md b/docs/src/Climate.md new file mode 100644 index 0000000..fa67492 --- /dev/null +++ b/docs/src/Climate.md @@ -0,0 +1,36 @@ +```@contents +Pages = ["index.md"] +Depth = 3 +``` + +## Thresholds + +```@docs +cal_anomaly_quantile +cal_anomaly_clim +cal_threshold +``` + +```@docs +_cal_anomaly_3d +``` + +```@docs +cal_mTRS_base +cal_mTRS_full + +cal_climatology_base +cal_climatology_full +``` + + +```@docs +Ipaper.cal_mTRS_base! +Ipaper.cal_mTRS_base3! +``` + + +```@docs +cal_yearly_Tair +cal_warming_level +``` diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 0000000..521ae4a --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,15 @@ +# ClimateStats in Julia (R base for Julia) + +[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://jl-pkgs.github.io/ClimateStats.jl/stable) +[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://jl-pkgs.github.io/ClimateStats.jl/dev) +[![CI](https://github.com/jl-pkgs/ClimateStats.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/jl-pkgs/ClimateStats.jl/actions/workflows/CI.yml) +[![Codecov](https://codecov.io/gh/jl-pkgs/ClimateStats.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/jl-pkgs/ClimateStats.jl) + +> Dongdong Kong + +## Installation + +``` +using Pkg +Pkg.add(url="https://github.com/jl-pkgs/ClimateStats.jl") +``` diff --git a/src/ClimateStats.jl b/src/ClimateStats.jl new file mode 100644 index 0000000..389f3ee --- /dev/null +++ b/src/ClimateStats.jl @@ -0,0 +1,27 @@ +module ClimateStats + +using Ipaper +using Printf +using Dates +using NaNStatistics + +using Reexport +@reexport using DocStringExtensions: TYPEDSIGNATURES, METHODLIST + +include("base.jl") +include("threshold.jl") +include("climatology.jl") +include("anomaly.jl") +include("warming_level.jl") + +# export cal_mTRS_base, cal_mTRS_season, cal_mTRS_full +export format_md +export cal_climatology_base, cal_climatology_full, cal_climatology_full, + _cal_anomaly, + cal_anomaly_quantile, cal_anomaly_clim, cal_threshold, + cal_warming_level, cal_yearly_Tair + + +export filter_mds + +end # module ClimateStats diff --git a/src/anomaly.jl b/src/anomaly.jl new file mode 100644 index 0000000..9d29dd0 --- /dev/null +++ b/src/anomaly.jl @@ -0,0 +1,240 @@ +## 计算anomaly的3种方法，考虑每年的升温幅度 +export _cal_anomaly_3d, _cal_anomaly + + +""" +$(TYPEDSIGNATURES) +""" +function _cal_anomaly_3d( + A::AbstractArray{T,3}, + TRS::AbstractArray{T,3}, + dates; + T_wl::Union{AbstractArray{T,3},Nothing}=nothing, + fun_anom=_exceed, ignored... +) where {T<:Real} + + mmdd = format_md.(dates) + mds = mmdd |> unique_sort + + years = year.(dates) + year_grps = years |> unique_sort + + ind_d = indexin(mmdd, mds) + ind_y = indexin(years, year_grps) + + _wl = T_wl === nothing ? T(0) : T_wl[:, :, ind_y] + fun_anom.(A, TRS[:, :, ind_d], _wl) +end + + +function _cal_anomaly( + A::AbstractArray{T,Na}, + TRS::AbstractArray{T,Nt}, + dates; + T_wl::Union{AbstractArray{T,Na},Nothing}=nothing, + fun_anom=_exceed, + deep=true, + ignored... +) where {T<:Real,Na,Nt} + + mmdd = format_md.(dates) + mds = mmdd |> unique_sort + + years = year.(dates) + year_grps = years |> unique_sort + + ind_d = indexin(mmdd, mds) + ind_y = indexin(years, year_grps) + + _wl = T_wl === nothing ? T(0) : selectdim_deep(T_wl, Na, ind_y; deep) + fun_anom.(A, selectdim_deep(TRS, Na, ind_d; deep), _wl) + # broadcast(fun_anom, A, deepcopy(selectdim(TRS, Na, ind_d)), _wl) +end + + +""" +$(TYPEDSIGNATURES) + +Calculate the anomaly of a 3D array of temperature data. + +# Arguments + +- `A` : the 3D array of temperature data +- `dates` : an array of dates corresponding to the temperature data +- `parallel` : whether to use parallel processing (default `true`) +- `use_mov` : whether to use a moving window to calculate the threshold + (default `true`) +- `method` : the method to use for calculating the threshold, one of `["full", + "season", "base", "pTRS"]` (default `"full"`) +- `probs` : default `[0.5]` +- `p1` : the start year for the reference period (default `1981`) +- `p2` : the end year for the reference period (default `2010`) +- `fun` : the function used to calculate the anomaly (default `_exceed`) + +# Returns + +An array of the same shape as `A` containing the temperature anomaly. + +# References + +1. Vogel, M. M., Zscheischler, J., Fischer, E. M., & Seneviratne, S. I. (2020). + Development of Future Heatwaves for Different Hazard Thresholds. Journal of + Geophysical Research: Atmospheres, 125(9). + + +""" +function cal_anomaly_quantile( + A::AbstractArray{T}, dates; + parallel::Bool=true, + use_mov=true, na_rm=true, + method="full", + p1=1981, p2=2010, + fun=_exceed, + probs=[0.5], + options... +) where {T<:Real} + + kw = (; probs, use_mov, na_rm, parallel, options...) + if method == "base" + mTRS = cal_mTRS_base(A, dates; p1, p2, kw...) |> squeeze_tail + anom = _cal_anomaly(A, mTRS, dates) + elseif method == "season" + mTRS = cal_mTRS_base(A, dates; p1, p2, kw...) |> squeeze_tail + T_wl = cal_warming_level(A, dates; p1, p2) + anom = _cal_anomaly(A, mTRS, dates; T_wl) + elseif method == "full" + TRS_full = cal_mTRS_full(A, dates; kw...) |> squeeze_tail + anom = fun.(A, TRS_full) + + elseif method == "pTRS" + # 最基础的方法 + years = year.(dates) + inds = p1 .<= years .<= p2 + data = selectdim(A, ndims(A), inds) + anom = map(prob -> begin + TRS = NanQuantile(data, [prob]) + fun.(A, TRS) + end, probs) + anom = cat(anom..., dims=ndims(A) + 1) |> squeeze_tail + end + anom +end + + +""" +$(TYPEDSIGNATURES) + +Calculate the anomaly of an array relative to its climatology. + +# Arguments + +- `A::AbstractArray{T}`: The input array to calculate the anomaly of. +- `dates`: The dates corresponding to the input array. + +- `parallel::Bool=true`: Whether to use parallel processing. +- `use_mov=true`: Whether to use a moving window to calculate the climatology. +- `method="full"`: The method to use for calculating the climatology. Can be "base", "season", or "full". +- `p1=1981`: The start year for the period to use for calculating the climatology. +- `p2=2010`: The end year for the period to use for calculating the climatology. +- `fun_clim=nanmean`: The function to use for calculating the climatology. +- `fun_anom=_exceed`: The function to use for calculating the anomaly. + +# Returns +- `anom`: The anomaly of the input array relative to its climatology. + +# Example +```julia +using Ipaper + +# Generate some sample data +A = rand(365, 10) +dates = Date(2000, 1, 1):Day(1):Date(2000, 12, 31) + +# Calculate the anomaly relative to the climatology +anom = cal_anomaly_clim(A, dates; method="base") +``` +""" +function cal_anomaly_clim( + A::AbstractArray{T}, + dates; + parallel::Bool=true, + use_mov=true, + method="full", + p1=1981, p2=2010, + fun_clim=nanmean, + fun_anom=_exceed +) where {T<:Real} + + kw = (; use_mov, parallel, fun=fun_clim) + + if method == "base" + mTRS = cal_climatology_base(A, dates; p1, p2, kw...) |> squeeze_tail + anom = _cal_anomaly(A, mTRS, dates; fun_anom) + elseif method == "season" + mTRS = cal_climatology_base(A, dates; p1, p2, kw...) |> squeeze_tail + T_wl = cal_warming_level(A, dates; p1, p2) + anom = _cal_anomaly(A, mTRS, dates; T_wl, fun_anom) + elseif method == "full" + TRS_full = cal_climatology_full(A, dates; kw...) |> squeeze_tail + anom = fun_anom.(A, TRS_full) + end + anom +end + + +""" +$(TYPEDSIGNATURES) + +Calculate the threshold value for a given dataset `A` and dates. The threshold value is calculated based on the specified method. + +# Arguments +- `A::AbstractArray{T}`: The input data array. +- `dates`: The dates corresponding to the input data array. + +- `parallel::Bool=true`: Whether to use parallel computation. +- `use_mov::Bool=true`: Whether to use moving window. +- `na_rm::Bool=true`: Whether to remove missing values. +- `method::String="full"`: Possible values are "base", "season", and "full". +- `p1::Int=1981`: The start year for the reference period. +- `p2::Int=2010`: The end year for the reference period. +- `probs::Vector{Float64}=[0.5]`: The probability levels to use for calculating the threshold value. +- `options...`: Additional options to pass to the underlying functions. + +# Returns + +For different methods: + +- `full`: Array with the dimension of `(dims..., ntime, nprob)` +- `base`: Array with the dimension of `(dims..., 366, nprob)` +- `season`: Array with the dimension of `(dims..., nyear)` + +# Examples +```julia +dates = Date(2010, 1):Day(1):Date(2020, 12, 31); +ntime = length(dates) +data = rand(10, ntime); +cal_threshold(data, dates; p1=2010, p2=2015, method="full") +``` +""" +function cal_threshold( + A::AbstractArray{T}, dates; + parallel::Bool=true, + use_mov=true, na_rm=true, + method="full", + p1=1981, p2=2010, + probs::Vector=[0.90, 0.95, 0.99, 0.999, 0.9999], + options... +) where {T<:Real} + + kw = (; probs, use_mov, na_rm, parallel, options...) + if method == "base" + cal_mTRS_base(A, dates; p1, p2, kw...) |> squeeze_tail + elseif method == "season" + # mTRS = cal_mTRS_base(A, dates; p1, p2, kw...) |> squeeze_tail + cal_warming_level(A, dates; p1, p2) + elseif method == "full" + cal_mTRS_full(A, dates; kw...) |> squeeze_tail # full + elseif method == "pTRS" + cal_pTRS(A, dates; p1, p2, probs, options...) + end +end diff --git a/src/base.jl b/src/base.jl new file mode 100644 index 0000000..169a4dc --- /dev/null +++ b/src/base.jl @@ -0,0 +1,56 @@ +using Dates +using Printf: @sprintf + +# DateType = Union{Date,DateTime,AbstractCFDateTime,Nothing} +_gte(x::T, trs::T, wl::T=T(0)) where {T<:Real} = x >= (trs + wl) +_gt(x::T, trs::T, wl::T=T(0)) where {T<:Real} = x > (trs + wl) +_exceed(x::T, trs::T, wl::T=T(0)) where {T<:Real} = x - (trs + wl) + + +# format_md(date) = @sprintf("%02d-%02d", month(date), day(date)) +format_md(date) = month(date)*100 + day(date) + + +# function find_adjacent_doy2(doy::Int; doy_max::Int=366, halfwin::Int=7) +# ind = (-halfwin:halfwin) .+ doy + +# if ind[end] > doy_max +# [doy-halfwin:doy_max, 1:ind[end]-doy_max] +# elseif ind[1] < 1 +# [1:ind[end], ind[1]+doy_max:doy_max] +# else +# [ind] +# end +# end + +function find_adjacent_doy(doy::Int; doy_max::Int=366, halfwin::Int=7) + ind = collect(-halfwin:halfwin) .+ doy + for i = eachindex(ind) + if ind[i] > doy_max + ind[i] = ind[i] - doy_max + end + if ind[i] <= 0 + ind[i] = ind[i] + doy_max + end + end + ind +end + + +function filter_mds(mmdd::AbstractVector, doy::Int; doy_max::Integer=366, halfwin::Int=7, use_mov=true) + !use_mov && (return mmdd .== doy) + + ind = (-halfwin:halfwin) .+ doy + if ind[end] > doy_max + # ind1, ind2 = doy-halfwin:doy_max, 1:ind[end]-doy_max + @.(doy - halfwin <= mmdd <= doy_max || 1 <= mmdd <= ind[end] - doy_max) + elseif ind[1] < 1 + # ind1, ind2 = 1:ind[end], ind[1]+doy_max:doy_max + @.(1 <= mmdd <= ind[end] || ind[1] + doy_max <= mmdd <= doy_max) + else + @.(ind[1] <= mmdd <= ind[end]) + end +end + + +export format_md, find_adjacent_doy, filter_mds diff --git a/src/climatology.jl b/src/climatology.jl new file mode 100644 index 0000000..0c168fd --- /dev/null +++ b/src/climatology.jl @@ -0,0 +1,86 @@ +# """ +# $(TYPEDSIGNATURES) +# """ +# function cal_climatology_base3!(Q::AbstractArray{T,3}, A::AbstractArray{T,3}, mmdd; +# use_mov=true, halfwin::Int=7, +# parallel::Bool=true, fun=nanmean, +# ignored...) where {T<:Real} + +# doy_max = maximum(mmdd) + +# nlon, nlat, _ = size(A) +# @inbounds @par parallel for doy = 1:doy_max +# ind = filter_mds(mmdd, doy; doy_max, halfwin, use_mov) +# for i = 1:nlon, j = 1:nlat +# x = @view A[i, j, ind] +# Q[i, j, doy] = fun(x) +# end +# end +# Q +# end + + +function cal_climatology_base!(Q::AbstractArray{T,N}, A::AbstractArray{T,N}, mmdd; + use_mov=true, halfwin::Int=7, + dims=N, + parallel::Bool=true, fun=nanmean, + ignored...) where {T<:Real,N} + + doy_max = maximum(mmdd) + + @par parallel for doy = 1:doy_max + ind = filter_mds(mmdd, doy; doy_max, halfwin, use_mov) + + _data = selectdim(A, dims, ind) + _Q = selectdim(Q, dims, doy) + _Q .= mapslices(fun, _data; dims) + end + Q +end + + +""" +$(TYPEDSIGNATURES) + +Calculate the climatology of a dataset `A` based on the `dates`. + +The climatology is the long-term average of a variable over a specific period of time. +This function calculates the climatology of the input dataset `A` based on the dates `dates`. +The calculation is performed by applying a function `fun!` to a moving window of the data. + +Arguments: +- `A : :AbstractArray{T}`: the input dataset, where `T` is a subtype of `Real`. +- `dates` : the dates associated with the input dataset, as a vector of `Date` objects. +- `fun!` : the function to apply to the moving window of the data. It should take an input array and return a scalar. +- `use_quantile`: default false, a boolean indicating whether to use a quantile-based filter to remove outliers. +- `p1`, `p2` : the references period + +Returns: +- a matrix of the same size as `A`, containing the climatology values. + +Example: +```julia +using Dates +using Statistics +# using NaNStatistics + +A = rand(365, 10) |> transpose # simulate a year of daily data for 10 variables +dates = Date(2022, 1, 1):Day(1):Date(2022, 12, 31) +clim = cal_climatology_base(A, dates; p1=2022, p2=2022) +``` +""" +function cal_climatology_base(A::AbstractArray{T}, dates; + (fun!)=cal_climatology_base!, kw...) where {T<:Real} + + cal_mTRS_base(A, dates; use_quantile=false, (fun!), kw...) +end + + +""" +$(TYPEDSIGNATURES) +""" +function cal_climatology_full(A::AbstractArray{T}, dates; + (fun!)=cal_climatology_base!, kw...) where {T<:Real} + + cal_mTRS_full(A, dates; use_quantile=false, (fun!), kw...) +end diff --git a/src/threshold.jl b/src/threshold.jl new file mode 100644 index 0000000..206da85 --- /dev/null +++ b/src/threshold.jl @@ -0,0 +1,211 @@ +export cal_mTRS_base, cal_mTRS_full + + +""" +$(TYPEDSIGNATURES) +""" +function cal_mTRS_base3!(Q::AbstractArray, data::AbstractArray{T,3}, mmdd; + dims=3, + probs::Vector=[0.90, 0.95, 0.99, 0.999, 0.9999], + parallel::Bool=true, halfwin::Int=7, use_mov=true, + method_q="base", na_rm=false, + ignore...) where {T<:Real} + + doy_max = maximum(mmdd) + + @inbounds @par parallel for doy = 1:doy_max + ind = filter_mds(mmdd, doy; doy_max, halfwin, use_mov) + + q = @view Q[:, :, doy, :] + x = @view data[:, :, ind] # 这一步耗费内存 + if method_q == "base" + NanQuantile_3d!(q, x; probs, dims, na_rm) + # NanQuantile_low!(q, x; probs, dims=3, na_rm) + elseif method_q == "mapslices" + q .= NanQuantile(x; probs, dims, na_rm) # mapslices is suppressed for 3d `NanQuantile` + end + end + Q +end + + +""" +$(TYPEDSIGNATURES) +""" +function cal_mTRS_base!(Q::AbstractArray{T}, + arr::AbstractArray{T,N}, mmdd; + dims=N, + fun=NanQuantile, + probs::Vector=[0.90, 0.95, 0.99, 0.999, 0.9999], + parallel=true, halfwin=7, use_mov=true, + ignore...) where {T<:Real,N} + + doy_max = maximum(mmdd) + + @inbounds @par parallel for doy = 1:doy_max + ind = filter_mds(mmdd, doy; doy_max, halfwin, use_mov) + # idx = ntuple(i -> i == dims ? ind : Colon(), N) + # ridx = ntuple(i -> i == dims ? doy : Colon(), N + 1) + # x = @view(arr[idx...]) + # q = @view(Q[ridx...]) + x = selectdim(arr, dims, ind) + q = selectdim(Q, dims, doy) + q .= fun(x; probs, dims) # NanQuantile, mapslices + end + Q +end + + + +""" + $(TYPEDSIGNATURES) + +Moving Threshold for Heatwaves Definition + +# Arguments +- `arr` : `time` should be in the last dimension. + +- `type`: The matching type of the moving `doys`, "md" (default) or "doy". + +# Return +- `TRS`: in the dimension of `[nlat, nlon, ndoy, nprob]` + +# References +1. Vogel, M. M., Zscheischler, J., Fischer, E. M., & Seneviratne, S. I. (2020). + Development of Future Heatwaves for Different Hazard Thresholds. Journal of + Geophysical Research: Atmospheres, 125(9). + +""" +function cal_mTRS_base(arr::AbstractArray{T,N}, dates; + dims=N, + use_quantile=true, (fun!)=cal_mTRS_base!, + probs::Vector=[0.90, 0.95, 0.99, 0.999, 0.9999], + dtype=nothing, + p1::Int=1961, p2::Int=1990, kw...) where {T<:Real,N} + + mmdd = factor(format_md.(dates)).refs + doy_max = maximum(mmdd) + + nprob = length(probs) + dtype = dtype === nothing ? eltype(arr) : dtype + + if use_quantile + dims_r = map(d -> d in dims ? [doy_max, nprob] : size(arr, d), 1:N) |> x -> vcat(x...) + else + dims_r = map(d -> d in dims ? doy_max : size(arr, d), 1:N) |> x -> vcat(x...) + end + Q = zeros(dtype, dims_r...) + + # constrain date in [p1, p2] + inds = p1 .<= year.(dates) .<= p2 + _data = selectdim(arr, dims, inds) + _mmdd = @view mmdd[inds] + + fun!(Q, _data, _mmdd; dims, probs, kw...) +end + +cal_mTRS = cal_mTRS_base + + + +""" +$(TYPEDSIGNATURES) + +Moving Threshold for Heatwaves Definition + +# Arguments + +- `use_mov`: Boolean (default true). + + if `true`, 31*15 values will be used to calculate threshold for each grid; + + if `false`, the input `arr` is smoothed first, then only 15 values will be + used to calculate threshold. + +!!! 必须是完整的年份，不然会出错 + +# References +1. Vogel, M. M., Zscheischler, J., Fischer, E. M., & Seneviratne, S. I. (2020). + Development of Future Heatwaves for Different Hazard Thresholds. Journal of + Geophysical Research: Atmospheres, 125(9). + +""" +function cal_mTRS_full(arr::AbstractArray{T,N}, dates; + dims=N, + width=15, verbose=true, + use_quantile=true, (fun!)=cal_mTRS_base!, + use_mov=true, + probs=[0.90, 0.95, 0.99, 0.999, 0.9999], kw...) where {T<:Real,N} + + years = year.(dates) + grps = unique_sort(years) + + YEAR_MIN = minimum(grps) + YEAR_MAX = maximum(grps) + + mmdd = factor(format_md.(dates)).refs + mds = unique_sort(mmdd) + doy_max = length(mds) + + # 滑动平均两种不同的做法 + if !use_mov + printstyled("running: 15d moving average first ... ") + @time arr = movmean(arr, 7; dims, FT=Float32) + end + + if use_quantile + nprob = length(probs) + dim_full = map(d -> d in dims ? [size(arr, d), nprob] : size(arr, d), 1:N) |> x -> vcat(x...) + dim_mTRS = map(d -> d in dims ? [doy_max, nprob] : size(arr, d), 1:N) |> x -> vcat(x...) + else + dim_full = size(arr) + dim_mTRS = map(d -> d in dims ? doy_max : size(arr, d), 1:N) |> x -> vcat(x...) + end + + mTRS_full = zeros(T, dim_full...) + mTRS = zeros(T, dim_mTRS...) + + TRS_head = zeros(T, dim_mTRS...) + TRS_tail = zeros(T, dim_mTRS...) + inds_head = findall(YEAR_MIN .<= years .<= (YEAR_MIN + width * 2)) + inds_tail = findall((YEAR_MAX - width * 2) .<= years .<= YEAR_MAX) + fun!(TRS_head, selectdim(arr, dims, inds_head), @view(mmdd[inds_head]); use_mov, probs, kw...) + fun!(TRS_tail, selectdim(arr, dims, inds_tail), @view(mmdd[inds_tail]); use_mov, probs, kw...) + + for year = grps + verbose && mod(year, 20) == 0 && println("running [year=$year]") + + inds_year = years .== year + md = @view(mmdd[inds_year]) |> unique + inds_md = findall(r_in(mds, md)) # this for mTRS + + year_beg = max(year - width, YEAR_MIN) + year_end = min(year + width, YEAR_MAX) + + if year <= YEAR_MIN + width + _mTRS = TRS_head + elseif year >= YEAR_MAX - width + _mTRS = TRS_tail + else + inds_data = year_beg .<= years .<= year_end + _data = selectdim(arr, dims, inds_data) + _mmdd = @view mmdd[inds_data] + + fun!(mTRS, _data, _mmdd; use_mov, probs, kw...) + _mTRS = mTRS # 366, 后面统一取ind + end + idx_year = ntuple(d -> d in dims ? inds_year : Colon(), N + 1) + idx_md = ntuple(d -> d in dims ? inds_md : Colon(), N + 1) + @views copy!(mTRS_full[idx_year...], _mTRS[idx_md...]) + # 两种方式表现差别不大 + # copy!(selectdim(mTRS_full, dims, inds_year), selectdim(_mTRS, dims, inds_md)) + end + mTRS_full +end + + +cal_pTRS(A; kw...) = NanQuantile(A; kw...) + +function cal_pTRS(A, dates; p1=1981, p2=2010, probs=[0.90, 0.95, 0.99, 0.999, 0.9999], kw...) + inds = p1 .<= year.(dates) .<= p2 + data = selectdim(A, ndims(A), inds) + NanQuantile(data; probs, kw...) +end diff --git a/src/warming_level.jl b/src/warming_level.jl new file mode 100644 index 0000000..38f05d0 --- /dev/null +++ b/src/warming_level.jl @@ -0,0 +1,67 @@ +""" +Calculate yearly air temperature. + +# Description +we use the fixed thresholds and add the seasonal warming signal. Thus, +thresholds are defined as a fixed baseline (such as for the fixed threshold) +plus seasonally moving mean warming of the corresponding future climate based on +the 31-year moving mean of the warmest three months. + +# Details +This function calculates the yearly air temperature based on the input +temperature data and dates. If `only_summer` is true, it only calculates the +temperature for summer months. The function applies the calculation along the +specified dimensions. + +# Arguments +- `A::AbstractArray{T,N}`: input array of temperature data. +- `dates`: array of dates corresponding to the temperature data. +- `dims=N`: dimensions to apply the function along. +- `only_summer=false`: if true, only calculate temperature for summer months. + +# Returns +- `T_year`: array of yearly temperature data. +""" +function cal_yearly_Tair(A::AbstractArray{T,N}, dates; + dims=N, only_summer=false) where {T<:Real,N} + + if only_summer + years = year.(dates) + months = month.(dates) + yms = years .* 100 .+ months + # yms = format.(dates, "yyyy-mm") + T_mon = apply(A, dims; by=yms) + T_mon = movmean(T_mon, 1; dims) #3个月滑动平均 + + ys = fld.(unique(yms), 100) # convert `year * 100 + month` to `year` + T_year = apply(T_mon, dims; by=ys, fun=maximum) # 最热的3个月，作为每年的升温幅度 + else + ys = year.(dates) + T_year = apply(A, dims; by=ys) + end + T_year +end + +cal_climatology_season(A::AbstractArray, dates) = cal_yearly_Tair(A, dates; only_summer=false) + +cal_mTRS_season(A::AbstractArray, dates) = cal_yearly_Tair(A, dates; only_summer=true) + + +""" +$(TYPEDSIGNATURES) +""" +function cal_warming_level(A::AbstractArray{T,N}, dates; + p1=1981, p2=2010, dims=N, only_summer=false) where {T<:Real,N} + + T_year = cal_yearly_Tair(A, dates; dims, only_summer) + + # yms = format.(dates, "yyyy-mm") + # ys = SubString.(unique(yms), 1, 4) + grps = year.(dates) |> unique_sort + inds_clim = @.(p1 <= grps <= p2) + + T_year_clim = selectdim(T_year, dims, inds_clim) + T_clim = apply(T_year_clim, dims; fun=nanmean) + + T_year .- T_clim +end diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..d248427 --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,9 @@ +using Test +using ClimateStats +using NaNStatistics +using Dates +using Ipaper: set_seed, obj_size + +include("test-stat_anomaly.jl") +include("test-stat_threshold.jl") +include("test-stat_warmingLevel.jl") diff --git a/test/test-stat_anomaly.jl b/test/test-stat_anomaly.jl new file mode 100644 index 0000000..39daa0a --- /dev/null +++ b/test/test-stat_anomaly.jl @@ -0,0 +1,120 @@ +# using Ipaper +using Test + +@testset "filter_mds" begin + function check_mds(doy=1) + findall(filter_mds(1:366, doy)) == sort(find_adjacent_doy(doy)) + end + + @test check_mds(1) + @test check_mds(10) + @test check_mds(366) +end + + +@testset "_cal_anomaly_3d" begin + ny = 10 + dates = Date(2010):Day(1):Date(2010 + ny - 1, 12, 31) + ntime = length(dates) + + dims = (100, 100) + set_seed(1) + nprob = () + + set_seed(1) + A = rand(dims..., ntime) + TRS = rand(dims..., 366, nprob...) + T_wl = rand(dims..., ny) + + r1 = _cal_anomaly(A, TRS, dates; T_wl) + r3 = _cal_anomaly_3d(A, TRS, dates; T_wl) + + @time r1 = _cal_anomaly(A, TRS, dates; T_wl) + @time r3 = _cal_anomaly_3d(A, TRS, dates; T_wl) + @test r1 == r3 +end + + +@testset "cal_anomaly_clim 3d" begin + dates = Date(1961, 1, 1):Day(1):Date(2000, 12, 31) |> collect + n = length(dates) + set_seed(1) + + A = rand(Float32, 4, 4, n) + obj_size(A) + + ## 采用`quantile`计算 + kw = (; parallel=true, p1=1961, p2=1980, na_rm=false) + @time anom_pTRS = cal_anomaly_quantile(A, dates; kw..., method="pTRS") + @time anom_base = cal_anomaly_quantile(A, dates; kw..., method="base") + @time anom_season = cal_anomaly_quantile(A, dates; kw..., method="season") + @time anom_full = cal_anomaly_quantile(A, dates; kw..., method="full") + + @test size(anom_pTRS) == size(anom_base) + + ## 采用`fun_clim`:`nanmean`计算 + kw = (; parallel=true, p1=1961, p2=1980, fun_clim=nanmedian) + @time anom_base2 = cal_anomaly_clim(A, dates; kw..., method="base") + @time anom_season2 = cal_anomaly_clim(A, dates; kw..., method="season") + @time anom_full2 = cal_anomaly_clim(A, dates; kw..., method="full") + + @test anom_base == anom_base2 + @test anom_season == anom_season2 + @test anom_full ≈ anom_full2 +end + + + +## Test for multi-dimension array --------------------------------------------- +function test_climatology(; dims=(), T=Float32) + A = rand(T, dims..., ntime) + kw = (; p1=2010, p2=2015, parallel=true, use_mov=true, fun=nanmean) # + + @time r_base = cal_climatology_base(A, dates; kw...) + @time r_full = cal_climatology_full(A, dates; kw...) + + @test size(r_base) == (dims..., 366) + @test size(r_full) == size(A) +end + +function test_anomaly(; dims=(), T=Float32) + A = rand(T, dims..., ntime) + kw = (; p1=2010, p2=2015, parallel=true, use_mov=true, fun_clim=nanmean) + + r_base = cal_anomaly_clim(A, dates; kw..., method="base") + r_seas = cal_anomaly_clim(A, dates; kw..., method="season") + r_full = cal_anomaly_clim(A, dates; kw..., method="full") + + @test size(r_base) == size(A) + @test size(r_seas) == size(A) + @test size(r_full) == size(A) +end + +function test_anomaly_quantile(; T=Float32, dims=(4,)) + A = rand(T, dims..., ntime) + kw = (; parallel=true, p1=2010, p2=2015, na_rm=false, probs=[0.5, 0.9]) + + anom_season = cal_anomaly_quantile(A, dates; kw..., method="season") + anom_base = cal_anomaly_quantile(A, dates; kw..., method="base") + anom_full = cal_anomaly_quantile(A, dates; kw..., method="full") + + @test size(anom_base) == (dims..., ntime, length(kw.probs)) + @test size(anom_base) == size(anom_full) + @test size(anom_base) == size(anom_season) +end + +ny = 10 +dates = Date(2010):Day(1):Date(2010 + ny - 1, 12, 31) +ntime = length(dates) +dims = (4, 4) + +@testset "climatology and anomaly in multiple dimension" begin + l_dims = [(), (4,), (4, 4), (4, 4, 4)] + for T in (Float32, Float64) + for dims = l_dims + test_climatology(; T, dims) + test_anomaly(; T, dims) + test_anomaly_quantile(; T, dims) + end + end +end diff --git a/test/test-stat_threshold.jl b/test/test-stat_threshold.jl new file mode 100644 index 0000000..2fd88d7 --- /dev/null +++ b/test/test-stat_threshold.jl @@ -0,0 +1,94 @@ +@testset "cal_mTRS_base" begin + # import Ipaper: cal_mTRS_base, cal_mTRS_full + dates = Date(1961, 1, 1):Day(1):Date(2000, 12, 31) |> collect + n = length(dates) + set_seed(1) + + A = rand(Float32, 4, 4, n) + + probs = [0.90, 0.95, 0.99, 0.999, 0.9999] + kw = (; probs, na_rm=true, parallel=false, (fun!)=Ipaper.cal_mTRS_base3!) + + @time r1 = cal_mTRS_base(A, dates; kw..., method_q="mapslices") + @time r2 = cal_mTRS_base(A, dates; kw..., method_q="base") + @test r1 == r2 +end + + +@testset "Threshold_nd" begin + + kw = (; parallel=true, p1=1961, p2=1965, na_rm=true) + dates = Date(1961, 1, 1):Day(1):Date(2000, 12, 31) #|> collect + n = length(dates) + set_seed(1) + + A = rand(Float32, 4, 4, n) + + @time r1 = cal_mTRS_base(A, dates; kw..., method_q="mapslices") + @time r3 = cal_mTRS_full(A, dates; kw..., method_q="mapslices") + + @test_nowarn begin + # set_seed(1) + # 1d + A = rand(Float32, n) + @time r2 = cal_mTRS_base(A, dates; kw...) + @time r4 = cal_mTRS_full(A, dates; kw...) + + # 2d + A = rand(Float32, 4, n) + @time r2 = cal_mTRS_base(A, dates; kw...) + @time r4 = cal_mTRS_full(A, dates; kw...) + + # 3d + A = rand(Float32, 4, 4, n) + @time r2 = cal_mTRS_base(A, dates; kw...) + @time r4 = cal_mTRS_full(A, dates; kw...) + + # 4d + A = rand(Float32, 4, 4, 4, n) + @time r2 = cal_mTRS_base(A, dates; kw...) + @time r4 = cal_mTRS_full(A, dates; kw...) + end +end + + +@testset "mTRS and climatology" begin + # p1 = 1961, p2 = 1965, + kw = (; parallel=false, na_rm=true) + dates = Date(1961, 1, 1):Day(1):Date(2000, 12, 31) #|> collect + n = length(dates) + set_seed(1) + + A = rand(Float32, 4, 4, n) + + # mTRS_base + @time r1 = cal_climatology_base(A, dates; fun=nanmedian, kw...) + @time r2 = cal_mTRS_base(A, dates; probs=[0.5], kw...)[:, :, :, 1] + # @time r3 = cal_mTRS_base(A, dates; + # use_quantile=false, (fun!)=Ipaper.cal_climatology_base3!, fun=nanmedian, kw...) + @test r1 == r2 + # @test r2 == r3 + + # mTRS_full + @time r1 = cal_climatology_full(A, dates; fun=nanmedian, kw...) + @time r2 = cal_mTRS_full(A, dates; probs=[0.5], kw...)[:, :, :, 1] + # @time r3 = cal_mTRS_full(A, dates; + # use_quantile=false, (fun!)=Ipaper.cal_climatology_base3!, fun=nanmedian, kw...) + + @test r1 == r2 + # @test r2 == r3 +end + + +@testset "cal_threshold" begin + dates = Date(2010, 1):Day(1):Date(2020, 12, 31) + ntime = length(dates) + data = rand(10, ntime) + r_full = cal_threshold(data, dates; p1=2010, p2=2015, method="full") + r_base = cal_threshold(data, dates; p1=2010, p2=2015, method="base") + r_seas = cal_threshold(data, dates; p1=2010, p2=2015, method="season") + + @test size(r_full) == (10, ntime, 5) + @test size(r_base) == (10, 366, 5) + @test size(r_seas) == (10, 11) +end diff --git a/test/test-stat_warmingLevel.jl b/test/test-stat_warmingLevel.jl new file mode 100644 index 0000000..47ed51d --- /dev/null +++ b/test/test-stat_warmingLevel.jl @@ -0,0 +1,17 @@ +import ClimateStats: cal_mTRS_season, cal_climatology_season + + +@testset "warmLevel" begin + dates = Date(1961):Day(1):Date(2000, 12, 31) |> collect + n = length(dates) + set_seed(1) + + arr = rand(Float32, 4, 4, n) + obj_size(arr) + + @test_nowarn T1 = cal_warming_level(arr, dates; only_summer=true, p1=1961, p2=2000) + @test_nowarn T2 = cal_warming_level(arr, dates; only_summer=false, p1=1961, p2=2000) + + # r1 =cal_mTRS_season(arr, dates) + # r2 = cal_climatology_season(arr, dates) +end