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, 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="GeoIO.jl", save="ArchGDAL.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/gpkg.jl

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

src/extra/gpkg/read.jl

@@ -0,0 +1,159 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+function gpkgread(fname; layer=1)
+  db = gpkgdatabase(fname)
+  table, geoms = gpkgextract(db; layer)
+  DBInterface.close!(db)
+  georef(table, geoms)
+end
+
+function gpkgdatabase(fname)
+  # connect to SQLite database on disk
+  db = SQLite.DB(fname)
+
+  # According to https://www.geopackage.org/spec/#r6
+  # PRAGMA integrity_check returns a single row with the value 'ok'
+  if first(DBInterface.execute(db, "PRAGMA integrity_check;")).integrity_check != "ok"
+    throw(ErrorException("database integrity at risk"))
+  end
+  # According to https://www.geopackage.org/spec/#r7
+  # PRAGMA foreign_key_check (w/ no parameter value) returns an empty result set
+  if !(isempty(DBInterface.execute(db, "PRAGMA foreign_key_check;")))
+    throw(ErrorException("foreign key violation(s)"))
+  end
+
+  # According to https://www.geopackage.org/spec/#r10
+  # A GeoPackage SHALL include a 'gpkg_spatial_ref_sys' table
+  if isnothing(SQLite.tableinfo(db, "gpkg_spatial_ref_sys"))
+    throw(ErrorException("missing required metadata tables in the GeoPackage SQL database"))
+  end
+  # According to https://www.geopackage.org/spec/#r13
+  # A GeoPackage SHALL include a 'gpkg_contents` table
+  if isnothing(SQLite.tableinfo(db, "gpkg_contents"))
+    throw(ErrorException("missing required metadata tables in the GeoPackage SQL database"))
+  end
+  db
+end
+
+function gpkgextract(db; layer=1)
+  # get the feature table or Any[] returned in sqlite query results
+  metadata = first(
+    DBInterface.execute(
+      db,
+      """
+      SELECT g.table_name AS tablename, g.column_name AS geomcolumn, 
+      c.srs_id AS srsid, srs.organization AS org, srs.organization_coordsys_id AS code
+      FROM gpkg_geometry_columns g, gpkg_spatial_ref_sys srs
+      """ *
+      # According to https://www.geopackage.org/spec/#r24
+      # The column_name column value in a gpkg_geometry_columns row SHALL be the name of a column in the table or view specified by the table_name column value for that row.
+      """
+      JOIN gpkg_contents c ON ( g.table_name = c.table_name )
+      """ *
+      # According to https://www.geopackage.org/spec/#r23
+      # gpkg_geometry_columns table_name column SHALL reference values in the gpkg_contents table_name column for rows with a data_type of 'features'
+      """
+      WHERE c.data_type = 'features'
+      """ *
+      # According to https://www.geopackage.org/spec/#r146
+      # The srs_id value in a gpkg_geometry_columns table row SHALL match the srs_id column value from the corresponding row in the gpkg_contents table.
+      """
+      AND g.srs_id = c.srs_id
+      LIMIT 1 OFFSET ($layer-1);
+      """
+    )
+  )
+
+  # According to https://www.geopackage.org/spec/#r33, feature table geometry columns
+  # SHALL contain geometries with the srs_id specified for the column by the gpkg_geometry_columns table srs_id column value.
+  org = metadata.org
+  code = metadata.code
+  srsid = metadata.srsid
+  if srsid == 0
+    crs = LatLon{NoDatum}
+  elseif srsid == 4326
+    crs = LatLon{WGS84Latest}
+  elseif srsid == -1
+    crs = Cartesian{NoDatum}
+  else
+    if org == "EPSG"
+      crs = CoordRefSystems.get(EPSG{code})
+    elseif org == "ESRI"
+      crs = CoordRefSystems.get(ESRI{code})
+    end
+  end
+
+  # According to https://www.geopackage.org/spec/#r14
+  # The table_name column value in a gpkg_contents table row 
+  # SHALL contain the name of a SQLite table or view.
+  tablename = metadata.tablename
+  geomcolumn = metadata.geomcolumn
+  tableinfo = SQLite.tableinfo(db, tablename)
+
+  # "pk" (either zero for columns that are not part of the primary key, or the 1-based index of the column within the primary key)
+  columns = [name for (name, pk) in zip(tableinfo.name, tableinfo.pk) if pk == 0]
+  gpkgbinary = DBInterface.execute(db, "SELECT  $(join(columns, ',')) FROM $tablename;")
+  table = map(gpkgbinary) do row
+    # According to https://www.geopackage.org/spec/#r30
+    # A feature table or view SHALL have only one geometry column.
+    values = [(key, getproperty(row, key)) for key in keys(row) if key != Symbol(geomcolumn)]
+
+    # create IOBuffer and seek geometry binary data
+    buff = wkbgeombuffer(row, geomcolumn)
+
+    geom = wkb2geom(buff, crs)
+
+    # returns a tuple of the corresponding aspatial attributes and the geometries for each row in the feature table
+    NamedTuple(values), geom
+  end
+  # separate aspatial attributes into its own vector
+  aspatial = getindex.(table, 1)
+
+  # check if type of aspatial elements generated are calls to the macro for declaring `NamedTuple` types
+  eltype(aspatial) != @NamedTuple{} || return nothing, getindex.(table, 2)
+  # aspatial attributes and corresponding geometries
+  aspatial, getindex.(table, 2)
+end
+
+function wkbgeombuffer(row, geomcolumn)
+  # get the column of SQL Geometry Binary specified by gpkg_geometry_columns table in column_name field
+  buff = IOBuffer(getproperty(row, Symbol(geomcolumn)))
+
+  # According to https://www.geopackage.org/spec/#r19
+  # A GeoPackage SHALL store feature table geometries in SQL BLOBs using the Standard GeoPackageBinary format
+  # check the GeoPackageBinaryHeader for the first byte[2] to be 'GP' in ASCII
+  read(buff, UInt16) == 0x5047 || @warn "Missing magic 'GP' string in GPkgBinaryGeometry"
+
+  # byte[1] version: 8-bit unsigned integer, 0 = version 1
+  read(buff, UInt8)
+
+  # bit layout of GeoPackageBinary flags byte
+  # https://www.geopackage.org/spec/#flags_layout
+  # ---------------------------------------
+  # bit # 7 # 6 # 5 # 4 # 3 # 2 # 1 # 0 #
+  # use # R # R # X # Y # E # E # E # B #
+  # ---------------------------------------
+  # R: reserved for future use; set to 0
+  # X: GeoPackageBinary type
+  # Y: empty geometry flag
+  # E: envelope contents indicator code (3-bit unsigned integer)
+  # B: byte order for SRS_ID and envelope values in header
+  flag = read(buff, UInt8)
+
+  # 0x07 is a 3-bit mask 0x00001110
+  # left-shift moves the 3-bit mask by one to align with E bits in flag layout
+  # bitwise AND operation isolates the E bits
+  # right-shift moves the E bits by one to align with the least significant bits
+  # results in a 3-bit unsigned integer
+  envelope = (flag & (0x07 << 1)) >> 1
+
+  # calculate GeoPackageBinaryHeader size in byte stream given extent of envelope:
+  # envelope is [minx, maxx, miny, maxy, minz, maxz], 48 bytes or envelope is [minx, maxx, miny, maxy], 32 bytes or no envelope, 0 bytes
+  # byte[2] magic + byte[1] version + byte[1] flag + byte[4] srs_id + byte[(8*2)×(x,y{,z})] envelope
+  skiplen = iszero(envelope) ? 4 : 4 + 8 * 2 * (envelope + 1)
+
+  # Skip reading the double[] envelope and start reading Well-Known Binary geometry 
+  skip(buff, skiplen)
+end

src/extra/gpkg/wkb.jl

@@ -0,0 +1,106 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+function wkb2geom(buff, crs)
+  byteswap = isone(read(buff, UInt8)) ? ltoh : ntoh
+  wkbtype = read(buff, UInt32)
+  # Input variants of WKB supported are standard, extended, and ISO WKB geometry with Z dimensions (M/ZM not supported)
+  # SQL/MM Part 3 and SFSQL 1.2 use offsets of 1000 (Z), 2000 (M) and 3000 (ZM) 
+  # to indicate the present of higher dimensional coordinates in a WKB geometry
+  if wkbtype >= 1001 && wkbtype <= 1007
+    # the SFSQL 1.2 offset of 1000 (Z) is present and subtracting a round number of 1000 gives the standard WKB type
+    wkbtype -= UInt32(1000)
+    # 99-402 was a short-lived extension to SFSQL 1.1 that used a high-bit flag to indicate the presence of Z coordinates in a WKB geometry
+    # the high-bit flag 0x80000000 for Z (or 0x40000000 for M) is set and masking it off gives the standard WKB type
+  elseif wkbtype > 0x80000000
+    # the SFSQL 1.1  high-bit flag 0x80000000 (Z) is present and removing the flag reveals the standard WKB type
+    wkbtype -= 0x80000000
+  end
+  if wkbtype <= 3
+    # 0 - 3 [Geometry, Point, Linestring, Polygon]
+    wkb2single(buff, crs, wkbtype, byteswap)
+  else
+    # 4 - 7 [MultiPoint, MultiLinestring, MultiPolygon, GeometryCollection]
+    wkb2multi(buff, crs, byteswap)
+  end
+end
+
+# read single features from Well-Known Binary IO Buffer and return Concrete Geometry
+function wkb2single(buff, crs, wkbtype, byteswap)
+  if wkbtype == 1
+    wkb2point(buff, crs, byteswap)
+  elseif wkbtype == 2
+    wkb2chain(buff, crs, byteswap)
+  elseif wkbtype == 3
+    wkb2poly(buff, crs, byteswap)
+  else
+    error("Unsupported WKB Geometry Type: $wkbtype")
+  end
+end
+
+function wkb2point(buff, crs, byteswap)
+  coordinates = wkb2coords(buff, crs, byteswap)
+  Point(referencecoords(coordinates, crs))
+end
+
+function wkb2coords(buff, crs, byteswap)
+  if CoordRefSystems.ncoords(crs) == 2
+    x = byteswap(read(buff, Float64))
+    y = byteswap(read(buff, Float64))
+    return (x, y)
+  elseif CoordRefSystems.ncoords(crs) == 3
+    x = byteswap(read(buff, Float64))
+    y = byteswap(read(buff, Float64))
+    z = byteswap(read(buff, Float64))
+    return (x, y, z)
+  end
+end
+
+function referencecoords(coordinates, crs)
+  if crs <: LatLon
+    crs(coordinates[2], coordinates[1])
+  elseif crs <: LatLonAlt
+    crs(coordinates[2], coordinates[1], coordinates[3])
+  else
+    crs(coordinates...)
+  end
+end
+
+function wkb2points(buff, npoints, crs, byteswap)
+  map(1:npoints) do _
+    coordinates = wkb2coords(buff, crs, byteswap)
+    Point(referencecoords(coordinates, crs))
+  end
+end
+
+function wkb2chain(buff, crs, byteswap)
+  npoints = byteswap(read(buff, UInt32))
+  chain = wkb2points(buff, npoints, crs, byteswap)
+  if length(chain) >= 2 && first(chain) == last(chain)
+    Ring(chain[1:(end - 1)])
+  elseif length(chain) >= 2
+    Rope(chain)
+  else
+    # single point or closed single point
+    Ring(chain)
+  end
+end
+
+function wkb2poly(buff, crs, byteswap)
+  nrings = byteswap(read(buff, UInt32))
+  rings = map(1:nrings) do _
+    wkb2chain(buff, crs, byteswap)
+  end
+  PolyArea(rings)
+end
+
+function wkb2multi(buff, crs, byteswap)
+  ngeoms = byteswap(read(buff, UInt32))
+  geoms = map(1:ngeoms) do _
+    wkbbswap = isone(read(buff, UInt8)) ? ltoh : ntoh
+    wkbtype = read(buff, UInt32)
+    wkb2single(buff, crs, wkbtype, wkbbswap)
+  end
+  Multi(geoms)
+end

src/load.jl

@@ -156,6 +156,11 @@ function load(fname; repair=true, layer=1, lenunit=nothing, numtype=Float64, kwa
     return cdmread(fname; kwargs...)
   end
 
+  # GeoPackage formats
+  if endswith(fname, ".gpkg")
+    return gpkgread(fname; layer, kwargs...)
+  end
+
   # GIS formats
   gisread(fname; repair, layer, numtype, kwargs...)
 end