write to shapefile (#15)

* write to shapefile

* remove initial conversion to MeshData

* Update src/GeoTables.jl

* add information about keyword arguments

* better description of examples

* docstring for keyargs

* Update README.md

Co-authored-by: Júlio Hoffimann <julio.hoffimann@gmail.com>

* remove temporal fix of Shapefile.jl

* tests for writing

* use datadir for writing

---------

Co-authored-by: Júlio Hoffimann <julio.hoffimann@gmail.com>

Author: ErickChacon

README.md

@@ -27,6 +27,18 @@ table = GeoTables.load("file.shp")
 GeoTables.save("file.geojson", table)
 ```
 
+Additional keyword arguments can be passed to `load` and `save` functions. Valid
+arguments are those accepted by `GeoJSON.read`, `GeoJSON.write`, `Shapefile.Table`,
+`Shapefile.write` and `ArchGDAL.read`. See below some examples:
+
+```julia
+# read `.geojson` geometries with Float64 precision
+table = GeoTables.load("file.geojson", numbertype = Float64)
+
+# force writing on existing `.shp` file
+GeoTables.save("file.shp", table, force = true)
+```
+
 ### Loading data from GADM
 
 The `gadm` function (down)loads data from the GADM dataset:

src/GeoTables.jl

@@ -19,44 +19,53 @@ include("conversion.jl")
 include("geotable.jl")
 
 """
-    load(fname, layer=0)
+    load(fname, layer=0, kwargs...)
 
 Load geospatial table from file `fname` and convert the
 `geometry` column to Meshes.jl geometries. Optionally,
-specify the layer of geometries to read within the file.
+specify the layer of geometries to read within the file and
+keyword arguments accepted by `Shapefile.Table`, `GeoJSON.read`
+and `ArchGDAL.read`. For example, use `numbertype = Float64` to
+read `.geojson` geometries with Float64 precision.
 
 ## Supported formats
 
 - `*.shp` via Shapefile.jl
 - `*.geojson` via GeoJSON.jl
 - Other formats via ArchGDAL.jl
 """
-function load(fname, layer=0)
+function load(fname; layer=0, kwargs...)
   if endswith(fname, ".shp")
-    table = SHP.Table(fname)
+    table = SHP.Table(fname; kwargs...)
   elseif endswith(fname, ".geojson")
     data  = Base.read(fname)
-    table = GJS.read(data)
+    table = GJS.read(data; kwargs...)
   else # fallback to GDAL
-    data  = AG.read(fname)
+    data  = AG.read(fname; kwargs...)
     table = AG.getlayer(data, layer)
   end
   GeoTable(table)
 end
 
 """
-    save(fname, geotable)
+    save(fname, geotable; kwargs...)
 
 Save geospatial table to file `fname` using the
 appropriate format based on the file extension.
+Optionally, specify keyword arguments accepted by
+`Shapefile.write` and `GeoJSON.write`. For example, use
+`force = true` to force writing on existing `.shp` file.
 
 ## Supported formats
 
+- `*.shp` via Shapefile.jl
 - `*.geojson` via GeoJSON.jl
 """
-function save(fname, geotable)
-  if endswith(fname, ".geojson")
-    GJS.write(fname, geotable)
+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"))
   end

src/conversion.jl

@@ -20,15 +20,19 @@ GI.ncoord(::GI.PointTrait, p::Point) = embeddim(p)
 GI.getcoord(::GI.PointTrait, p::Point) = coordinates(p)
 GI.getcoord(::GI.PointTrait, p::Point, i) = coordinates(p)[i]
 
-GI.ngeom(::Any, s::Segment) = nvertices(s)
-GI.getgeom(::Any, s::Segment, i) = vertices(s)[i]
+GI.ncoord(::GI.LineTrait, s::Segment) = embeddim(s)
+GI.ngeom(::GI.LineTrait, s::Segment) = nvertices(s)
+GI.getgeom(::GI.LineTrait, s::Segment, i) = vertices(s)[i]
 
-GI.ngeom(::Any, c::Chain) = nvertices(c)
-GI.getgeom(::Any, c::Chain, i) = vertices(c)[i]
+GI.ncoord(::GI.LineStringTrait, c::Chain) = embeddim(c)
+GI.ngeom(::GI.LineStringTrait, c::Chain) = nvertices(c)
+GI.getgeom(::GI.LineStringTrait, c::Chain, i) = vertices(c)[i]
 
-GI.ngeom(::Any, p::Polygon) = length(chains(p))
-GI.getgeom(::Any, p::Polygon, i) = chains(p)[i]
+GI.ncoord(::GI.PolygonTrait, p::Polygon) = embeddim(p)
+GI.ngeom(::GI.PolygonTrait, p::Polygon) = length(chains(p))
+GI.getgeom(::GI.PolygonTrait, p::Polygon, i) = chains(p)[i]
 
+GI.ncoord(::Any, m::Multi) = embeddim(m)
 GI.ngeom(::Any, m::Multi) = length(collect(m))
 GI.getgeom(::Any, m::Multi, i) = collect(m)[i]
 

test/runtests.jl

@@ -201,10 +201,62 @@ datadir = joinpath(@__DIR__,"data")
   end
 
   @testset "save" begin
-    table = GeoTables.load(joinpath(datadir,"path.shp"))
-    GeoTables.save(joinpath(datadir,"path.geojson"), table)
-    table = GeoTables.load(joinpath(datadir,"zone.shp"))
-    GeoTables.save(joinpath(datadir,"path.geojson"), table)
+
+    @testset "points" begin
+      table = GeoTables.load(joinpath(datadir, "points.geojson"), )
+      GeoTables.save(joinpath(datadir, "tpoints.geojson"), table)
+      GeoTables.save(joinpath(datadir, "tpoints.shp"), table, force = true)
+
+      table = GeoTables.load(joinpath(datadir, "points.gpkg"))
+      GeoTables.save(joinpath(datadir, "tpoints.geojson"), table)
+      GeoTables.save(joinpath(datadir, "tpoints.shp"), table, force = true)
+
+      table = GeoTables.load(joinpath(datadir, "points.shp"))
+      GeoTables.save(joinpath(datadir, "tpoints.geojson"), table)
+      GeoTables.save(joinpath(datadir, "tpoints.shp"), table, force = true)
+    end
+
+    @testset "lines" begin
+      table = GeoTables.load(joinpath(datadir, "lines.geojson"))
+      GeoTables.save(joinpath(datadir, "tlines.geojson"), table)
+      # GeoTables.save(joinpath(datadir, "tlines.shp"), table, force = true)
+
+      table = GeoTables.load(joinpath(datadir, "lines.gpkg"))
+      GeoTables.save(joinpath(datadir, "tlines.geojson"), table)
+      # GeoTables.save(joinpath(datadir, "tlines.shp"), table, force = true)
+
+      table = GeoTables.load(joinpath(datadir, "lines.shp"))
+      GeoTables.save(joinpath(datadir, "tlines.geojson"), table)
+      GeoTables.save(joinpath(datadir, "tlines.shp"), table, force = true)
+    end
+
+    @testset "polygons" begin
+      table = GeoTables.load(joinpath(datadir, "polygons.geojson"))
+      GeoTables.save(joinpath(datadir, "tpolygons.geojson"), table)
+      # GeoTables.save(joinpath(datadir, "tpolygons.shp"), table, force = true)
+
+      table = GeoTables.load(joinpath(datadir, "polygons.gpkg"))
+      GeoTables.save(joinpath(datadir, "tpolygons.geojson"), table)
+      # GeoTables.save(joinpath(datadir, "tpolygons.shp"), table, force = true)
+
+      table = GeoTables.load(joinpath(datadir, "polygons.shp"))
+      GeoTables.save(joinpath(datadir, "tpolygons.geojson"), table)
+      GeoTables.save(joinpath(datadir, "tpolygons.shp"), table, force = true)
+    end
+
+    @testset "multipolygons" begin
+      table = GeoTables.load(joinpath(datadir,"path.shp"))
+      GeoTables.save(joinpath(datadir, "tpath.geojson"), table)
+      GeoTables.save(joinpath(datadir, "tpath.shp"), table, force = true)
+
+      table = GeoTables.load(joinpath(datadir,"zone.shp"))
+      GeoTables.save(joinpath(datadir, "tzone.geojson"), table)
+      GeoTables.save(joinpath(datadir, "tzone.shp"), table, force = true)
+
+      table = GeoTables.load(joinpath(datadir,"ne_110m_land.shp"))
+      GeoTables.save(joinpath(datadir, "tne_110m_land.geojson"), table)
+      GeoTables.save(joinpath(datadir, "tne_110m_land.shp"), table, force = true)
+    end
   end
 
   @testset "gadm" begin