generalise routines for VanWessem2018

scripts/vanwessem2018/main-looped.jl

@@ -0,0 +1,64 @@
+include("../intro.jl")
+
+filepaths = String[]
+for (root, dirs, files) in walkdir(datadir("BCs/VanWessem-2018/3D"))
+    for file in files
+        if occursin(".nc", file)
+            push!(filepaths, joinpath(root, file))
+        end
+    end
+end
+
+varnames = String[]
+for filepath in filepaths
+    i1 = findlast("/", filepath)[1]+1
+    i2 = findfirst("_", filepath)[1]-1
+    push!(varnames, lowercase(filepath[i1:i2]))
+end
+varnames = Tuple(varnames)
+
+# Define regridding config
+source_dimnames = ("rlon", "rlat", "time")
+target_dimnames = ("xc", "yc", "time")
+source_gridname = "+proj=stere +lat_0=-90 +lat_ts=-80"
+target_gridname = "+proj=stere +lat_0=-90 +lat_ts=-80"
+extrapolation_boundary_conditions = (Flat(), Flat(), Flat())
+
+aggregate_over_months(X) = cat([mean(X[:, :, i:12:end], dims = 3) for i in 1:12]..., dims = 3)
+aggregate_dims((x1, x2, x3)) = return(x1, x2, collect(1:12))
+
+# Regrid
+x = range(-3040f3, stop = 3040f3, step = 32f3)
+y = copy(x)
+t = collect(1:12)
+target_dims = (x, y, t)
+target_grid = ndgrid(target_dims...)
+target_vars = Vector{Array{Float32, 3}}(undef, length(filepaths))
+target_atts = Vector{Dict{String, String}}(undef, length(filepaths))
+for i in eachindex(filepaths)
+    regrid = Regrid(
+        filepaths[i],
+        source_dimnames,
+        target_dimnames,
+        source_gridname,
+        target_gridname,
+        varnames[i:i],
+        extrapolation_boundary_conditions;
+        scale_dims_by = 1f5,
+        aggregate_vars = aggregate_over_months,
+        aggregate_dims = aggregate_dims,
+        attributes2extract = ("units", "long_name")
+    )
+    target_vars[i] = regrid(target_grid)[1]
+    target_atts[i] = regrid.attributes[1]
+end
+
+# heatmap(target_vars[4][:, :, 7])
+
+# Save regridded data
+x_atts = Dict("units" => "m", "long_name" => "x-coordinate")
+y_atts = Dict("units" => "m", "long_name" => "y-coordinate")
+t_atts = Dict("units" => "month", "long_name" => "Month number")
+dim_atts = (x_atts, y_atts, t_atts)
+fn = datadir("ANT-32KM_RACMO-VW18.nc")
+save2nc(fn, target_dimnames, target_dims ./ 1f3, dim_atts, varnames, target_vars, target_atts)

scripts/vanwessem2018/main-scalar.jl

@@ -0,0 +1,54 @@
+include("../intro.jl")
+
+filepaths = String[]
+for (root, dirs, files) in walkdir(datadir("BCs/VanWessem-2018"))
+    for file in files
+        if occursin(".nc", file)
+            push!(filepaths, joinpath(root, file))
+        end
+    end
+end
+
+varnames = String[]
+for filepath in filepaths
+    i1 = findlast("/", filepath)[1]+1
+    i2 = findfirst("_", filepath)[1]-1
+    push!(varnames, lowercase(filepath[i1:i2]))
+end
+varnames = Tuple(varnames)
+
+# Define regridding config
+source_dimnames = ("rlon", "rlat", "time")
+target_dimnames = ("xc", "yc", "time")
+source_gridname = "+proj=stere +lat_0=-90 +lat_ts=-80"
+target_gridname = "+proj=stere +lat_0=-90 +lat_ts=-80"
+extrapolation_boundary_conditions = (Flat(), Flat(), Flat())
+
+aggregate_over_months(X) = cat([mean(X[:, :, i:12:end], dims = 3) for i in 1:12]..., dims = 3)
+aggregate_dims((x1, x2, x3)) = return(x1, x2, collect(1:12))
+
+i = 2
+# Define mean GHF regridding
+regrid = Regrid(
+    filepaths[i],
+    source_dimnames,
+    target_dimnames,
+    source_gridname,
+    target_gridname,
+    varnames[i:i],
+    extrapolation_boundary_conditions;
+    scale_dims_by = 1f5,
+    aggregate_vars = aggregate_over_months,
+    aggregate_dims = aggregate_dims,
+)
+
+# Regrid
+x = range(-3040f3, stop = 3040f3, step = 32f3)
+y = copy(x)
+t = collect(1:12)
+target_dims = (x, y, t)
+target_grid = ndgrid(target_dims...)
+
+# Filter out missing values and pack into a tuple.
+target_var = regrid(target_grid)[1]
+heatmap(target_var[:, :, 2])

src/regrid.jl

@@ -56,6 +56,7 @@ struct Regrid
     vars
     extrapolation_boundary_conditions#::NTuple{N, <:Interpolations.BoundaryCondition}
     interpolators
+    attributes
 end
 
 function Regrid(
@@ -65,11 +66,23 @@ function Regrid(
     source_gridname,
     target_gridname,
     varnames,
-    extrapolation_boundary_conditions,
+    extrapolation_boundary_conditions;
+    scale_dims_by = 1,
+    aggregate_vars = nothing,
+    aggregate_dims = nothing,
+    attributes2extract = nothing,
 )
     sanitycheck_lon_lat(source_dimnames)
     source2target, target2source = get_projections(source_gridname, target_gridname)
     source_dims, vars = load_data(source_file, source_dimnames, varnames)
+    attributes = get_attributes(source_file, varnames, attributes2extract)
+
+    source_dims = scale_dims_by .* source_dims
+    if !isnothing(aggregate_vars)
+        vars = aggregate_vars.(vars)
+        source_dims = aggregate_dims(source_dims)
+    end
+
     interpolators = [linear_interpolation(source_dims, var,
         extrapolation_bc = extrapolation_boundary_conditions) for var in vars]
     return Regrid(
@@ -85,6 +98,7 @@ function Regrid(
         vars,
         extrapolation_boundary_conditions,
         interpolators,
+        attributes,
     )
 end
 

src/utils.jl

@@ -13,10 +13,32 @@ respectively referred to as `dimnames` and `varnames`.
 """
 function load_data(filepath, dimnames, varnames)
     dims = Tuple(ncread(filepath, dimname) for dimname in dimnames)
-    vars = Tuple(ncread(filepath, varname) for varname in varnames)
+    vars = Tuple(dropped_ncread(filepath, varname) for varname in varnames)
     return dims, vars
 end
 
+function dropped_ncread(filepath, varname)
+    var = ncread(filepath, varname)
+    dimnumbers = Tuple(1:length(size(var)))
+
+    if 1 ∈ size(var)
+        return dropdims(var, dims = Tuple([dimnumbers[i] for i in eachindex(dimnumbers)
+            if size(var, i) == 1]))
+    else
+        return var
+    end
+end
+
+function get_attributes(source_file, varnames, attributes2extract)
+    if isnothing(attributes2extract)
+        return [Dict() for varname in varnames]
+    else
+        return Tuple([Dict(att => ncgetatt(source_file, varname, att) for
+            att in attributes2extract) for varname in varnames])
+    end
+end
+
+
 """
     missings2nans(X::Array{Union{T, Missing}})