From ec470e77ad3552d8f31b0eca415009f74b1d1c5a Mon Sep 17 00:00:00 2001 From: Hongyang Zhou Date: Mon, 17 Jun 2024 16:58:15 -0400 Subject: [PATCH] Handle #49 --- Project.toml | 1 + src/Batsrus.jl | 3 ++- src/plot/utility.jl | 37 ++++++++++++++++++++++++++++--------- test/runtests.jl | 21 +++++++++------------ 4 files changed, 40 insertions(+), 22 deletions(-) diff --git a/Project.toml b/Project.toml index 481f3f9f..3f766b9f 100644 --- a/Project.toml +++ b/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" diff --git a/src/Batsrus.jl b/src/Batsrus.jl index 6e795483..1f31e3cb 100644 --- a/src/Batsrus.jl +++ b/src/Batsrus.jl @@ -3,7 +3,8 @@ module Batsrus # # Hongyang Zhou, hyzhou@umich.edu -using Printf, Reexport, Requires +using Printf, Reexport, Requires +using NaturalNeighbours: interpolate, Triangle export BATLData, load, readlogdata, readtecdata, showhead, # io diff --git a/src/plot/utility.jl b/src/plot/utility.jl index 7b2cda3d..225b5cfc 100644 --- a/src/plot/utility.jl +++ b/src/plot/utility.jl @@ -7,11 +7,12 @@ 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!" @@ -19,7 +20,8 @@ function getdata2d(bd::BATLData, var::AbstractString, 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 \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 7ee10687..39f6f83b 100644 --- a/test/runtests.jl +++ b/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