VTK Format: Add support for reading ".vtu" files

Project.toml

@@ -16,8 +16,10 @@ Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
 PlyIO = "42171d58-473b-503a-8d5f-782019eb09ec"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
+ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
+VTKBase = "4004b06d-e244-455f-a6ce-a5f9919cc534"
 
 [compat]
 ArchGDAL = "0.10"
@@ -32,6 +34,8 @@ Meshes = "0.35"
 PlyIO = "1.1"
 PrecompileTools = "1.2"
 PrettyTables = "2.2"
+ReadVTK = "0.1"
 Shapefile = "0.10"
 Tables = "1.7"
+VTKBase = "1.0"
 julia = "1.9"

src/GeoIO.jl

@@ -18,6 +18,10 @@ import PlyIO
 # geostats formats
 import GslibIO
 
+# VTK formats
+import ReadVTK
+import VTKBase.VTKCellTypes
+
 # GIS formats
 import Shapefile as SHP
 import GeoJSON as GJS
@@ -26,11 +30,15 @@ import GeoParquet as GPQ
 import GeoInterface as GI
 
 # image extensions
-const IMGEXT = (".png", ".jpg", ".jpeg", ".tif", ".tiff")
+const IMGEXTS = [".png", ".jpg", ".jpeg", ".tif", ".tiff"]
+
+# VTK extensions
+const VTKEXTS = [".vtu"]
 
 # supported formats
 const FORMATS = [
   (extension=".ply", load="PlyIO.jl", save=""),
+  (extension=".vtu", load="ReadVTK.jl", save=""),
   (extension=".kml", load="ArchGDAL.jl", save=""),
   (extension=".gslib", load="GslibIO.jl", save="GslibIO.jl"),
   (extension=".shp", load="Shapefile.jl", save="Shapefile.jl"),
@@ -68,6 +76,7 @@ include("conversion.jl")
 # extra code for backends
 include("extra/ply.jl")
 include("extra/gdal.jl")
+include("extra/vtk.jl")
 
 # user functions
 include("load.jl")

src/extra/vtk.jl

@@ -0,0 +1,92 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+const GEOMTYPE = Dict(
+  VTKCellTypes.VTK_LINE => Segment,
+  VTKCellTypes.VTK_TRIANGLE => Triangle,
+  VTKCellTypes.VTK_PIXEL => Quadrangle,
+  VTKCellTypes.VTK_QUAD => Quadrangle,
+  VTKCellTypes.VTK_POLYGON => Ngon,
+  VTKCellTypes.VTK_TETRA => Tetrahedron,
+  VTKCellTypes.VTK_VOXEL => Hexahedron,
+  VTKCellTypes.VTK_HEXAHEDRON => Hexahedron,
+  VTKCellTypes.VTK_PYRAMID => Pyramid
+)
+
+function vtkread(fname)
+  if endswith(fname, ".vtu")
+    vturead(fname)
+  else
+    error("unsupported VTK file format")
+  end
+end
+
+function vturead(fname)
+  vtk = ReadVTK.VTKFile(fname)
+
+  # get points
+  coords = ReadVTK.get_points(vtk)
+  points = [Point(Tuple(c)) for c in eachcol(coords)]
+
+  # get connectivity info
+  cells = ReadVTK.get_cells(vtk)
+  offsets = cells.offsets
+  connectivity = cells.connectivity
+  vtktypes = [VTKCellTypes.VTKCellType(id) for id in cells.types]
+
+  # list of connectivity indices
+  inds = map(eachindex(offsets), vtktypes) do i, vtktype
+    start = i == 1 ? 1 : (offsets[i - 1] + 1)
+    connec = Tuple(connectivity[start:offsets[i]])
+    _adjustconnec(connec, vtktype)
+  end
+
+  # list of mapped geometry types
+  types = [GEOMTYPE[vtktype] for vtktype in vtktypes]
+
+  # construct connectivity elements
+  connec = [connect(ind, G) for (ind, G) in zip(inds, types)]
+
+  # construct mesh
+  mesh = SimpleMesh(points, connec)
+
+  # extract point data
+  vtable = try
+    vtkdata = ReadVTK.get_point_data(vtk)
+    _astable(vtkdata)
+  catch
+    nothing
+  end
+
+  # extract element data
+  etable = try
+    vtkdata = ReadVTK.get_cell_data(vtk)
+    _astable(vtkdata)
+  catch
+    nothing
+  end
+
+  # georeference
+  GeoTable(mesh; vtable, etable)
+end
+
+function _astable(vtkdata)
+  pairs = map(keys(vtkdata)) do name
+    column = ReadVTK.get_data(vtkdata[name])
+    Symbol(name) => column
+  end
+  (; pairs...)
+end
+
+# in the case of VTK_PIXEL and VTK_VOXEL
+# we need to flip vertices in the connectivity list 
+function _adjustconnec(connec, vtktype)
+  if vtktype == VTKCellTypes.VTK_PIXEL
+    (connec[1], connec[2], connec[4], connec[3])
+  elseif vtktype == VTKCellTypes.VTK_VOXEL
+    (connec[1], connec[2], connec[4], connec[3], connec[5], connec[6], connec[8], connec[7])
+  else
+    connec
+  end
+end

src/load.jl

@@ -20,7 +20,7 @@ To see supported formats, use the [`formats`](@ref) function.
 """
 function load(fname; layer=0, fix=true, kwargs...)
   # image formats
-  if any(ext -> endswith(fname, ext), IMGEXT)
+  if any(ext -> endswith(fname, ext), IMGEXTS)
     data = FileIO.load(fname) |> rotr90
     dims = size(data)
     values = (; color=vec(data))
@@ -38,6 +38,11 @@ function load(fname; layer=0, fix=true, kwargs...)
     return plyread(fname; kwargs...)
   end
 
+  # vtk formats
+  if any(ext -> endswith(fname, ext), VTKEXTS)
+    return vtkread(fname)
+  end
+
   # GIS formats
   table = if endswith(fname, ".shp")
     SHP.Table(fname; kwargs...)

src/save.jl

@@ -15,7 +15,7 @@ To see supported formats, use the [`formats`](@ref) function.
 """
 function save(fname, geotable; kwargs...)
   # image formats
-  if any(ext -> endswith(fname, ext), IMGEXT)
+  if any(ext -> endswith(fname, ext), IMGEXTS)
     grid = domain(geotable)
     @assert grid isa Grid "grid not found"
     table = values(geotable)

test/Project.toml

@@ -6,6 +6,7 @@ GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
 GeoTables = "e502b557-6362-48c1-8219-d30d308dcdb0"
 Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3"
 ReferenceTests = "324d217c-45ce-50fc-942e-d289b448e8cf"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"

test/runtests.jl

@@ -5,6 +5,7 @@ using GeoTables
 using Test, Random
 using ReferenceTests
 using Dates
+import ReadVTK
 import GeoInterface as GI
 import Shapefile as SHP
 import ArchGDAL as AG
@@ -150,6 +151,17 @@ end
       @test isnan(table."Water Saturation"[end])
     end
 
+    @testset "VTK" begin
+      file = ReadVTK.get_example_file("celldata_appended_binary_compressed.vtu")
+      table = GeoIO.load(file)
+      @test table.geometry isa SimpleMesh
+      @test eltype(table.cell_ids) <: Int
+      @test eltype(table.element_ids) <: Int
+      @test eltype(table.levels) <: Int
+      @test eltype(table.indicator_amr) <: Float64
+      @test eltype(table.indicator_shock_capturing) <: Float64
+    end
+
     @testset "Shapefile" begin
       table = GeoIO.load(joinpath(datadir, "points.shp"))
       @test length(table.geometry) == 5