Handle #49

Project.toml

@@ -9,6 +9,7 @@ Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LightXML = "9c8b4983-aa76-5018-a973-4c85ecc9e179"
+NaturalNeighbours = "f16ad982-4edb-46b1-8125-78e5a8b5a9e6"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"

src/Batsrus.jl

@@ -3,7 +3,8 @@ module Batsrus
 #
 # Hongyang Zhou, [email protected]
 
-using Printf, Reexport, Requires 
+using Printf, Reexport, Requires
+using NaturalNeighbours: interpolate, Triangle
 
 export BATLData,
    load, readlogdata, readtecdata, showhead, # io

src/plot/utility.jl

@@ -7,19 +7,21 @@
 Return 2D slices of data `var` from `bd`. If `plotrange` is not set, output data resolution
 is the same as the original. If `innermask==true`, then the inner boundary cells are set to
 NaN. If `innermask == true` but the rbody parameter is not found in the header, we use the
-keyword `rbody` as the inner radius.
+keyword `rbody` as the inner radius. If `useMatplotlib==false`, a native Julia scattered
+interpolation is used (but is somehow slower than Matplotlib).
 """
 function getdata2d(bd::BATLData, var::AbstractString,
    plotrange::Vector=[-Inf, Inf, -Inf, Inf], plotinterval::Real=Inf;
-   innermask::Bool=false, rbody=1.0)
+   innermask::Bool=false, rbody::Real=1.0, useMatplotlib::Bool=true)
    x, w, ndim = bd.x, bd.w, bd.head.ndim
    @assert ndim == 2 "data must be in 2D!"
 
    varIndex_ = findindex(bd, var)
 
    if bd.head.gencoord # Generalized coordinates
       X, Y = eachslice(x, dims=3)
-      W = @view w[:,:,varIndex_]
+      X, Y = vec(X), vec(Y)
+      W = @views w[:,:,varIndex_] |> vec
 
       adjust_plotrange!(plotrange, extrema(X), extrema(Y))
       # Set a heuristic value if not set
@@ -28,11 +30,16 @@ function getdata2d(bd::BATLData, var::AbstractString,
       end
       xi = range(plotrange[1], stop=plotrange[2], step=plotinterval)
       yi = range(plotrange[3], stop=plotrange[4], step=plotinterval)
-      # Perform linear interpolation of the data (x,y) on grid(xi,yi)
-      triang = @views matplotlib.tri.Triangulation(X[:], Y[:])
-      interpolator = @views matplotlib.tri.LinearTriInterpolator(triang, W[:])
-      Xi, Yi = meshgrid(xi, yi)
-      Wi = interpolator(Xi, Yi)
+
+      if useMatplotlib
+         triang = matplotlib.tri.Triangulation(X, Y)
+         # Perform linear interpolation on the triangle mesh
+         interpolator = matplotlib.tri.LinearTriInterpolator(triang, W)
+         Xi, Yi = meshgrid(xi, yi)
+         Wi = interpolator(Xi, Yi)
+      else
+         xi, yi, Wi = interpolate2d_generalized_coords(X, Y, W, plotrange, plotinterval)
+      end
    else # Cartesian coordinates
       xrange = range(x[1,1,1], x[end,1,1], length=size(x,1))
       yrange = range(x[1,1,2], x[1,end,2], length=size(x,2))
@@ -73,7 +80,7 @@ function getdata2d(bd::BATLData, var::AbstractString,
             end
          end
       end
-      
+
    end
 
    xi, yi, Wi
@@ -111,4 +118,16 @@ function adjust_plotrange!(plotrange, xlimit, ylimit)
    plotrange[4] = ifelse(isinf(plotrange[4]), ylimit[2], plotrange[4])
 
    return
+end
+
+"Perform Triangle interpolation of 2D data `W` on grid `X`, `Y`."
+function interpolate2d_generalized_coords(X::T, Y::T, W::T,
+   plotrange::Vector{<:AbstractFloat}, plotinterval::Real) where
+   {T<:AbstractVector}
+   xi = range(plotrange[1], stop=plotrange[2], step=plotinterval)
+   yi = range(plotrange[3], stop=plotrange[4], step=plotinterval)
+   itp = interpolate(X, Y, W)
+   Wi = [itp(x, y; method=Triangle()) for y in yi, x in xi]::Matrix{eltype(W)}
+
+   xi, yi, Wi
 end

test/runtests.jl

@@ -55,10 +55,12 @@ end
       @test bd["B"][:,end,end] == Float32[1.118034, -0.559017, 0.0]
    end
 
-   @testset "Reading 2D unstructured binary" begin
-      #file = "z=0_raw_1_t25.60000_n00000258.out"
-      #bd = load(file)
-      #TODO test getdata2d on gencoord
+   @testset "Reading 2D unstructured" begin
+      file = "bx0_mhd_6_t00000100_n00000352.out"
+      bd = load(joinpath(datapath, file))
+      plotrange = [-Inf, Inf, -Inf, Inf]
+      x, y, w = Batsrus.getdata2d(bd, "rho", plotrange, useMatplotlib=false)
+      @test w[1,2] == 5.000018304080387
    end
 
    @testset "Reading 3D structured binary" begin
@@ -149,18 +151,14 @@ end
          else
             @test c.get_array()[end] == 0.9750000000000002
          end
-         c = PyPlot.contourf(bd, "rho", innermask=true)
+         c = @suppress_err PyPlot.contourf(bd, "rho", innermask=true)
          @static if matplotlib.__version__ < "3.8"
             @test c.get_array()[end] == 1.0500000000000003
          else
             @test c.get_array()[end] == 0.9750000000000002
          end
          c = PyPlot.contour(bd, "rho")
-         @static if matplotlib.__version__ < "3.8"
-            @test c.get_array()[end] == 1.0500000000000003
-         else
-            @test c.get_array()[end] == 0.9750000000000002
-         end
+         @test c.get_array()[end] == 1.0500000000000003
          c=  PyPlot.contour(bd, "rho"; levels=[1.0])
          @test c.get_array()[end] == 1.0
          c = PyPlot.tricontourf(bd, "rho")
@@ -186,8 +184,7 @@ end
          file = "bx0_mhd_6_t00000100_n00000352.out"
          bd = load(joinpath(datapath, file))
          plotdata(bd, "P", plotmode="contbar")
-         ax = gca()
-         @test isa(ax, PyPlot.PyObject)
+         @test isa(gca(), PyPlot.PyObject)
       end
    end