Add agwrite function

src/GeoTables.jl

@@ -17,6 +17,7 @@ import GeoInterface as GI
 
 include("conversion.jl")
 include("geotable.jl")
+include("agwrite.jl")
 
 """
     load(fname, layer=0, lazy=false, kwargs...)
@@ -66,14 +67,15 @@ Optionally, specify keyword arguments accepted by
 
 - `*.shp` via Shapefile.jl
 - `*.geojson` via GeoJSON.jl
+- Other formats via ArchGDAL.jl
 """
 function save(fname, geotable; kwargs...)
   if endswith(fname, ".shp")
     SHP.write(fname, geotable; kwargs...)
   elseif endswith(fname, ".geojson")
     GJS.write(fname, geotable; kwargs...)
-  else
-    throw(ErrorException("file format not supported"))
+  else # fallback to GDAL
+    agwrite(fname, geotable; kwargs...)
   end
 end
 

src/agwrite.jl

@@ -0,0 +1,57 @@
+# adapted from https://github.com/evetion/GeoDataFrames.jl/blob/master/src/io.jl
+# and from https://github.com/yeesian/ArchGDAL.jl/blob/master/test/test_tables.jl#L264
+
+const DRIVER = AG.extensions()
+
+asstrings(options::Dict{<:AbstractString,<:AbstractString}) =
+  [uppercase(String(k)) * "=" * String(v) for (k, v) in options]
+
+function agwrite(fname, geotable; layername="data", options=Dict("geometry_name" => "geometry"))
+  geoms = domain(geotable)
+  table = values(geotable)
+  rows = Tables.rows(table)
+  schema = Tables.schema(table)
+
+  # Set geometry name in options
+  if !haskey(options, "geometry_name")
+    options["geometry_name"] = "geometry"
+  end
+
+  ext = last(splitext(fname))
+  driver = AG.getdriver(DRIVER[ext])
+  optionlist = asstrings(options)
+  agtypes = map(schema.types) do type
+    try
+      T = nonmissingtype(type)
+      convert(AG.OGRFieldType, T)
+    catch
+      error("type $type not supported")
+    end
+  end
+
+  AG.create(fname; driver) do dataset
+    AG.createlayer(; dataset, name=layername, options=optionlist) do layer
+      for (name, type) in zip(schema.names, agtypes)
+        AG.addfielddefn!(layer, String(name), type)
+      end
+
+      for (row, geom) in zip(rows, geoms)
+        AG.addfeature(layer) do feature
+          for name in schema.names
+            x = Tables.getcolumn(row, name)
+            i = AG.findfieldindex(feature, name)
+            if ismissing(x)
+              AG.setfieldnull!(feature, i)
+            else
+              AG.setfield!(feature, i, x)
+            end
+          end
+
+          AG.setgeom!(feature, GI.convert(AG.IGeometry, geom))
+        end
+      end
+    end
+  end
+
+  fname
+end

test/runtests.jl

@@ -309,6 +309,23 @@ savedir = mktempdir()
         newtable = GeoTables.load(joinpath(savedir, "t$file.shp"))
       end
     end
+
+    @testset "agwrite" begin
+      table = GeoTables.load(joinpath(datadir, "lines.geojson"), numbertype=Float64)
+      GeoTables.save(joinpath(savedir, "ag-lines.gpkg"), table)
+      agtable = GeoTables.load(joinpath(savedir, "ag-lines.gpkg"))
+      @test agtable == table
+
+      table = GeoTables.load(joinpath(datadir, "points.geojson"), numbertype=Float64)
+      GeoTables.save(joinpath(savedir, "ag-points.gpkg"), table)
+      agtable = GeoTables.load(joinpath(savedir, "ag-points.gpkg"))
+      @test agtable == table
+
+      table = GeoTables.load(joinpath(datadir, "polygons.geojson"), numbertype=Float64)
+      GeoTables.save(joinpath(savedir, "ag-polygons.gpkg"), table)
+      agtable = GeoTables.load(joinpath(savedir, "ag-polygons.gpkg"))
+      @test agtable == table
+    end
   end
 
   @testset "gadm" begin