Add support for saving/loading geotables without attributes in format…

…s that support data (#77)

* Add support for saving/loading geotables without attributes in formats that support data

* Update dependencies

Project.toml

@@ -42,8 +42,8 @@ GRIBDatasets = "0.3"
 GeoInterface = "1.0"
 GeoJSON = "0.7, 0.8"
 GeoParquet = "0.1, 0.2"
-GeoTables = "1.7"
-GslibIO = "1.4"
+GeoTables = "1.18.5"
+GslibIO = "1.4.10"
 ImageIO = "0.6"
 Meshes = "0.41"
 NCDatasets = "0.13, 0.14"
@@ -52,7 +52,7 @@ PrecompileTools = "1.2"
 PrettyTables = "2.2"
 ReadVTK = "0.2"
 Rotations = "1.6"
-Shapefile = "0.10, 0.11, 0.12"
+Shapefile = "0.12.1"
 StaticArrays = "1.6"
 Tables = "1.7"
 TransformsBase = "1.4"

src/extra/csv.jl

@@ -22,16 +22,14 @@ end
 function csvwrite(fname, geotable; coords=nothing, floatformat=nothing, kwargs...)
   dom = domain(geotable)
   tab = values(geotable)
-  cols = Tables.columns(tab)
-  names = Tables.columnnames(tab)
   Dim = embeddim(dom)
 
   if Dim > 3
     throw(ArgumentError("embedding dimensions greater than 3 are not supported"))
   end
 
   cnames = if isnothing(coords)
-    _cnames(Dim, names)
+    _cnames(Dim)
   else
     if length(coords) ≠ Dim
       throw(ArgumentError("the number of coordinate names must be equal to $Dim (embedding dimension)"))
@@ -40,12 +38,20 @@ function csvwrite(fname, geotable; coords=nothing, floatformat=nothing, kwargs..
   end
 
   points = [centroid(dom, i) for i in 1:nelements(dom)]
-  cpairs = map(cnames, 1:Dim) do name, d
-    name => [coordinates(p)[d] for p in points]
+  ccolumns = map(1:Dim) do d
+    [coordinates(p)[d] for p in points]
+  end
+
+  newtab = if isnothing(tab)
+    (; zip(cnames, ccolumns)...)
+  else
+    cols = Tables.columns(tab)
+    names = Tables.columnnames(tab)
+    ucnames = uniquenames(names, cnames)
+    pairs = (nm => Tables.getcolumn(cols, nm) for nm in names)
+    (; zip(ucnames, ccolumns)..., pairs...)
   end
 
-  pairs = (nm => Tables.getcolumn(cols, nm) for nm in names)
-  newtab = (; cpairs..., pairs...)
 
   transform(col, val) = _floatformat(val, floatformat)
   CSV.write(fname, newtab; transform, kwargs...)
@@ -54,20 +60,12 @@ end
 _floatformat(val, format) = val
 _floatformat(val::AbstractFloat, format) = isnothing(format) ? val : generate_formatter(format)(val)
 
-function _cnames(Dim, names)
-  cnames = if Dim == 1
+function _cnames(Dim)
+  if Dim == 1
     [:x]
   elseif Dim == 2
     [:x, :y]
   else
     [:x, :y, :z]
   end
-
-  # make unique
-  map(cnames) do name
-    while name ∈ names
-      name = Symbol(name, :_)
-    end
-    name
-  end
 end

src/extra/geotiff.jl

@@ -28,6 +28,10 @@ function geotiffwrite(fname, geotable; kwargs...)
   dims = size(grid)
 
   table = values(geotable)
+  if isnothing(table)
+    throw(ArgumentError("GeoTiff format needs data to save"))
+  end
+
   cols = Tables.columns(table)
   names = Tables.columnnames(cols)
 

src/extra/vtkwrite.jl

@@ -80,27 +80,29 @@ _extractvals(gtb) = _extractvals(domain(gtb), values(gtb), values(gtb, 0))
 _extractvals(dom::Domain, etable, vtable) = dom, etable, vtable
 function _extractvals(subdom::SubDomain, etable, vtable)
   dom = parent(subdom)
-  inds = parentindices(subdom)
-  nelems = nelements(dom)
-
-  cols = Tables.columns(etable)
-  names = Tables.columnnames(cols)
-  pairs = map(names) do name
-    x = Tables.getcolumn(cols, name)
-    y = fill(NaN, nelems)
-    y[inds] .= x
-    name => y
-  end
 
-  mask = :MASK
-  # make unique
-  while mask ∈ names
-    mask = Symbol(mask, :_)
+  newtable = if isnothing(etable)
+    nothing
+  else
+    inds = parentindices(subdom)
+    nelems = nelements(dom)
+
+    cols = Tables.columns(etable)
+    names = Tables.columnnames(cols)
+    pairs = map(names) do name
+      x = Tables.getcolumn(cols, name)
+      y = fill(NaN, nelems)
+      y[inds] .= x
+      name => y
+    end
+
+    mask = uniquename(names, :MASK)
+    maskcol = zeros(UInt8, nelems)
+    maskcol[inds] .= 1
+
+    (; pairs..., mask => maskcol)
   end
-  maskcol = zeros(UInt8, nelems)
-  maskcol[inds] .= 1
 
-  newtable = (; pairs..., mask => maskcol)
   dom, newtable, nothing
 end
 

src/save.jl

@@ -17,11 +17,18 @@ function save(fname, geotable; kwargs...)
   # IMG formats
   if any(ext -> endswith(fname, ext), IMGEXTS)
     grid = domain(geotable)
-    @assert grid isa Grid "grid not found"
+    if !(grid isa Grid)
+      throw(ArgumentError("image formats only support grids"))
+    end
     table = values(geotable)
+    if isnothing(table)
+      throw(ArgumentError("image formats need data to save"))
+    end
     cols = Tables.columns(table)
     names = Tables.columnnames(cols)
-    @assert :color ∈ names "colors not found"
+    if :color ∉ names
+      throw(ArgumentError("color column not found"))
+    end
     colors = Tables.getcolumn(cols, :color)
     img = reshape(colors, size(grid))
     return FileIO.save(fname, img)

src/utils.jl

@@ -26,3 +26,19 @@ function geomcolumn(names)
     error("geometry column not found")
   end
 end
+
+# add "_" to `name` until it is unique compared to the table `names`
+function uniquename(names, name)
+  uname = name
+  while uname ∈ names
+    uname = Symbol(uname, :_)
+  end
+  uname
+end
+
+# make `newnames` unique compared to the table `names`
+function uniquenames(names, newnames)
+  map(newnames) do name
+    uniquename(names, name)
+  end
+end

test/runtests.jl

@@ -511,6 +511,17 @@ end
       gtb2 = GeoIO.load(img2)
       @test gtb1.geometry == gtb2.geometry
       @test psnr_equality()(gtb1.color, gtb2.color)
+
+      # error: image formats only support grids
+      file = joinpath(savedir, "error.jpg")
+      gtb = georef((; a=rand(10)), rand(Point2, 10))
+      @test_throws ArgumentError GeoIO.save(file, gtb)
+      # error: image formats need data to save
+      gtb = georef(nothing, CartesianGrid(2, 2))
+      @test_throws ArgumentError GeoIO.save(file, gtb)
+      # error: color column not found
+      gtb = georef((; a=rand(4)), CartesianGrid(2, 2))
+      @test_throws ArgumentError GeoIO.save(file, gtb)
     end
 
     @testset "STL" begin
@@ -611,11 +622,19 @@ end
     @testset "PLY" begin
       file1 = joinpath(datadir, "beethoven.ply")
       file2 = joinpath(savedir, "beethoven.ply")
-      table1 = GeoIO.load(file1)
-      GeoIO.save(file2, table1)
-      table2 = GeoIO.load(file2)
-      @test table1 == table2
-      @test values(table1, 0) == values(table2, 0)
+      gtb1 = GeoIO.load(file1)
+      GeoIO.save(file2, gtb1)
+      gtb2 = GeoIO.load(file2)
+      @test gtb1 == gtb2
+      @test values(gtb1, 0) == values(gtb2, 0)
+
+      mesh = gtb1.geometry
+      gtb1 = georef((; a=rand(nelements(mesh))), mesh)
+      file = joinpath(savedir, "plywithdata.ply")
+      GeoIO.save(file, gtb1)
+      gtb2 = GeoIO.load(file)
+      @test gtb1 == gtb2
+      @test values(gtb1, 0) == values(gtb2, 0)
     end
 
     @testset "CSV" begin
@@ -857,6 +876,9 @@ end
       file = joinpath(savedir, "error.tif")
       gtb = georef((; a=rand(8)), CartesianGrid(2, 2, 2))
       @test_throws ArgumentError GeoIO.save(file, gtb)
+      # error: GeoTiff format needs data to save
+      gtb = georef(nothing, CartesianGrid(2, 2))
+      @test_throws ArgumentError GeoIO.save(file, gtb)
       # error: all variables must have the same type
       gtb = georef((a=rand(1:9, 25), b=rand(25)), CartesianGrid(5, 5))
       @test_throws ArgumentError GeoIO.save(file, gtb)
@@ -916,6 +938,13 @@ end
     pset = PointSet(rand(Point2, 10))
     gtb1 = georef(nothing, pset)
 
+    # Shapefile
+    file = joinpath(savedir, "noattribs.shp")
+    GeoIO.save(file, gtb1)
+    gtb2 = GeoIO.load(file)
+    @test all(ismissing, gtb2[:, 1])
+    @test gtb2.geometry == gtb1.geometry
+
     # GeoJSON
     file = joinpath(savedir, "noattribs.geojson")
     GeoIO.save(file, gtb1)
@@ -936,5 +965,37 @@ end
     gtb2 = GeoIO.load(file)
     @test isnothing(values(gtb2))
     @test gtb2 == gtb1
+
+    # CSV
+    file = joinpath(savedir, "noattribs.csv")
+    GeoIO.save(file, gtb1)
+    gtb2 = GeoIO.load(file, coords=[:x, :y])
+    @test isnothing(values(gtb2))
+    @test gtb2 == gtb1
+
+    # GSLIB
+    file = joinpath(savedir, "noattribs.gslib")
+    GeoIO.save(file, gtb1)
+    gtb2 = GeoIO.load(file)
+    @test isnothing(values(gtb2))
+    @test gtb2 == gtb1
+
+    # VTK
+    file = joinpath(savedir, "noattribs.vts")
+    gtb1 = georef(nothing, CartesianGrid(10, 10))
+    GeoIO.save(file, gtb1)
+    gtb2 = GeoIO.load(file)
+    @test isnothing(values(gtb2))
+    @test gtb2 == gtb1
+
+    # MSH
+    gtb = GeoIO.load(joinpath(datadir, "tetrahedron1.msh"))
+    mesh = gtb.geometry
+    file = joinpath(savedir, "noattribs.msh")
+    gtb1 = georef(nothing, mesh)
+    GeoIO.save(file, gtb1)
+    gtb2 = GeoIO.load(file)
+    @test isnothing(values(gtb2))
+    @test gtb2 == gtb1
   end
 end