bilinear passed test

Project.toml

@@ -10,7 +10,6 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 Ipaper = "e58298cb-69f7-4186-aecd-5834b6793426"
-LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
@@ -32,8 +31,7 @@ RbaseTerraExt = "RCall"
 CFTime = "0.1"
 DataFrames = "1.6"
 DocStringExtensions = "0.9"
-Ipaper = "0.1.16, 0.1.17, 0.2"
-LoopVectorization = "0.12"
+Ipaper = "0.2"
 NCDatasets = "0.13, 0.14"
 Parameters = "0.12"
 PrecompileTools = "1.2"

src/BiasCorrection/QDM.jl

@@ -62,7 +62,7 @@ function QDM_chunk!(y_adj::AbstractVector{T},
   y_obs::AbstractVector{T}, y_calib::AbstractVector{T}, y_pred::AbstractVector{T}, dates; 
   ny_win=10, type_adj="add") where {T<:Real}
 
-  lst_index = split_date(dates; ny_win, merge_small=0.7)
+  lst_index = split_date(dates; ny_win, ratio_small=0.7)
   for inds in lst_index
     y_adj[inds] .= QDM(y_obs, y_calib, y_pred[inds]; type_adj)
   end

src/Interpolation/bilinear.jl

@@ -1,5 +1,7 @@
 """
-    bilinear(x, y, z::AbstractArray{T,3}, xx, yy; na_rm=true) where {T<:Real}
+    bilinear(x, y, z::AbstractArray{T,3}, xx, yy; na_rm=true, parallel=true, progress=true) where {T<:Real}
+    bilinear(x, y, z::AbstractArray{T,3}; b::bbox, reverse_lat=true, cellsize=1, na_rm=true)
+    bilinear(ra::SpatRaster{T,3}; cellsize=1, na_rm=true, kw...) where {T<:Real}
 
 Suppose that the location, (locx, locy) lies in between the first two grid points in both x an y. That is locx is between x1 and x2 and locy is between y1 and y2. Let `ex= (l1-x1)/(x2-x1)` `ey= (l2-y1)/(y2-y1)`. The interpolant is
 
@@ -29,9 +31,9 @@ Z = rand(T, length(lon), length(lat), 2)
 r = bilinear(lon, lat, Z, Lon, Lat; na_rm=true)
 ```
 """
-function bilinear(x, y, z::AbstractArray{T,3}, xx, yy; na_rm=true, 
+function bilinear(x, y, z::AbstractArray{T,3}, xx, yy; na_rm=true,
   parallel=true, progress=true) where {T<:Real}
-  
+
   nx = length(x)
   ny = length(y)
 
@@ -76,7 +78,7 @@ function bilinear(x, y, z::AbstractArray{T,3}, xx, yy; na_rm=true,
         z12 = z[I, J2, k]
         z21 = z[I2, J, k]
         z22 = z[I2, J2, k]
-        
+
         if na_rm
           zmean = nanmean4(z11, z12, z21, z22) # mean
           isnan(z11) && (z11 = z12)
@@ -89,8 +91,8 @@ function bilinear(x, y, z::AbstractArray{T,3}, xx, yy; na_rm=true,
           isnan(z21) && (z21 = zmean)
           isnan(z22) && (z22 = zmean)
         end
-        
-        @fastmath res[i, j, k] = 
+
+        @fastmath res[i, j, k] =
           z11 * (1 - _ex) * (1 - _ey) +
           z12 * (1 - _ex) * _ey +
           z21 * _ex * (1 - _ey) +
@@ -102,15 +104,17 @@ function bilinear(x, y, z::AbstractArray{T,3}, xx, yy; na_rm=true,
 end
 
 function bilinear(x, y, z::AbstractArray{T,3};
-  range=[70, 140, 15, 55], cellsize=1, na_rm=true) where {T<:Real}
-
-  delta = cellsize / 2
-  rlon = range[1:2]
-  rlat = range[3:4]
-  xx = rlon[1]+delta:cellsize:rlon[2] #|> collect
-  yy = rlat[1]+delta:cellsize:rlat[2] #|> collect
+  b::bbox, reverse_lat=true, cellsize=1, na_rm=true) where {T<:Real}
+  xx, yy = bbox2dims(b; cellsize, reverse_lat)
   bilinear(x, y, z, xx, yy; na_rm)
 end
 
+function bilinear(ra::SpatRaster{T,3}; cellsize=1, na_rm=true, kw...) where {T<:Real}
+  (; lon, lat, b, time, name, bands) = ra
+  xx, yy = bbox2dims(b; cellsize)
+  Z = bilinear(lon, lat, ra.A, xx, yy; na_rm)
+  rast(Z, b; time, name, bands, kw...)
+end
+
 
 export bilinear

src/Interpolation/utilize.jl

@@ -1,13 +1,11 @@
 export rm_empty, weighted_nanmean!, weighted_nanmean, earth_dist
-using LoopVectorization
 
 function weighted_nanmean(x::AbstractVector{T1}, w::AbstractVector{T2}) where {T1,T2}
   T = promote_type(T1, T2)
   ∑ = ∅ = T(0)
   ∑w = ∅w = T2(0)
 
-  # @inbounds 
-  @turbo for i = eachindex(x)
+  @inbounds for i = eachindex(x)
     # if !isnan(x[i]); end
     xᵢ = x[i]
     notnan = xᵢ == xᵢ

src/nc_write.jl

@@ -78,7 +78,11 @@ $(METHODLIST)
 """
 function nc_write!(f::AbstractString, varname::AbstractString, val,
   dims::Vector{<:Union{NcDim,AbstractString}}, attrib::Dict=Dict();
+  units=nothing, longname=nothing,
   compress=1, kw...)
+
+  !isnothing(units) && (attrib["units"] = units)
+  !isnothing(longname) && (attrib["longname"] = longname)
 
   mode = check_file(f) ? "a" : "c"
   ds = nc_open(f, mode)

src/tools.jl

@@ -1,22 +1,29 @@
 # https://github.com/pacificclimate/ClimDown/blob/master/R/QDM.R
 # mk.factor.set
-function split_chunk(n::Int, nchunk=4; chunk = nothing, merge_small=0.5)
+"""
+    split_chunk(n::Int, nchunk=4; chunk=nothing, ratio_small::Real=0.5)
+
+# Arguments
+- `ratio_small`: The minimum ratio of the last chunk to the maximum length of chunks.
+如果最后一个块的小于该比例，则将最后两个块合并。
+"""
+function split_chunk(n::Int, nchunk=4; chunk=nothing, ratio_small::Real=0.0)
   if chunk === nothing
     chunk = cld(n, nchunk)
   else
     nchunk = cld(n, chunk)
   end
-  
+
   lst = map(i -> begin
-    i_beg = (i - 1) * chunk + 1
-    i_end = min(i * chunk, n)
-    i_beg:i_end
-  end, 1:nchunk)
-  
+      i_beg = (i - 1) * chunk + 1
+      i_end = min(i * chunk, n)
+      i_beg:i_end
+    end, 1:nchunk)
+
   lens = map(length, lst)
-  len_max = maximum(lens)
+  nmin = ratio_small * chunk
 
-  if lens[end] < merge_small * len_max
+  if lens[end] < nmin
     i_beg = lst[end-1][1]
     i_end = lst[end][end]
     lst[end-1] = i_beg:i_end
@@ -25,12 +32,11 @@ function split_chunk(n::Int, nchunk=4; chunk = nothing, merge_small=0.5)
   lst
 end
 
-function split_chunk(x::Union{UnitRange, AbstractVector}, nchunk=4; kw...)
+function split_chunk(x::Union{UnitRange,AbstractVector}, nchunk=4; kw...)
   lst = split_chunk(length(x), nchunk; kw...)
   map(i -> x[i], lst)
 end
 
-
 # dates = make_date(2015, 1, 1):Day(1):make_date(2100, 12, 31) |> collect
 # - chunk: nyear
 function split_date(dates; ny_win=10, kw...)

test/interp/test-bilinear.jl

@@ -1,5 +1,5 @@
 # major update, fix the error of bilinear nanmean4
-using Test
+using Test, NetCDFTools
 
 @testset "nanmean" begin
   x = [1.0, 2, 3, 4]
@@ -9,23 +9,28 @@ end
 
 
 @testset "bilinear" begin
-  cell = 10
-  lon = 70:cell:140
-  lat = 15:cell:55
-  ntime = 2
-
+  set_seed(1)
+  b = bbox(70, 15, 140, 55)
+  lon, lat = bbox2dims(b; cellsize=10)
+
   cell2 = 5
-  Lon = 70+cell2/2:cell2:140
-  Lat = 15+cell2/2:cell2:55
+  ntime = 2
   Z = rand(Float32, length(lon), length(lat), ntime)
-  r1 = bilinear(lon, lat, Z, Lon, Lat; na_rm=true)
-  r2 = bilinear(lon, lat, Z; range=[70, 140, 15, 55], cellsize=cell2)
+
+  dates = Date(2010):Day(1):Date(2010, 1, 1) |> collect
+  ra = rast(deepcopy(Z), b; time=dates, name="tasmax")
+  ra_5 = bilinear(ra; cellsize=5)
 
+  Lon, Lat = bbox2dims(b; cellsize=5)
+  r1 = bilinear(lon, lat, Z, Lon, Lat; na_rm=true)
+  r2 = bilinear(lon, lat, Z; b, cellsize=cell2)
+
   @test size(r1) == (length(Lon), length(Lat), ntime)
   @test r1 == r2
-  # 如何检测结果是否正确？
+  @test ra_5.A == r1
 end
 
+
 @testset "bilinear vs. cdo" begin
   f_cdo = proj_path("data/HI_tasmax_resampled_cdo.nc")
   # cdo_bilinear(f, fout, fgrid; verbose=true)

test/runtests.jl

@@ -1,13 +1,15 @@
 using Test
 using NetCDFTools
-using Ipaper
+using Ipaper, Dates
+import NetCDFTools: weighted_nanmean
 
 dir_root = dirname(dirname(@__FILE__))
 proj_path(f) = dirname(dirname(@__FILE__)) * "/" *  f
 # println(dirname(@__FILE__))
 # println(pwd())
-include("interp/test-angle.jl")
+include("test-Ipaper.jl")
 include("interp/test-weighted_nanmean.jl")
+include("interp/test-angle.jl")
 include("interp/test-spInterp.jl")
 include("interp/test-bilinear.jl")
 

test/test-Ipaper.jl

@@ -0,0 +1,25 @@
+using Test, NetCDFTools, Ipaper
+
+@testset "updateMask" begin
+  A = array(1.0:16, (4, 4)) |> collect
+  mask = A .> 10
+  updateMask!(A, mask)
+  length(findall(A .> 10)) == 6
+
+  A = rand(4, 4, 2)
+  mask = A[:, :, 1] .> 0.5
+  updateMask!(A, mask)
+  @test isempty(findall(A[:, :, 1] .<= 0.5))
+end
+
+
+@testset "split_chunk" begin
+  @test split_chunk(1:10, 4; ratio_small=0.5) == [1:3, 4:6, 7:10]
+  @test split_chunk(1:10, 4; ratio_small=0.0) == [1:3, 4:6, 7:9, 10:10]
+
+  dates = make_date(2015, 1, 1):Day(1):make_date(2100, 12, 31) #|> collect
+  lst = split_date(dates; ny_win=10, ratio_small=0.0)
+  @test length(lst) == 9
+  lst = split_date(dates; ny_win=10, ratio_small=0.7)
+  @test length(lst) == 8
+end