GeoPackage writer

Project.toml

@@ -9,6 +9,7 @@ CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 CommonDataModel = "1fbeeb36-5f17-413c-809b-666fb144f157"
 CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb"
+DBInterface = "a10d1c49-ce27-4219-8d33-6db1a4562965"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8"
 GRIBDatasets = "82be9cdb-ee19-4151-bdb3-b400788d9abc"
@@ -27,6 +28,7 @@ PlyIO = "42171d58-473b-503a-8d5f-782019eb09ec"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3"
+SQLite = "0aa819cd-b072-5ff4-a722-6bc24af294d9"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
@@ -41,6 +43,7 @@ CSV = "0.10"
 Colors = "0.12, 0.13"
 CommonDataModel = "0.2, 0.3, 0.4"
 CoordRefSystems = "0.19"
+DBInterface = "2.6"
 FileIO = "1.16"
 Format = "1.3"
 GRIBDatasets = "0.3, 0.4"
@@ -59,6 +62,7 @@ PlyIO = "1.1"
 PrecompileTools = "1.2"
 PrettyTables = "3.0"
 ReadVTK = "0.2"
+SQLite = "1.6"
 Shapefile = "0.13"
 StaticArrays = "1.6"
 Tables = "1.7"

src/GeoIO.jl

@@ -38,6 +38,10 @@ import PlyIO
 # CSV format
 import CSV
 
+# Database interfaces
+import DBInterface
+import SQLite
+
 # geostats formats
 import GslibIO
 
@@ -72,7 +76,7 @@ const CDMEXTS = [".grib", ".nc"]
 const FORMATS = [
   (extension=".csv", load="CSV.jl", save="CSV.jl"),
   (extension=".geojson", load="GeoJSON.jl", save="GeoJSON.jl"),
-  (extension=".gpkg", load="ArchGDAL.jl", save="ArchGDAL.jl"),
+  (extension=".gpkg", load="ArchGDAL.jl", save="GeoIO.jl"),
   (extension=".grib", load="GRIBDatasets.jl", save=""),
   (extension=".gslib", load="GslibIO.jl", save="GslibIO.jl"),
   (extension=".jpeg", load="ImageIO.jl", save="ImageIO.jl"),
@@ -127,6 +131,7 @@ include("extra/csv.jl")
 include("extra/gdal.jl")
 include("extra/geotiff.jl")
 include("extra/gis.jl")
+include("extra/gpkg.jl")
 include("extra/img.jl")
 include("extra/msh.jl")
 include("extra/obj.jl")

src/extra/gis.jl

@@ -50,6 +50,8 @@ function giswrite(fname, geotable; warn, kwargs...)
   elseif endswith(fname, ".parquet")
     CRS = crs(domain(geotable))
     GPQ.write(fname, geotable, (:geometry,), projjson(CRS); kwargs...)
+  elseif endswith(fname, ".gpkg")
+    gpkgwrite(fname, geotable;)
   else # fallback to GDAL
     agwrite(fname, geotable; kwargs...)
   end

src/extra/gpkg.jl

@@ -0,0 +1,6 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+include("gpkg/wkb.jl")
+include("gpkg/write.jl")

src/extra/gpkg/wkb.jl

@@ -0,0 +1,112 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+function writewkbgeom(io, geom)
+  wkbtype = _wkbtype(geom)
+  write(io, htol(one(UInt8)))
+  write(io, htol(UInt32(wkbtype)))
+  geom2wkb(io, wkbtype, geom)
+end
+
+
+function _wkbtype(geometry)
+  if geometry isa Point
+    # wkbPoint 
+    return 1
+  elseif geometry isa Rope || geometry isa Ring
+    # wkbLineString
+    return 2
+  elseif geometry isa PolyArea
+    # wkbPolygon
+    return 3
+  elseif geometry isa Multi
+    # wkbMulti
+    fg = first(parent(geometry))
+    return _wkbtype(fg) + 3
+  end
+end
+
+function geom2wkb(io, wkbtype, geom)
+  if wkbtype > 3
+    multi2wkb(io, wkbtype, geom)
+  else
+    single2wkb(io, wkbtype, geom)
+  end
+end
+
+
+function single2wkb(io, wkbtype, geom)
+  # wkbPolygon
+  if wkbtype == 3
+    poly2wkb(io, [boundary(geom::PolyArea)])
+  elseif wkbtype == 2
+    coordlist = vertices(geom)
+    if typeof(geom) <: Ring
+      # wkbLineString[length+1]
+      return chainring2wkb(io, coordlist)
+    end
+    # wkbLineString
+    chainrope2wkb(io, coordlist)
+  elseif wkbtype == 1
+    # wkbPoint
+    point2wkb(io, geom)
+  else
+    throw(ErrorException("Well-Known Binary Geometry unknown: $wkbtype"))
+  end
+end
+
+function point2wkb(io, geom)
+  coordinates = CoordRefSystems.raw(coords(geom))
+  if typeof(geom) <: LatLon
+    write(io, htol(coordinates[2]))
+    write(io, htol(coordinates[1]))
+  elseif typeof(geom) <: LatLonAlt
+    write(io, htol(coordinates[2]))
+    write(io, htol(coordinates[1]))
+    write(io, htol(coordinates[3]))
+  else
+    write(io, htol(coordinates[1]))
+    write(io, htol(coordinates[2]))
+    if length(coordinates) == 3
+      write(io, htol(coordinates[3]))
+    end
+  end
+end
+
+function chainrope2wkb(io, rope)
+  write(io, htol(UInt32(length(rope))))
+  for point in rope
+    point2wkb(io, point)
+  end
+end
+
+function chainring2wkb(io, ring)
+  # add a point to close linestring
+  write(io, htol(UInt32(length(ring) + 1)))
+  for point in ring
+    point2wkb(io, point)
+  end
+  # write a point to close linestring
+  point2wkb(io, first(ring))
+end
+
+function poly2wkb(io, rings)
+  write(io, htol(UInt32(length(rings))))
+  for ring in rings
+    points = vertices(ring)
+    chainrope2wkb(io, points)
+  end
+end
+
+function multi2wkb(io, multiwkbtype, geoms)
+  # `geoms` is treated as a single [`Geometry`]
+  # `parent(geoms)` returns the collection of geometries with the same types
+  write(io, htol(UInt32(length(parent(geoms)))))
+  for geom in parent(geoms)
+    write(io, one(UInt8))
+    # wkbGeometryType + 3 = Multi-wkbGeometryType
+    write(io, multiwkbtype - 3)
+    single2wkb(io, multiwkbtype - 3, geom)
+  end
+end