Add support for saving and loading geotables in CSV format (#38)

* Add support for saving and loading geotables in CSV format

* Apply suggestions

Project.toml

@@ -5,6 +5,7 @@ version = "1.5.6"
 
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 CommonDataModel = "1fbeeb36-5f17-413c-809b-666fb144f157"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 GRIBDatasets = "82be9cdb-ee19-4151-bdb3-b400788d9abc"
@@ -28,6 +29,7 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [compat]
 ArchGDAL = "0.10"
+CSV = "0.10"
 CommonDataModel = "0.2, 0.3"
 FileIO = "1.16"
 GRIBDatasets = "0.3"

src/GeoIO.jl

@@ -29,6 +29,9 @@ import NCDatasets
 # mesh formats
 import PlyIO
 
+# CSV format
+import CSV
+
 # geostats formats
 import GslibIO
 
@@ -51,6 +54,7 @@ const CDMEXTS = [".grib", ".nc"]
 # supported formats
 const FORMATS = [
   (extension=".ply", load="PlyIO.jl", save="PlyIO.jl"),
+  (extension=".csv", load="CSV.jl", save="CSV.jl"),
   (extension=".vtu", load="ReadVTK.jl", save="WriteVTK.jl"),
   (extension=".vtp", load="ReadVTK.jl", save="WriteVTK.jl"),
   (extension=".vtr", load="ReadVTK.jl", save="WriteVTK.jl"),
@@ -94,6 +98,7 @@ include("conversion.jl")
 
 # extra code for backends
 include("extra/ply.jl")
+include("extra/csv.jl")
 include("extra/cdm.jl")
 include("extra/gdal.jl")
 include("extra/vtkread.jl")

src/extra/csv.jl

@@ -0,0 +1,43 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+csvread(fname; coords, kwargs...) = georef(CSV.File(fname; kwargs...), coords)
+
+function csvwrite(fname, geotable; coords=nothing, kwargs...)
+  dom = domain(geotable)
+  tab = values(geotable)
+  cols = Tables.columns(tab)
+  names = Tables.columnnames(tab)
+  Dim = embeddim(dom)
+
+  cnames = if isnothing(coords)
+    _cnames(Dim, names)
+  else
+    if length(coords) ≠ Dim
+      throw(ArgumentError("the number of coordinates names must be equal to $Dim (embedding dimension)"))
+    end
+    Symbol.(coords)
+  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]
+  end
+
+  pairs = (nm => Tables.getcolumn(cols, nm) for nm in names)
+  newtab = (; cpairs..., pairs...)
+
+  CSV.write(fname, newtab; kwargs...)
+end
+
+function _cnames(Dim, names)
+  map(1:Dim) do d
+    name = Symbol(:X, d)
+    # make unique
+    while name ∈ names
+      name = Symbol(name, :_)
+    end
+    name
+  end
+end

src/load.jl

@@ -43,6 +43,11 @@ function load(fname; layer=0, fix=true, kwargs...)
     return plyread(fname; kwargs...)
   end
 
+  # CSV format
+  if endswith(fname, ".csv")
+    return csvread(fname; kwargs...)
+  end
+
   # GSLIB format
   if endswith(fname, ".gslib")
     return GslibIO.load(fname; kwargs...)

src/save.jl

@@ -37,6 +37,11 @@ function save(fname, geotable; kwargs...)
     return plywrite(fname, geotable; kwargs...)
   end
 
+  # CSV format
+  if endswith(fname, ".csv")
+    return csvwrite(fname, geotable; kwargs...)
+  end
+
   # GSLIB format
   if endswith(fname, ".gslib")
     return GslibIO.save(fname, geotable; kwargs...)

test/data/points.csv

@@ -0,0 +1,6 @@
+X1,X2,variable,code,name
+0.0,0.0,0.07336635446929285,1,word1
+1.0,1.1,0.34924148955718615,2,word2
+2.0,2.2,0.6988266836914685,3,word3
+3.0,3.3,0.6282647403425017,4,word4
+4.0,4.4,0.9149290036628314,5,word5

test/runtests.jl

@@ -142,6 +142,15 @@ end
       @test isnothing(values(table, 2))
     end
 
+    @testset "CSV" begin
+      table = GeoIO.load(joinpath(datadir, "points.csv"), coords=["X1", "X2"])
+      @test eltype(table.code) <: Integer
+      @test eltype(table.name) <: AbstractString
+      @test eltype(table.variable) <: Real
+      @test table.geometry isa PointSet
+      @test length(table.geometry) == 5
+    end
+
     @testset "GSLIB" begin
       table = GeoIO.load(joinpath(datadir, "grid.gslib"))
       @test table.geometry isa CartesianGrid
@@ -429,6 +438,39 @@ end
       @test values(table1, 0) == values(table2, 0)
     end
 
+    @testset "CSV" begin
+      file1 = joinpath(datadir, "points.csv")
+      file2 = joinpath(savedir, "points.csv")
+      gtb1 = GeoIO.load(file1, coords=[:X1, :X2])
+      GeoIO.save(file2, gtb1, coords=["X1", "X2"])
+      gtb2 = GeoIO.load(file2, coords=[:X1, :X2])
+      @test gtb1 == gtb2
+      @test values(gtb1, 0) == values(gtb2, 0)
+
+      grid = CartesianGrid(2, 2, 2)
+      gtb1 = georef((; a=rand(8)), grid)
+      file = joinpath(savedir, "grid.csv")
+      GeoIO.save(file, gtb1)
+      gtb2 = GeoIO.load(file, coords=[:X1, :X2, :X3])
+      @test gtb1.a == gtb2.a
+      @test nelements(gtb1.geometry) == nelements(gtb2.geometry)
+      @test collect(gtb2.geometry) == centroid.(gtb1.geometry)
+
+      # make coordinate names unique
+      pset = PointSet(rand(Point2, 10))
+      gtb1 = georef((X1=rand(10), X2=rand(10)), pset)
+      file = joinpath(savedir, "pset.csv")
+      GeoIO.save(file, gtb1)
+      gtb2 = GeoIO.load(file, coords=[:X1_, :X2_])
+      @test propertynames(gtb1) == propertynames(gtb2)
+      @test gtb1.X1 == gtb2.X1
+      @test gtb1.X2 == gtb2.X2
+      @test gtb1.geometry == gtb2.geometry
+
+      # throw: invalid number of coordinate names
+      @test_throws ArgumentError GeoIO.save(file, gtb1, coords=["X", "Y", "Z"])
+    end
+
     @testset "VTK" begin
       file1 = ReadVTK.get_example_file("celldata_appended_binary_compressed.vtu", output_directory=savedir)
       file2 = joinpath(savedir, "unstructured.vtu")