Add support for loading CRS in GeoJSON & Update to Meshes.jl v0.46

Project.toml

@@ -8,9 +8,11 @@ ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 CommonDataModel = "1fbeeb36-5f17-413c-809b-666fb144f157"
+CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8"
 GRIBDatasets = "82be9cdb-ee19-4151-bdb3-b400788d9abc"
+GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
 GeoParquet = "e99870d8-ce00-4fdd-aeee-e09192881159"
@@ -36,16 +38,18 @@ ArchGDAL = "0.10"
 CSV = "0.10"
 Colors = "0.12"
 CommonDataModel = "0.2, 0.3"
+CoordRefSystems = "0.9"
 FileIO = "1.16"
 Format = "1.3"
 GRIBDatasets = "0.3"
+GeoFormatTypes = "0.4"
 GeoInterface = "1.0"
 GeoJSON = "0.7, 0.8"
 GeoParquet = "0.1, 0.2"
 GeoTables = "1.21"
-GslibIO = "1.4.10"
+GslibIO = "1.5"
 ImageIO = "0.6"
-Meshes = "0.45"
+Meshes = "0.46"
 NCDatasets = "0.13, 0.14"
 PlyIO = "1.1"
 PrecompileTools = "1.2"

src/GeoIO.jl

@@ -11,6 +11,7 @@ using Unitful
 using GeoTables
 using StaticArrays
 using PrettyTables
+using CoordRefSystems
 using Meshes: SubDomain
 using Format: generate_formatter
 using TransformsBase: →
@@ -45,6 +46,7 @@ import GeoJSON as GJS
 import ArchGDAL as AG
 import GeoParquet as GPQ
 import GeoInterface as GI
+import GeoFormatTypes as GFT
 
 # image extensions
 const IMGEXTS = [".png", ".jpg", ".jpeg"]

src/conversion.jl

@@ -17,23 +17,25 @@
 GI.geomtrait(::MultiRing) = GI.MultiLineStringTrait()
 GI.geomtrait(::MultiPolygon) = GI.MultiPolygonTrait()
 
-GI.ncoord(::GI.PointTrait, p::Point) = embeddim(p)
+GI.ncoord(::GI.PointTrait, p::Point) = CoordRefSystems.ncoords(Meshes.crs(p))
 GI.getcoord(::GI.PointTrait, p::Point) = ustrip.(to(p))
 GI.getcoord(::GI.PointTrait, p::Point, i) = ustrip(to(p)[i])
+GI.getcoord(::GI.PointTrait, p::Point{3,<:LatLon}) = reverse(CoordRefSystems.rawvalues(coords(p)))
+GI.getcoord(trait::GI.PointTrait, p::Point{3,<:LatLon}, i) = GI.getcoord(trait, p)[i]
 
-GI.ncoord(::GI.LineTrait, s::Segment) = embeddim(s)
+GI.ncoord(::GI.LineTrait, s::Segment) = CoordRefSystems.ncoords(Meshes.crs(s))
 GI.ngeom(::GI.LineTrait, s::Segment) = nvertices(s)
 GI.getgeom(::GI.LineTrait, s::Segment, i) = vertex(s, i)
 
-GI.ncoord(::GI.LineStringTrait, c::Chain) = embeddim(c)
+GI.ncoord(::GI.LineStringTrait, c::Chain) = CoordRefSystems.ncoords(Meshes.crs(c))
 GI.ngeom(::GI.LineStringTrait, c::Chain) = nvertices(c) + isclosed(c)
 GI.getgeom(::GI.LineStringTrait, c::Chain, i) = vertex(c, i)
 
-GI.ncoord(::GI.PolygonTrait, p::Polygon) = embeddim(p)
+GI.ncoord(::GI.PolygonTrait, p::Polygon) = CoordRefSystems.ncoords(Meshes.crs(p))
 GI.ngeom(::GI.PolygonTrait, p::Polygon) = length(rings(p))
 GI.getgeom(::GI.PolygonTrait, p::Polygon, i) = rings(p)[i]
 
-GI.ncoord(::GI.AbstractGeometryTrait, m::Multi) = embeddim(m)
+GI.ncoord(::GI.AbstractGeometryTrait, m::Multi) = CoordRefSystems.ncoords(Meshes.crs(m))
 GI.ngeom(::GI.AbstractGeometryTrait, m::Multi) = length(parent(m))
 GI.getgeom(::GI.AbstractGeometryTrait, m::Multi, i) = parent(m)[i]
 
@@ -46,16 +48,20 @@
 # Convert geometries to Meshes.jl types
 # --------------------------------------
 
-function topoints(geom, is3d::Bool)
-  if is3d
-    [Point(GI.x(p), GI.y(p), GI.z(p)) for p in GI.getpoint(geom)]
-  else
-    [Point(GI.x(p), GI.y(p)) for p in GI.getpoint(geom)]
-  end
-end
+getcrs(geom) = getcrs(GI.crs(geom), GI.is3d(geom) ? 3 : 2)
+getcrs(code::GFT.EPSG, _) = CoordRefSystems.get(EPSG{GFT.val(code)})
+getcrs(_, Dim) = Cartesian{NoDatum,Dim}
+
+topoint(geom, ::Type{<:Cartesian{Datum,2}}) where {Datum} = Point(Cartesian{Datum}(GI.x(geom), GI.y(geom)))
+
+topoint(geom, ::Type{<:Cartesian{Datum,3}}) where {Datum} = Point(Cartesian{Datum}(GI.x(geom), GI.y(geom), GI.z(geom)))
+
+topoint(geom, ::Type{<:LatLon{Datum}}) where {Datum} = Point(LatLon{Datum}(GI.y(geom), GI.x(geom)))
 
-function tochain(geom, is3d::Bool)
-  points = topoints(geom, is3d)
+topoints(geom, CRS) = [topoint(p, CRS) for p in GI.getpoint(geom)]
+
+function tochain(geom, CRS)
+  points = topoints(geom, CRS)
   if first(points) == last(points)
     # fix backend issues: https://github.com/JuliaEarth/GeoTables.jl/issues/32
     while first(points) == last(points) && length(points) ≥ 2
@@ -67,9 +73,9 @@
   end
 end
 
-function topolygon(geom, is3d::Bool, fix::Bool)
+function topolygon(geom, CRS, fix::Bool)
   # fix backend issues: https://github.com/JuliaEarth/GeoTables.jl/issues/32
-  toring(g) = close(tochain(g, is3d))
+  toring(g) = close(tochain(g, CRS))
   outer = toring(GI.getexterior(geom))
   if GI.nhole(geom) == 0
     PolyArea(outer; fix)
@@ -79,36 +85,28 @@
   end
 end
 
-function GI.convert(::Type{Point}, ::GI.PointTrait, geom)
-  if GI.is3d(geom)
-    Point(GI.x(geom), GI.y(geom), GI.z(geom))
-  else
-    Point(GI.x(geom), GI.y(geom))
-  end
-end
+GI.convert(::Type{Point}, ::GI.PointTrait, geom) = topoint(geom, getcrs(geom))
 
-GI.convert(::Type{Segment}, ::GI.LineTrait, geom) = Segment(topoints(geom, GI.is3d(geom))...)
+GI.convert(::Type{Segment}, ::GI.LineTrait, geom) = Segment(topoints(geom, getcrs(geom))...)
 
-GI.convert(::Type{Chain}, ::GI.LineStringTrait, geom) = tochain(geom, GI.is3d(geom))
+GI.convert(::Type{Chain}, ::GI.LineStringTrait, geom) = tochain(geom, getcrs(geom))
 
 GI.convert(::Type{Polygon}, trait::GI.PolygonTrait, geom) = _convert_with_fix(trait, geom, true)
 
-function GI.convert(::Type{Multi}, ::GI.MultiPointTrait, geom)
-  Multi(topoints(geom, GI.is3d(geom)))
-end
+GI.convert(::Type{Multi}, ::GI.MultiPointTrait, geom) = Multi(topoints(geom, getcrs(geom)))
 
 function GI.convert(::Type{Multi}, ::GI.MultiLineStringTrait, geom)
-  is3d = GI.is3d(geom)
-  Multi([tochain(g, is3d) for g in GI.getgeom(geom)])
+  CRS = getcrs(geom)
+  Multi([tochain(g, CRS) for g in GI.getgeom(geom)])
 end
 
 GI.convert(::Type{Multi}, trait::GI.MultiPolygonTrait, geom) = _convert_with_fix(trait, geom, true)
 
-_convert_with_fix(::GI.PolygonTrait, geom, fix) = topolygon(geom, GI.is3d(geom), fix)
+_convert_with_fix(::GI.PolygonTrait, geom, fix) = topolygon(geom, getcrs(geom), fix)
 
 function _convert_with_fix(::GI.MultiPolygonTrait, geom, fix)
-  is3d = GI.is3d(geom)
-  Multi([topolygon(g, is3d, fix) for g in GI.getgeom(geom)])
+  CRS = getcrs(geom)
+  Multi([topolygon(g, CRS, fix) for g in GI.getgeom(geom)])
 end
 
 # -----------------------------------------

test/Project.toml

@@ -1,6 +1,7 @@
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
+CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"

test/convert.jl

@@ -53,8 +53,8 @@
   @test GeoIO.geom2meshes(chain) == Ring((2.2, 2.2))
 
   # GeoJSON.jl
-  points = Point.([(0.0f0, 0.0f0), (2.2f0, 2.2f0), (0.5f0, 2.0f0)])
-  outer = Point.([(0.0f0, 0.0f0), (2.2f0, 2.2f0), (0.5f0, 2.0f0)])
+  points = Point.([LatLon(0.0f0, 0.0f0), LatLon(2.2f0, 2.2f0), LatLon(2.0f0, 0.5f0)])
+  outer = Point.([LatLon(0.0f0, 0.0f0), LatLon(2.2f0, 2.2f0), LatLon(2.0f0, 0.5f0)])
   point = GJS.read("""{"type":"Point","coordinates":[1,1]}""")
   chain = GJS.read("""{"type":"LineString","coordinates":[[0,0],[2.2,2.2],[0.5,2]]}""")
   poly = GJS.read("""{"type":"Polygon","coordinates":[[[0,0],[2.2,2.2],[0.5,2],[0,0]]]}""")
@@ -64,13 +64,13 @@
   multipoly = GJS.read(
     """{"type":"MultiPolygon","coordinates":[[[[0,0],[2.2,2.2],[0.5,2],[0,0]]],[[[0,0],[2.2,2.2],[0.5,2],[0,0]]]]}"""
   )
-  @test GeoIO.geom2meshes(point) == Point(1.0f0, 1.0f0)
+  @test GeoIO.geom2meshes(point) == Point(LatLon(1.0f0, 1.0f0))
   @test GeoIO.geom2meshes(chain) == Rope(points)
   @test GeoIO.geom2meshes(poly) == PolyArea(outer)
   @test GeoIO.geom2meshes(multipoint) == Multi(points)
   @test GeoIO.geom2meshes(multichain) == Multi([Rope(points), Rope(points)])
   @test GeoIO.geom2meshes(multipoly) == Multi([PolyArea(outer), PolyArea(outer)])
   # degenerate chain with 2 equal points
   chain = GJS.read("""{"type":"LineString","coordinates":[[2.2,2.2],[2.2,2.2]]}""")
-  @test GeoIO.geom2meshes(chain) == Ring(Point(2.2f0, 2.2f0))
+  @test GeoIO.geom2meshes(chain) == Ring(Point(LatLon(2.2f0, 2.2f0)))
 end

test/gis.jl

@@ -31,7 +31,13 @@
     # compare domain and values
     d1 = domain(gt1)
     d2 = domain(gt2)
-    @test _isequal(d1, d2)
+    fmt1 = last(splitext(f1))
+    fmt2 = fmt
+    if fmt1 != fmt2 && (fmt1 == ".geojson" || fmt2 == ".geojson")
+      @test !_isequal(d1, d2)
+    else
+      @test _isequal(d1, d2)
+    end
     t1 = values(gt1)
     t2 = values(gt2)
     c1 = Tables.columns(t1)

test/noattrs.jl

@@ -9,13 +9,6 @@
   @test all(ismissing, gtb2[:, 1])
   @test gtb2.geometry == gtb1.geometry
 
-  # GeoJSON
-  file = joinpath(savedir, "noattribs.geojson")
-  GeoIO.save(file, gtb1)
-  gtb2 = GeoIO.load(file, numbertype=Float64)
-  @test isnothing(values(gtb2))
-  @test gtb2 == gtb1
-
   # GeoParquet
   file = joinpath(savedir, "noattribs.parquet")
   GeoIO.save(file, gtb1)
@@ -61,4 +54,13 @@
   gtb2 = GeoIO.load(file)
   @test isnothing(values(gtb2))
   @test gtb2 == gtb1
+
+  # GeoJSON
+  pset = PointSet([Point(LatLon(30, 60)), Point(LatLon(30, 61)), Point(LatLon(31, 60))])
+  gtb1 = georef(nothing, pset)
+  file = joinpath(savedir, "noattribs.geojson")
+  GeoIO.save(file, gtb1)
+  gtb2 = GeoIO.load(file, numbertype=Float64)
+  @test isnothing(values(gtb2))
+  @test gtb2 == gtb1
 end

test/runtests.jl

@@ -2,6 +2,7 @@ using GeoIO
 using Tables
 using Meshes
 using GeoTables
+using CoordRefSystems
 using Test, Random
 using ReferenceTests
 using Colors