Add point & line interpolations; Refactor (#57)

* Bake dimensionality into BATLData

* Bump minor version

Project.toml

@@ -1,7 +1,7 @@
 name = "Batsrus"
 uuid = "e74ebddf-6ac1-4047-a0e5-c32c99e57753"
 authors = ["Hongyang Zhou <[email protected]>"]
-version = "0.5.11"
+version = "0.6.0"
 
 [deps]
 FortranFiles = "c58ffaec-ab22-586d-bfc5-781a99fd0b10"

docs/src/man/manual.md

@@ -27,6 +27,12 @@ head, data = readlogdata(logfilename)
 
 ### Data Extraction
 
+- Checking variable range
+
+```julia
+get_var_range(bd, "rho")
+```
+
 - Raw variables
 
 ```julia
@@ -41,6 +47,21 @@ v = getvars(bd, ["Bx", "By", "Bz"])
 Bmag = bd["Bmag"]
 ```
 
+- Extracting data at a given location
+
+```julia
+loc = Float32[0.0, 0.0] # The type determines the output type
+d = interp1d(bd, "rho", loc)
+```
+
+- Extracting data along a given line
+
+```julia
+point1 = Float32[-10.0, -1.0]
+point2 = Float32[10.0, 1.0]
+w = slice1d(bd, "rho", point1, point2)
+```
+
 Here is a full list of predefined derived quantities:
 
 | Derived variable name | Meaning                          | Required variable |

ext/BatsrusMakieExt/BatsrusMakieExt.jl

@@ -1,7 +1,7 @@
 module BatsrusMakieExt
 
 using Batsrus, Printf
-import Batsrus: findindex, hasunit, getunit, getdata2d
+import Batsrus: findindex, hasunit, getunit, slice2d
 import Batsrus.UnitfulBatsrus
 import Makie
 

ext/BatsrusMakieExt/typerecipe.jl

@@ -19,7 +19,7 @@ end
 "Conversion for 2D plots."
 function Makie.convert_arguments(P::Makie.GridBased, bd::BATLData, var::String;
    plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1)
-   x, y, w = getdata2d(bd, var, plotrange, plotinterval)
+   x, y, w = slice2d(bd, var, plotrange, plotinterval)
 
    unitx = getunit(bd, bd.head.variables[1])
    unity = getunit(bd, bd.head.variables[2])

src/Batsrus.jl

@@ -5,13 +5,14 @@ module Batsrus
 
 using Printf, Reexport, Requires
 using Parsers
-using NaturalNeighbours: interpolate, Triangle
+using Interpolations: cubic_spline_interpolation, BSpline, Linear, scale, interpolate
+import NaturalNeighbours as NN
 
 export BATLData,
    load, readlogdata, readtecdata, showhead, # io
    getvars, getvar, cutdata, subvolume, subsurface, # select
    Batl, convertTECtoVTU, convertIDLtoVTK, readhead, readtree, getConnectivity, # vtk
-   getdata2d # plot/utility
+   interp1d, slice1d, slice2d, get_var_range # plot/utility
 
 "Type for the file information."
 struct FileList
@@ -30,7 +31,7 @@ struct FileList
 end
 
 "Primary Batsrus data storage type."
-struct BATLData{T<:AbstractFloat}
+struct BATLData{dim, T<:AbstractFloat}
    "header information"
    head::NamedTuple
    "grid"

src/io.jl

@@ -21,31 +21,26 @@ function load(file::AbstractString; npict::Int=1, verbose::Bool=false)
    seekstart(fileID) # Rewind to start
 
    ## Read data from files
-   # Skip npict-1 snapshots (because we only want the npict-th snapshot)
+   # Skip npict-1 snapshots (since we only want the npict-th snapshot)
    skip(fileID, pictsize*(npict-1))
 
    filehead = getfilehead(fileID, filelist)
-
    # Read data
    if filelist.type == :ascii
       x, w = allocateBuffer(filehead, Float64) # why Float64?
-      getascii!(x, w, fileID, filehead)
+      getascii!(x, w, fileID)
    else
       skip(fileID, TAG) # skip record start tag
       T = filelist.type == :real4 ? Float32 : Float64
       x, w = allocateBuffer(filehead, T)
-      getbinary!(x, w, fileID, filehead)
+      getbinary!(x, w, fileID)
    end
 
    close(fileID)
 
    #setunits(filehead,"")
 
-	data = BATLData(filehead, x, w, filelist)
-
-   verbose && @info "Finished reading $(filelist.name)"
-
-   data
+	data = BATLData{Int(filehead.ndim), eltype(w)}(filehead, x, w, filelist)
 end
 
 "Read information from log file."
@@ -172,22 +167,22 @@ function readtecdata(file::AbstractString; verbose::Bool=false)
    data = Array{Float32,2}(undef, length(VARS), nNode)
 
    if nDim == 3
-	   connectivity = Array{Int32,2}(undef,8,nCell)
+	   connectivity = Array{Int32,2}(undef, 8, nCell)
    elseif nDim == 2
-	   connectivity = Array{Int32,2}(undef,4,nCell)
+	   connectivity = Array{Int32,2}(undef, 4, nCell)
    end
 
-   IsBinary = false
+   isBinary = false
    try
       Parsers.parse.(Float32, split(readline(f)))
    catch
-      IsBinary = true
+      isBinary = true
       verbose && @info "reading binary file"
    end
 
    seek(f, pt0)
 
-   if IsBinary
+   if isBinary
 	   @inbounds for i in 1:nNode
 		   read!(f, @view data[:,i])
 	   end
@@ -271,8 +266,7 @@ function getfilehead(fileID::IOStream, filelist::FileList)
 
    if type == :ascii
       headline = readline(fileID)
-      line = readline(fileID)
-      line = split(line)
+      line = readline(fileID) |> split
       it = Parsers.parse(Int, line[1])
       t = Parsers.parse(Float64, line[2])
       ndim = Parsers.parse(Int8, line[3])
@@ -326,8 +320,8 @@ end
 
 function skipline(s::IO)
    while !eof(s)
-       c = read(s, Char)
-       c == '\n' && break
+      c = read(s, Char)
+      c == '\n' && break
    end
 
    return
@@ -367,11 +361,9 @@ function getfilesize(fileID::IOStream, type::Symbol, lenstr::Int32)
       read(fileID, lenstr)
       skip(fileID, TAG)
    end
-
    # Header length
    pointer1 = position(fileID)
    headlen = pointer1 - pointer0
-
    # Calculate the snapshot size = header + data + recordmarks
    nxs = prod(nx)
    pictsize =
@@ -408,29 +400,59 @@ function allocateBuffer(filehead::NamedTuple, T::DataType)
 end
 
 "Read ascii format coordinates and data values."
-function getascii!(x, w, fileID::IOStream, filehead::NamedTuple)
-   ndim = filehead.ndim
-   Ids = CartesianIndices(size(x)[1:ndim])
-   @inbounds @views for ids in Ids
+function getascii!(x::Array{T, 2}, w, fileID::IOStream) where T
+   @inbounds @views for id in axes(x, 1)
+      temp = Parsers.parse.(Float64, split(readline(fileID)))
+      x[id] = temp[1]
+      w[id,:] = temp[2:end]
+   end
+end
+
+function getascii!(x::Array{T, 3}, w, fileID::IOStream) where T
+   @inbounds @views for id in CartesianIndices(size(x)[1:2])
       temp = Parsers.parse.(Float64, split(readline(fileID)))
-      x[ids,:] = temp[1:ndim]
-      w[ids,:] = temp[ndim+1:end]
+      x[id,:] = temp[1:2]
+      w[id,:] = temp[3:end]
    end
+end
 
-   return
+function getascii!(x::Array{T, 4}, w, fileID::IOStream) where T
+   @inbounds @views for id in CartesianIndices(size(x)[1:3])
+      temp = Parsers.parse.(Float64, split(readline(fileID)))
+      x[id,:] = temp[1:3]
+      w[id,:] = temp[4:end]
+   end
 end
 
 "Read binary format coordinates and data values."
-function getbinary!(x, w, fileID::IOStream, filehead::NamedTuple)
-   dimlast = filehead.ndim + 1
+function getbinary!(x::Array{T, 2}, w, fileID::IOStream) where T
    read!(fileID, x)
    skip(fileID, 2*TAG)
+   dimlast = 2
    @inbounds for iw in axes(w, dimlast)
       read!(fileID, selectdim(w, dimlast, iw))
       skip(fileID, 2*TAG)
    end
+end
 
-   return
+function getbinary!(x::Array{T, 3}, w, fileID::IOStream) where T
+   read!(fileID, x)
+   skip(fileID, 2*TAG)
+   dimlast = 3
+   @inbounds for iw in axes(w, dimlast)
+      read!(fileID, selectdim(w, dimlast, iw))
+      skip(fileID, 2*TAG)
+   end
+end
+
+function getbinary!(x::Array{T, 4}, w, fileID::IOStream) where T
+   read!(fileID, x)
+   skip(fileID, 2*TAG)
+   dimlast = 4
+   @inbounds for iw in axes(w, dimlast)
+      read!(fileID, selectdim(w, dimlast, iw))
+      skip(fileID, 2*TAG)
+   end
 end
 
 """

src/plot/plots.jl

@@ -3,30 +3,31 @@
 using RecipesBase
 
 # Build a recipe which acts on a custom type.
-@recipe function f(bd::BATLData, var::AbstractString;
-   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1)
-   ndim = bd.head.ndim
+@recipe function f(bd::BATLData{1, T}, var::AbstractString) where {T}
    hasunits = hasunit(bd)
 
-   if ndim == 1
-      VarIndex_ = findindex(bd, var)
-      if hasunits
-         unitx = getunit(bd, bd.head.variables[1])
-         unitw = getunit(bd, var)
-         x = bd.x .* unitx
-         y = bd.w[:,VarIndex_] .* unitw
-      else
-         x = bd.x
-         y = @view bd.w[:,VarIndex_]
-      end
+   if hasunits
+      unitx = getunit(bd, bd.head.variables[1])
+      unitw = getunit(bd, var)
+      x = bd.x .* unitx
+      y = getview(bd, var) .* unitw
+   else
+      x = bd.x
+      y = getview(bd, var)
+   end
 
-      @series begin
-         seriestype --> :path
-         x, y
-      end
-   elseif ndim == 2
-      x, y, w = getdata2d(bd, var, plotrange, plotinterval)
+   @series begin
+      seriestype --> :path
+      x, y
+   end
+end
+
+@recipe function f(bd::BATLData{2, T}, var::AbstractString;
+   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1) where {T}
+   hasunits = hasunit(bd)
 
+   x, y, w = slice2d(bd, var, plotrange, plotinterval)
+   if hasunits
       unitx = getunit(bd, bd.head.variables[1])
       unity = getunit(bd, bd.head.variables[2])
       unitw = getunit(bd, var)
@@ -36,10 +37,10 @@ using RecipesBase
          y *= unity
          w *= unitw
       end
+   end
 
-      @series begin
-         seriestype --> :contourf  # use := if you want to force it
-         x, y, w'
-      end
+   @series begin
+      seriestype --> :contourf  # use := if you want to force it
+      x, y, w'
    end
 end