add function processing MODIS

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

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

src/dates.jl → 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

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

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")

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

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

src/sf/sf.jl

@@ -9,6 +9,7 @@ const guanshan_flowdir_gis = "$dir_data/GuanShan_flowdir_gis.tif"
 export guanshan_dem, guanshan_flowdir_cpp, guanshan_flowdir_gis
 
 
+using ProgressMeter
 using Ipaper: file_ext, obj_size
 
 export bbox, in_bbox, bbox_overlap