From c51c2311c870fcb40b82450e4b166ea0864ecad0 Mon Sep 17 00:00:00 2001 From: Dongdong Kong Date: Wed, 2 Oct 2024 12:16:34 +0800 Subject: [PATCH] add function processing MODIS --- src/Dates/Dates.jl | 15 +++++++++++++++ src/Dates/interval_intersect.jl | 21 +++++++++++++++++++++ src/{dates.jl => Dates/utilize.jl} | 16 ++++++---------- src/Dates/weight_d8mon.jl | 30 ++++++++++++++++++++++++++++++ src/Ipaper.jl | 2 +- src/Statistics/Statistics.jl | 5 +++-- src/Statistics/weighted_nansum.jl | 24 ++++++++++++++++++++++++ 7 files changed, 100 insertions(+), 13 deletions(-) create mode 100644 src/Dates/Dates.jl create mode 100644 src/Dates/interval_intersect.jl rename src/{dates.jl => Dates/utilize.jl} (82%) create mode 100644 src/Dates/weight_d8mon.jl create mode 100644 src/Statistics/weighted_nansum.jl diff --git a/src/Dates/Dates.jl b/src/Dates/Dates.jl new file mode 100644 index 0000000..948e0b2 --- /dev/null +++ b/src/Dates/Dates.jl @@ -0,0 +1,15 @@ +import Dates +import Dates: Date, DateTime, Year, Month, Day, year, month, day, format +# using CFTime + +include("interval_intersect.jl") +include("utilize.jl") +include("weight_d8mon.jl") + +export dates_miss, dates_nmiss, + DateTime, Date, year, month, day, Year, Month, Day, format, + make_datetime, make_date, + date_doy, + date_year, date_ym, date_dn +export weight_d8mon +export interval_intersect, datediff diff --git a/src/Dates/interval_intersect.jl b/src/Dates/interval_intersect.jl new file mode 100644 index 0000000..bac8b51 --- /dev/null +++ b/src/Dates/interval_intersect.jl @@ -0,0 +1,21 @@ +function interval_intersect(x::Tuple{T,T}, y::Tuple{T,T}) where {T<:Real} + left = max(x[1], y[1]) # intersect + right = min(x[2], y[2]) + + if left <= right + return right - left + else + return 0 + end +end + +function interval_intersect(x::Tuple{T,T}, y::Tuple{T,T}) where {T<:Union{Date,DateTime}} + left = max(x[1], y[1]) # intersect + right = min(x[2], y[2]) + + if left <= right + return convert(Dates.Day, right - left) + else + return Day(0) + end +end diff --git a/src/dates.jl b/src/Dates/utilize.jl similarity index 82% rename from src/dates.jl rename to src/Dates/utilize.jl index 2b5c6b9..e25d047 100644 --- a/src/dates.jl +++ b/src/Dates/utilize.jl @@ -1,6 +1,9 @@ -import Dates -import Dates: Date, DateTime, Year, Month, Day, year, month, day, format -# using CFTime +add_d8(date::Union{Date,DateTime}) = + min(date + Day(7), DateTime(year(date), 12, 31)) + +function datediff(x::T, y::T; unit=Day) where {T<:Union{Date,DateTime}} + convert(unit, x - y).value +end # only for daily scale function dates_miss(dates) @@ -57,10 +60,3 @@ end Dates.year(x::AbstractString) = parse(Int, x[1:4]) Dates.month(x::AbstractString) = parse(Int, x[6:7]) Dates.day(x::AbstractString) = parse(Int, x[9:10]) - - -export dates_miss, dates_nmiss, - DateTime, Date, year, month, day, Year, Month, Day, format, - make_datetime, make_date, - date_doy, - date_year, date_ym, date_dn diff --git a/src/Dates/weight_d8mon.jl b/src/Dates/weight_d8mon.jl new file mode 100644 index 0000000..8c57677 --- /dev/null +++ b/src/Dates/weight_d8mon.jl @@ -0,0 +1,30 @@ +""" + weight_d8_mon(dates_beg::Vector{T}, date::T) where {T<:Union{Date,DateTime}} + weight_d8_mon(dates_beg::Vector{T}, dates_end::Vector{T}, date::T) where {T<:Union{Date,DateTime}} + +convert MODIS 8-day to monthly +""" +function weight_d8mon(dates_beg::Vector{T}, date::T) where {T<:Union{Date,DateTime}} + dates_end = add_d8.(dates_beg) + weight_d8mon(dates_beg, dates_end, date) +end + +function weight_d8mon(dates_beg::Vector{T}, dates_end::Vector{T}, date::T) where {T<:Union{Date,DateTime}} + + date_beg = date + date_end = DateTime(year(date_beg), month(date_beg), daysinmonth(date_beg)) + interval = (date_beg, date_end) + + inds = findall(@.( + date_beg <= dates_beg <= date_end || + date_beg <= dates_end <= date_end)) + days_full = datediff.(dates_end[inds], dates_beg[inds]) .+ 1 + + days = map(i -> begin + date_beg, date_end = dates_beg[i], dates_end[i] + int2 = date_beg, date_end + interval_intersect(interval, int2) + Day(1) |> x -> x.value + end, inds) + (; date_beg=dates_beg[inds], date_end=dates_end[inds], index=inds, + days, days_full, w=days ./ days_full) +end diff --git a/src/Ipaper.jl b/src/Ipaper.jl index 8dfb937..37a13b5 100644 --- a/src/Ipaper.jl +++ b/src/Ipaper.jl @@ -35,7 +35,7 @@ include("cmd.jl") include("aria2c.jl") include("main_cdo.jl") -include("dates.jl") +include("Dates/Dates.jl") include("factor.jl") include("list.jl") include("par.jl") diff --git a/src/Statistics/Statistics.jl b/src/Statistics/Statistics.jl index 0e17704..e716e46 100644 --- a/src/Statistics/Statistics.jl +++ b/src/Statistics/Statistics.jl @@ -5,15 +5,16 @@ using Statistics: mean, median, quantile # weighted_mean(x, w) = mean(x, weights(w)) # weighted_sum(x, w) = sum(x, weights(w)) -weighted_sum(x::AbstractVector, w::AbstractVector) = sum(x .* w) weighted_mean(x::AbstractVector, w::AbstractVector) = sum(x .* w) / sum(w) include("movmean.jl") include("NanQuantile.jl") include("match2.jl") +include("weighted_nansum.jl") - +export weighhted_nansum export weighted_mean, weighted_sum + export mean, median, quantile, movmean, weighted_movmean export _nanquantile!, NanQuantile diff --git a/src/Statistics/weighted_nansum.jl b/src/Statistics/weighted_nansum.jl new file mode 100644 index 0000000..34e4140 --- /dev/null +++ b/src/Statistics/weighted_nansum.jl @@ -0,0 +1,24 @@ +weighted_sum(x::AbstractVector, w::AbstractVector) = sum(x .* w) + +function weighted_nansum(x::AbstractVector{T}, w::AbstractVector) where {T<:Real} + ∑ = T(0) + @inbounds for i in eachindex(x) + ∑ += ifelse(x[i] == x[i], x[i] * w[i], T(0)) + end + return ∑ +end + +function weighted_nansum(A::AbstractArray{T,3}, w::AbstractVector) where {T<:Real} + nlon, nlat, ntime = size(A) + R = zeros(nlon, nlat) #.* T(NaN) + + @inbounds for i = 1:nlon, j = 1:nlat + # ∑ = T(0) + for k = 1:ntime + xi = A[i, j, k] + xi == xi && (R[i, j] += xi * w[i]) + # R[i, j] += xi == xi ? xi * w[i] : T(0) + end + end + return R +end