Add support for loading CRS in GeoParquet

Project.toml

@@ -20,8 +20,10 @@ GeoTables = "e502b557-6362-48c1-8219-d30d308dcdb0"
 GslibIO = "4610876b-9b01-57c8-9ad9-06315f1a66a5"
 ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534"
 ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19"
+JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
 Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
+Parquet2 = "98572fba-bba0-415d-956f-fa77e587d26d"
 PlyIO = "42171d58-473b-503a-8d5f-782019eb09ec"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
@@ -51,8 +53,10 @@ GeoTables = "1.24.0"
 GslibIO = "1.5"
 ImageCore = "0.10"
 ImageIO = "0.6"
+JSON3 = "1.14"
 Meshes = "0.50 - 0.51"
 NCDatasets = "0.13, 0.14"
+Parquet2 = "0.2"
 PlyIO = "1.1"
 PrecompileTools = "1.2"
 PrettyTables = "2.2"

src/GeoIO.jl

@@ -18,6 +18,10 @@ using TransformsBase: Identity, →
 using ImageCore: channelview
 using Unitful: m
 
+# for GeoPackage utils
+import Parquet2
+import JSON3
+
 # image formats
 import FileIO
 

src/conversion.jl

@@ -64,6 +64,7 @@ crstype(crs::GFT.EPSG, _) = CoordRefSystems.get(EPSG{GFT.val(crs)})
 crstype(crs::GFT.WellKnownText, _) = CoordRefSystems.get(GFT.val(crs))
 crstype(crs::GFT.WellKnownText2, _) = CoordRefSystems.get(GFT.val(crs))
 crstype(crs::GFT.ESRIWellKnownText, _) = CoordRefSystems.get(GFT.val(crs))
+crstype(crs::GFT.ProjJSON, _) = CoordRefSystems.get(projsoncode(GFT.val(crs)))
 crstype(_, geom) = Cartesian{NoDatum,GI.is3d(geom) ? 3 : 2}
 
 function topoint(geom, CRS)

src/load.jl

@@ -147,7 +147,11 @@ function load(fname; repair=true, layer=0, lenunit=m, numbertype=Float64, kwargs
   end
 
   # construct geotable
-  geotable = asgeotable(table)
+  geotable = if endswith(fname, ".parquet")
+    asgeotable(table, gpqcrs(fname))
+  else
+    asgeotable(table)
+  end
 
   # repair pipeline
   pipeline = if repair

src/save.jl

@@ -164,7 +164,13 @@ function save(fname, geotable; kwargs...)
     end
     GJS.write(fname, geotable |> proj; kwargs...)
   elseif endswith(fname, ".parquet")
-    GPQ.write(fname, geotable, (:geometry,); kwargs...)
+    gpqcrs = try
+      code = CoordRefSystems.code(crs(domain(geotable)))
+      GFT.ProjJSON(projjson(code))
+    catch
+      nothing
+    end
+    GPQ.write(fname, geotable, (:geometry,), gpqcrs; kwargs...)
   else # fallback to GDAL
     agwrite(fname, geotable; kwargs...)
   end

src/utils.jl

@@ -6,8 +6,7 @@
 const Met{T} = Quantity{T,u"𝐋",typeof(u"m")}
 const Deg{T} = Quantity{T,NoDims,typeof(u"°")}
 
-function asgeotable(table)
-  crs = GI.crs(table)
+function asgeotable(table, crs=GI.crs(table))
   cols = Tables.columns(table)
   names = Tables.columnnames(cols)
   gcol = geomcolumn(names)
@@ -60,3 +59,28 @@ spatialref(code) = AG.importUserInput(codestring(code))
 
 codestring(::Type{EPSG{Code}}) where {Code} = "EPSG:$Code"
 codestring(::Type{ESRI{Code}}) where {Code} = "ESRI:$Code"
+
+function projjson(code)
+  spref = spatialref(code)
+  wktptr = Ref{Cstring}()
+  options = ["MULTILINE=NO"]
+  GDAL.osrexporttoprojjson(spref, wktptr, options)
+  str = unsafe_string(wktptr[])
+  JSON3.read(str, Dict)
+end
+
+projsoncode(crs::String) = projsoncode(JSON3.read(crs, Dict))
+
+function projsoncode(crs::Dict)
+  code = parse(Int, crs["id"]["code"])
+  type = crs["id"]["authority"]
+  type == "EPSG" ? EPSG{code} : ESRI{code}
+end
+
+function gpqcrs(fname)
+  ds = Parquet2.Dataset(fname)
+  meta = Parquet2.metadata(ds)["geo"]
+  json = JSON3.read(meta, Dict)
+  crs = json["columns"]["geometry"]["crs"]
+  isnothing(crs) ? nothing : GFT.ProjJSON(crs)
+end