Add native Julia GeoPackage implementation with improvements from rev…

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"
 CoordRefSystems = "0.18"
+DBInterface = "2.5"
 FileIO = "1.16"
 Format = "1.3"
 GRIBDatasets = "0.3, 0.4"
@@ -59,6 +62,7 @@ PlyIO = "1.1"
 PrecompileTools = "1.2"
 PrettyTables = "2.2"
 ReadVTK = "0.2"
+SQLite = "1.6"
 Shapefile = "0.13"
 StaticArrays = "1.6"
 Tables = "1.7"

src/GeoIO.jl

@@ -56,6 +56,10 @@ import ArchGDAL.GDAL
 # PROJJSON CRS
 import JSON3
 
+# SQLite for GeoPackage
+import SQLite
+import DBInterface
+
 # VTK extensions
 const VTKEXTS = [".vtu", ".vtp", ".vtr", ".vts", ".vti"]
 
@@ -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="GeoIO.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/conversion.jl

@@ -128,5 +128,8 @@ function togeometry(::GI.MultiPolygonTrait, geom, crs)
   Multi([topolygon(g, CRS) for g in GI.getgeom(geom)])
 end
 
+# Handle geometries that are already Meshes.jl types - return as-is
+geom2meshes(geom::Geometry, crs) = geom
+
 geom2meshes(geom, crs=GI.crs(geom)) = geom2meshes(GI.geomtrait(geom), geom, crs)
 geom2meshes(trait, geom, crs) = togeometry(trait, geom, crs)

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; kwargs...)
   else # fallback to GDAL
     agwrite(fname, geotable; kwargs...)
   end
@@ -63,6 +65,8 @@ function gistable(fname; layer, numtype, kwargs...)
     return GJS.read(fname; numbertype=numtype, kwargs...)
   elseif endswith(fname, ".parquet")
     return GPQ.read(fname; kwargs...)
+  elseif endswith(fname, ".gpkg")
+    return gpkgread(fname; layer, kwargs...)
   else # fallback to GDAL
     data = AG.read(fname; kwargs...)
     return AG.getlayer(data, layer - 1)

src/extra/gpkg.jl

@@ -0,0 +1,18 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+# Shared WKB type constants
+const WKB_POINT = 0x00000001
+const WKB_LINESTRING = 0x00000002
+const WKB_POLYGON = 0x00000003
+const WKB_MULTIPOINT = 0x00000004
+const WKB_MULTILINESTRING = 0x00000005
+const WKB_MULTIPOLYGON = 0x00000006
+
+const WKB_Z = 0x80000000
+const WKB_M = 0x40000000
+const WKB_ZM = 0xC0000000
+
+include("gpkg/read.jl")
+include("gpkg/write.jl")