Geopackage reader

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.6"
 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

@@ -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="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/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,101 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+# According to https://www.geopackage.org/spec/#r2
+# a GeoPackage should contain "GPKG" in ASCII in 
+# "application_id" field of SQLite db header
+const GP10_APPLICATION_ID = Int(0x47503130)
+const GP11_APPLICATION_ID = Int(0x47503131)
+const GPKG_APPLICATION_ID = Int(0x47504B47)
+const GPKG_1_2_VERSION = 10200
+const GPKG_1_3_VERSION = 10300
+const GPKG_1_4_VERSION = 10400
+
+# types used within the GeoPackageBinary SQL BLOBs
+@enum wkbGeometryType begin
+  wkbUnknown = 0
+  wkbPoint = 1
+  wkbLineString = 2
+  wkbPolygon = 3
+  wkbMultiPoint = 4
+  wkbMultiLineString = 5
+  wkbMultiPolygon = 6
+  wkbGeometryCollection = 7
+  wkbCircularString = 8
+  wkbCompoundCurve = 9
+  wkbCurvePolygon = 10
+  wkbMultiCurve = 11
+  wkbMultiSurface = 12
+  wkbCurve = 13
+  wkbSurface = 14
+
+  wkbNone = 100 # pure attribute records
+  wkbLinearRing = 101
+
+  # ISO SQL/MM Part 3: Spatial
+  # Z-aware types
+  wkbPointZ = 1001
+  wkbLineStringZ = 1002
+  wkbPolygonZ = 1003
+  wkbMultiPointZ = 1004
+  wkbMultiLineStringZ = 1005
+  wkbMultiPolygonZ = 1006
+  wkbGeometryCollectionZ = 1007
+  wkbCircularStringZ = 1008
+  wkbCompoundCurveZ = 1009
+  wkbCurvePolygonZ = 1010
+  wkbMultiCurveZ = 1011
+  wkbMultiSurfaceZ = 1012
+  wkbCurveZ = 1013
+  wkbSurfaceZ = 1014
+
+  # ISO SQL/MM Part 3.
+  # M-aware types
+  wkbPointM = 2001
+  wkbLineStringM = 2002
+  wkbPolygonM = 2003
+  wkbMultiPointM = 2004
+  wkbMultiLineStringM = 2005
+  wkbMultiPolygonM = 2006
+  wkbGeometryCollectionM = 2007
+  wkbCircularStringM = 2008
+  wkbCompoundCurveM = 2009
+  wkbCurvePolygonM = 2010
+  wkbMultiCurveM = 2011
+  wkbMultiSurfaceM = 2012
+  wkbCurveM = 2013
+  wkbSurfaceM = 2014
+
+  # ISO SQL/MM Part 3.
+  # ZM-aware types ... Meshes.jl doesn't generally support this?
+  wkbPointZM = 3001
+  wkbLineStringZM = 3002
+  wkbPolygonZM = 3003
+  wkbMultiPointZM = 3004
+  wkbMultiLineStringZM = 3005
+  wkbMultiPolygonZM = 3006
+  wkbGeometryCollectionZM = 3007
+  wkbCircularStringZM = 3008
+  wkbCompoundCurveZM = 3009
+  wkbCurvePolygonZM = 3010
+  wkbMultiCurveZM = 3011
+  wkbMultiSurfaceZM = 3012
+  wkbCurveZM = 3013
+  wkbSurfaceZM = 3014
+
+  # 2.5D extension as per 99-402
+  # https://lists.osgeo.org/pipermail/postgis-devel/2004-December/000702.html
+
+  # Julia throws InexactError: convert(Int32, 0x8000000x)
+  wkbPoint25D = -2147483647 # 0x80000001
+  wkbLineString25D = -2147483646 # 0x80000002
+  wkbPolygon25D = -2147483645 # 0x80000003
+  wkbMultiPoint25D = -2147483644 # 0x80000004
+  wkbMultiLineString25D = -2147483643 # 0x80000005
+  wkbMultiPolygon25D = -2147483642 # 0x80000006
+  wkbGeometryCollection25D = -2147483641 # 0x80000007
+end
+
+include("gpkg/read.jl")
+include("gpkg/write.jl")