Add conversion to zarr for variables at egdes

src/Particles.jl

@@ -46,7 +46,7 @@ export CartesianXYGrid, dump, in_bbox, find_index, find_index_and_weights, apply
 export WmsServer, get_map, plot_image
 
 # dflow.jl
-export load_nc_info, load_dflow_grid, load_nc_var, load_nc_map_slice, find_index, apply_index, get_times, get_reftime, as_DateTime, initialize_interpolation
+export load_nc_info, load_dflow_grid, load_nc_var, load_nc_map_slice, load_nc_map_slice_at_faces, find_index, apply_index, get_times, get_reftime, as_DateTime, initialize_interpolation
 
 # era5.jl
 export EraData, load_map_slice, get_reftime, get_times, as_DateTime, initialize_interpolation

src/dflow.jl

@@ -8,6 +8,34 @@ using Glob
 
 debuglevel = 1 # 0-nothing, larger more output
 
+# cache arrays
+global edge_faces_cache=[]
+global edge_type_cache=[]
+global nfaces_cache=[]
+
+function fill_cache_arrays(map)
+   #global edge_faces_cache
+   #global edge_type_cache
+   if length(edge_faces_cache)==0
+      for i in eachindex(map)
+         temp_mat=map[i].vars["mesh2d_edge_faces"][:,:]
+         push!(edge_faces_cache,temp_mat)
+      end
+   end
+   if length(edge_type_cache)==0
+      for i in eachindex(map)
+         temp_vec=map[i].vars["mesh2d_edge_type"][:]
+         push!(edge_type_cache,temp_vec)
+      end
+   end
+   if length(nfaces_cache)==0
+      for i in eachindex(map)
+         temp=Int64(map[i].dim["mesh2d_nFaces"].dimlen)
+         push!(nfaces_cache,temp)
+      end
+   end
+end
+
 
 """
    map=load_nc_info(path,filename)
@@ -93,17 +121,18 @@ end
    waterlevel = load_nc_map_slice(map,"s1",1)
 Load data for a time-dependent variable for a specific time and for all domains.
 """
-function load_nc_map_slice(map, varname, itime, sigma_layer_index=-1, domainnr = 0)
+function load_nc_map_slice(map, varname, itime, layer_index=-1, domainnr = 0)
    result = []
    domainnr != 0 || (domainnr = 1:length(map))
    for i in domainnr
       # push!(result,map[i][varname][:,itime])
       if ndims(map[i][varname]) == 2
          push!(result, map[i][varname][:,itime])
       elseif ndims(map[i][varname]) == 3
-         if sigma_layer_index > 0
+         if layer_index > 0
             # if sigma_layer_index <= size(map[i][varname])[1]
-               push!(result, map[i][varname][sigma_layer_index,:,itime])
+               @show size(map[i][varname])
+               push!(result, map[i][varname][layer_index,:,itime])
             # else
             #    throw(BoundsError(map[i][varname], sigma_layer_index))
             # end
@@ -118,6 +147,86 @@ function load_nc_map_slice(map, varname, itime, sigma_layer_index=-1, domainnr =
    return result
 end
 
+"""
+function map_edges_to_faces_3d!(at_faces,at_edges,edge_faces,edge_type)
+Average edges around each face to get values at faces
+The function works in place, so at_faces is modified.
+edge_faces : array of size (2,nedges) with the face numbers at the ends of each edge
+edge_type  : array of size (nedges) with the type of each edge (1=internal, ...)
+"""
+function map_edges_to_faces_3d!(at_faces,at_edges,edge_faces,edge_type)
+   # at_faces: array of size (nfaces,nlayers)
+   # at_edges: array of size (nedges,nlayers)
+   nfaces=size(at_faces,1)
+   nlayers=size(at_faces,2)
+   nedges=size(at_edges,1)
+   edge_count=zeros(nfaces) # number of edges per face
+   at_faces[:,:].=0.0 #initialize
+   for i=1:size(edge_faces,2)
+       if edge_type[i]==1 # internal edge
+           for j=1:2
+               if edge_faces[j,i]>0
+                   at_faces[edge_faces[j,i],:]+=at_edges[i,:]
+                   edge_count[edge_faces[j,i]]+=1
+               end
+           end
+       end
+   end
+   for i=1:nfaces
+       if edge_count[i]>0
+           at_faces[i,:]./=edge_count[i]
+       end
+   end
+end
+
+function map_edges_to_faces_2d!(at_faces,at_edges,edge_faces,edge_type)
+   # at_faces: array of size (nfaces)
+   # at_edges: array of size (nedges)
+   nfaces=size(at_faces,1)
+   nedges=size(at_edges,1)
+   edge_count=zeros(nfaces) # number of edges per face
+   at_faces[:].=0.0 #initialize
+   for i=1:size(edge_faces,2)
+       if edge_type[i]==1 # internal edge
+           for j=1:2
+               if edge_faces[j,i]>0
+                   at_faces[edge_faces[j,i]]+=at_edges[i]
+                   edge_count[edge_faces[j,i]]+=1
+               end
+           end
+       end
+   end
+   for i=1:nfaces
+       if edge_count[i]>0
+           at_faces[i]/=edge_count[i]
+       end
+   end
+end
+
+"""
+   waterlevel = load_nc_map_slice_at_faces(map,"vicwwu",1,5)
+Load data for a time-dependent variable for a specific time and for all domains.
+Apply simple interpolation from edges to faces if needed.
+"""
+function load_nc_map_slice_at_faces(map, varname, itime, layer_index=-1, domainnr = 0)
+   fill_cache_arrays(map) # fill global cache arrays if not already done
+   # first load data at edges or faces
+   values_temp=load_nc_map_slice(map,varname,itime,layer_index,domainnr)
+   location=map[1][varname].atts["location"]
+   if location=="face" #no interpolation needed, call load_nc_map_slice directly
+      return values_temp
+   end
+   #so interpolation is needed here
+   result = []
+   domainnr != 0 || (domainnr = 1:length(map))  
+   for i in domainnr
+      nfaces=nfaces_cache[i]
+      result_part=zeros(nfaces)
+      map_edges_to_faces_2d!(result_part,values_temp[i],edge_faces_cache[i],edge_type_cache[i])
+      push!(result,result_part)
+   end
+   return result
+end
 
 """
    times=get_times(map,Dates.DateTime(2019,1,1))

src/dflow_map_interp_to_zarr.jl

@@ -18,7 +18,7 @@ debuglevel=1
 #
 # these variables are added to the configuration by default
 try_vars = ["waterlevel","x_velocity","y_velocity","z_velocity","salinity","temperature","waterdepth",
-    "z_center_3d","z_iface_3d"] 
+    "eddy_visc_z", "z_center_3d","z_iface_3d"] 
 # default settings per variable
 defaults = Dict(
     "waterlevel" => Dict(
@@ -57,6 +57,12 @@ defaults = Dict(
         "add_offset" => 0.0,
         "data_type" => "Int16",
         "_FillValue" => 9999 ),
+    "eddy_visc_z" => Dict(
+        "name" => "eddy_visc_z",
+        "scale_factor" => 1.0E-6, 
+        "add_offset" => 0.0,
+        "data_type" => "Int16",
+        "_FillValue" => 9999 ),
     "waterdepth" => Dict(
         "name" => "waterdepth",
         "scale_factor" => 0.01, 
@@ -104,6 +110,7 @@ aliases=Dict{String,Vector{String}}(
     "salinity"    => ["sa1","mesh2d_sa1"], #sa1 not sal (one not L)?
     "waterdepth"  => ["waterdepth","mesh2d_waterdepth"],
     "temperature" => ["tem1","mesh2d_tem1"], 
+    "eddy_visc_z" => ["vicwwu","mesh2d_vicwwu"], 
     "x_center"    => ["FlowElem_xcc","mesh2d_face_x"],
     "y_center"    => ["FlowElem_ycc","mesh2d_face_y"],
     "z_center"    => ["mesh2d_layer_z","LayCoord_cc"], #1d
@@ -465,6 +472,7 @@ function scale_values(in_values,in_dummy,out_type,out_offset,out_scale,out_dummy
     return out_values
  end
 
+
 function interp_var(inputs::Vector{NcFile},interp::Interpolator,output::ZGroup,varname::String,xpoints,ypoints,config,dumval=NaN)
     println("interpolating variable name=$(varname)")
     globals=config["global"]
@@ -531,7 +539,7 @@ function interp_var(inputs::Vector{NcFile},interp::Interpolator,output::ZGroup,v
             var=out_temp
         end
     else #is 3D
-        if hastime
+        if hastime #x y z t
             nz=in_size[1]
             varatts["_ARRAY_DIMENSIONS"]=["time","z","y","x"]
             varatts["coordinates"]="time z_center y_center x_center"
@@ -551,13 +559,19 @@ function interp_var(inputs::Vector{NcFile},interp::Interpolator,output::ZGroup,v
                 print("|")
                 for ilayer in 1:nz
                     print(".")
-                    @time in_temp_uninterpolated=load_nc_map_slice(inputs,ncname,it,ilayer)
+                    #in_temp_uninterpolated=[]
+                    # check if variable is given on edges
+                    if varatts["location"]=="edge" #not face
+                        @time in_temp_uninterpolated=load_nc_map_slice_at_faces(inputs,ncname,it,ilayer)
+                    else
+                        @time in_temp_uninterpolated=load_nc_map_slice(inputs,ncname,it,ilayer)
+                    end
                     @time in_temp=interpolate(interp,xpoints,ypoints,in_temp_uninterpolated,dumval)
                     @time out_temp=scale_values(in_temp,in_dummy,out_type,out_offset,out_scale,out_dummy)
                     var[:,:,ilayer,it]=out_temp[:,:]
                 end
             end
-        else
+        else # x y z
             nz=in_size[1]
             varatts["_ARRAY_DIMENSIONS"]=["z","y","x"]
             varatts["coordinates"]="z_center y_center x_center"

test/test_delft3dfm_to_zarr.jl

@@ -165,6 +165,8 @@ function test_lock_exchange()
    @test haskey(config,"waterlevel")
    @test haskey(config,"x_velocity")
 
+   @test haskey(config,"z_velocity")
+
    # run interpolation with default config file
    result=Base.invokelatest(main,[configfile]) #invokelatest is needed for include in function
 
@@ -186,6 +188,14 @@ function test_lock_exchange()
    @test haskey(z,"y_velocity")
    v = z["y_velocity"][:,:,:,:]
    @test size(v) == (344,1,40,16)
+   @test haskey(z,"z_velocity")
+   w = z["z_velocity"][:,:,:,:]
+   @test size(w) == (344,1,41,16) # one more interface than layers
+
+   @test haskey(z,"eddy_visc_z")
+   nu = z["eddy_visc_z"][:,:,:,:]
+   @test size(nu) == (344,1,41,16) # one more interface than layers  
+
    @test haskey(z,"waterlevel")
    h = z["waterlevel"][:,:,:]
    @test size(h) == (344,1,16)
@@ -219,8 +229,3 @@ test_estuary()
 test_f34()
 test_lock_exchange()
 
-#test2 needs large input files that are not present in the repository.
-#only run tests if the files have been added (manually for now)
-# if isfile(joinpath("../test_data","DCSM-FM_0_5nm_0000_map.nc"))
-#    test2()
-# end

test/test_dflow.jl

@@ -110,6 +110,45 @@ function test3()
    @test u_value==0.04905390667374792
 end
 
+function test4()
+   # open file for 3d lock exchange
+   dflow_map=load_nc_info("../test_data","locxxz_fullgrid_map.nc")
+   @test length(dflow_map)==1
+
+   #load and inspect a layer of salinity
+   ilayer=5
+   itime=10
+   sal=load_nc_map_slice(dflow_map,"mesh2d_sa1",itime,ilayer)
+   @test length(sal)==1
+   @test length(sal[1])==400
+   @test sal[1][12]≈5.000000000000
+   # inspect vertical eddy viscosity
+   ilayer=2
+   itime=10
+   nu=load_nc_map_slice(dflow_map,"mesh2d_vicwwu",itime,ilayer)
+   @test length(nu)==1
+   @test length(nu[1])==1201
+   @test nu[1][100]≈2.181251740966235e-5
+   # inspect vertical eddy viscosity at FACES
+   ilayer=2
+   itime=10
+   nu_face=load_nc_map_slice_at_faces(dflow_map,"mesh2d_vicwwu",itime,ilayer)
+   @test length(nu_face)==1
+   @test length(nu_face[1])==400
+   @test nu_face[1][100]≈2.2144234788637896e-5
+   #load and inspect a layer of x_velocity
+   ilayer=1
+   itime=10
+   u=load_nc_map_slice(dflow_map,"mesh2d_ucx",itime,ilayer)
+   @test length(u)==1
+   @test length(u[1])==400
+   @test u[1][100]≈0.6020406834128101
+end
+
+#
+# run tests
+#
+
 test1()
 
 #test2 needs large input files that are not present in the repository.
@@ -119,3 +158,5 @@ if isfile(joinpath("../test_data","DCSM-FM_0_5nm_0000_map.nc"))
 end
 
 test3()
+
+test4()