From 4d5c91587e9cdca260518fb30398bf887f33b9f8 Mon Sep 17 00:00:00 2001 From: Dongdong Kong Date: Thu, 3 Oct 2024 12:48:45 +0800 Subject: [PATCH] need to test bilinear --- Project.toml | 2 +- src/Interpolation/bilinear.jl | 30 +++++++++++++++++------------- src/nc_write.jl | 4 ++++ test/interp/test-bilinear.jl | 27 ++++++++++++++++----------- test/runtests.jl | 2 +- 5 files changed, 39 insertions(+), 26 deletions(-) diff --git a/Project.toml b/Project.toml index 74dbe34..7e94c15 100644 --- a/Project.toml +++ b/Project.toml @@ -32,7 +32,7 @@ RbaseTerraExt = "RCall" CFTime = "0.1" DataFrames = "1.6" DocStringExtensions = "0.9" -Ipaper = "0.1.16, 0.1.17, 0.2" +Ipaper = "0.2" LoopVectorization = "0.12" NCDatasets = "0.13, 0.14" Parameters = "0.12" diff --git a/src/Interpolation/bilinear.jl b/src/Interpolation/bilinear.jl index 7d3c4d4..3d846f8 100644 --- a/src/Interpolation/bilinear.jl +++ b/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 diff --git a/src/nc_write.jl b/src/nc_write.jl index 6af1cd9..8e77817 100644 --- a/src/nc_write.jl +++ b/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) diff --git a/test/interp/test-bilinear.jl b/test/interp/test-bilinear.jl index 784fd98..480ffd8 100644 --- a/test/interp/test-bilinear.jl +++ b/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) diff --git a/test/runtests.jl b/test/runtests.jl index 08cc3fc..05a3248 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,6 @@ using Test using NetCDFTools -using Ipaper +using Ipaper, Dates dir_root = dirname(dirname(@__FILE__)) proj_path(f) = dirname(dirname(@__FILE__)) * "/" * f