add convert methods and reproject for GeoFormatTypes

Project.toml

@@ -9,12 +9,14 @@ version = "0.3.0"
 DataStreams = "9a8bc11e-79be-5b39-94d7-1ccc349a1a85"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 GDAL = "add2ef01-049f-52c4-9ee2-e494f65e021a"
+GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 
 [compat]
 DataStreams = "0.4.2"
 GDAL = "1.0.2"
 GeoInterface = "0.4"
+GeoFormatTypes = "0.1.1"
 julia = "1"
 
 [extras]

src/ArchGDAL.jl

@@ -1,9 +1,13 @@
 module ArchGDAL
 
-    import GDAL, GeoInterface
+    import GDAL, GeoInterface, GeoFormatTypes
     import DataStreams: Data
     import GeoInterface: coordinates, geotype
+    import Base: convert
+
     using Dates
+
+    const GFT = GeoFormatTypes
 
     include("utils.jl")
     include("types.jl")

src/ogr/geometry.jl

@@ -1,3 +1,32 @@
+
+# This is friendly type-pyracy on `convert` for GeoFormatTypes types
+Base.convert(::Type{<:Geometry}, source::GFT.AbstractWellKnownText) = 
+    fromWKT(GFT.val(source))
+Base.convert(::Type{<:Geometry}, source::GFT.WellKnownBinary) = 
+    fromWKB(GFT.val(source))
+Base.convert(::Type{<:Geometry}, source::GFT.GeoJSON) = 
+    fromJSON(GFT.val(source))
+Base.convert(::Type{<:Geometry}, source::GFT.GML) = 
+    fromGML(GFT.val(source))
+
+Base.convert(::Type{<:GFT.AbstractWellKnownText}, source::AbstractGeometry) = 
+    GFT.WellKnownText(GFT.Geom(), toWKT(source))
+Base.convert(::Type{<:GFT.WellKnownBinary}, source::AbstractGeometry) = 
+    GFT.WellKnownBinary(GFT.Geom(), toWKB(source))
+Base.convert(::Type{<:GFT.GeoJSON}, source::AbstractGeometry) = 
+    GFT.GeoJSON(toJSON(source))
+Base.convert(::Type{<:GFT.GML}, source::AbstractGeometry) = 
+    GFT.GML(GFT.Geom(), toGML(source))
+Base.convert(::Type{<:GFT.KML}, source::AbstractGeometry) = 
+    GFT.KML(toKML(source))
+
+# Convert to Geometry, then to the target format. There may be
+# a more optimal method, but this is a simple, high-level first pass.
+# The Geom trait is needed to separate out convert for CRS for WellKnownText
+# and GML, which may contain both. It is handled in GeoFormatTypes.
+Base.convert(format::Type{<:GFT.GeoFormat}, mode::GFT.Geom, source::GFT.GeoFormat) =
+    convert(format, convert(Geometry, source))
+
 """
     fromWKB(data)
 

src/spatialref.jl

@@ -1,3 +1,85 @@
+"""
+Import format and run function `f` on the result.
+
+The name could be improved?
+"""
+importcrs(f::Function, sourceformat::GFT.GeoFormat) = begin
+    obj = importcrs(sourceformat)
+    try f(obj) finally destroy(obj) end
+end
+importcrs(x::GFT.EPSG) = importEPSG(GFT.val(x))
+importcrs(x::GFT.AbstractWellKnownText) = importWKT(GFT.val(x))
+importcrs(x::GFT.ESRIWellKnownText) = importESRI(GFT.val(x))
+importcrs(x::GFT.ProjString) = importPROJ4(GFT.val(x))
+importcrs(x::GFT.GML) = importXML(GFT.val(x))
+importcrs(x::GFT.KML) = importEPSG(EPSG(4327))
+importcrs(x::GFT.GeoJSON) = importEPSG(EPSG(4327)) # TODO check GeoJSON standard
+
+Base.convert(target::Type{GFT.GeoFormat}, mode::GFT.CRS, input::GFT.GeoFormat) =
+    unsafe_convert(target, mode, importcrs(input))
+
+unsafe_convert(::Type{GFT.CoordSys}, mode::GFT.CRS, srs::Ref) = GFT.CoordSys(toMICoordSys(srs))
+unsafe_convert(::Type{GFT.ProjString}, mode::GFT.CRS, srs::Ref) = GFT.ProjString(toPROJ4(srs))
+unsafe_convert(::Type{GFT.WellKnownText}, mode::GFT.CRS, srs::Ref) = GFT.WellKnownText(toWKT(srs))
+unsafe_convert(::Type{GFT.ESRIWellKnownText}, mode::GFT.CRS, srs::Ref) =
+    GFT.ESRIWellKnownText(toWKT(morphtoESRI!(srs)))
+unsafe_convert(::Type{GFT.GML}, mode::GFT.CRS, srs::Ref) = GFT.GML(toXML(srs))
+
+# These should be better integrated with geometry packages or follow some standard
+const ReprojectCoord = Union{<:NTuple{2,<:Number},<:NTuple{3,<:Number},AbstractVector{<:Number}}
+
+"""
+    reproject(points, sourceproj::GeoFormat, destproj::GeoFormat)
+
+Reproject points to a different coordinate reference system and/or format.
+
+## Arguments
+-`points`: Vector of Geometry points
+-`crs`: The current coordinate reference system, using any GeoFormat
+-`targetcrs`: The coordinate reference system to transform to, using any CRS capable GeoFormat
+
+## Example
+
+reproject([(118, 34), (119, 35)], EPSG(2025), EPSG(4326))
+"""
+reproject(x, crs::GFT.GeoFormat, targetcrs::Nothing) = x
+
+# Vector/Tuple coordinate(s)
+reproject(coord::ReprojectCoord, crs::GFT.GeoFormat, targetcrs::GFT.GeoFormat) =
+    GeoInterface.coordinates(reproject(createpoint(coord...), crs, targetcrs))
+reproject(coords::AbstractArray{<:ReprojectCoord}, crs::GFT.GeoFormat, targetcrs::GFT.GeoFormat) =
+    GeoInterface.coordinates.(reproject([createpoint(c...) for c in coords], crs, targetcrs))
+
+# GeoFormat
+reproject(geom::GFT.GeoFormat, crs::GFT.GeoFormat, targetcrs::GFT.GeoFormat) =
+    convert(typeof(geom), reproject(convert(Geometry, geom), crs, targetcrs))
+
+# Geometries
+reproject(geom::AbstractGeometry, crs::GFT.GeoFormat, targetcrs::GFT.GeoFormat) =
+    crs2transform(crs, targetcrs) do transform
+        transform!(geom, transform)
+    end
+reproject(geoms::AbstractArray{<:AbstractGeometry}, crs::GFT.GeoFormat, targetcrs::GFT.GeoFormat) =
+    crs2transform(crs, targetcrs) do transform
+        transform!.(geoms, Ref(transform))
+    end
+
+"""
+    crs2transform(f, crs::GeoFormat, targetcrs::GeoFormat)
+
+Run the function on a coord transform generated from
+two crs definitions. These can be in any GeoFormatTypes format
+that holds a coordinate reference system.
+"""
+crs2transform(f, crs::GFT.GeoFormat, targetcrs::GFT.GeoFormat) =
+    importcrs(crs) do crs_ref
+        importcrs(targetcrs) do targetcrs_ref
+            createcoordtrans(crs_ref, targetcrs_ref) do transform
+                f(transform)
+            end
+        end
+    end
+
 """
     newspatialref(wkt::AbstractString = "")
 

test/test_cookbook_projection.jl

@@ -1,6 +1,8 @@
 using Test
-import GeoInterface
-import ArchGDAL; const AG = ArchGDAL
+import GeoInterface, GeoFormatTypes, ArchGDAL 
+const AG = ArchGDAL
+const GFT = GeoFormatTypes
+const GI = GeoFormatTypes
 
 @testset "Reproject a Geometry" begin
     @testset "Method 1" begin
@@ -25,6 +27,22 @@ import ArchGDAL; const AG = ArchGDAL
                 @test zs ≈ [0.0, 0.0]
         end end end
     end
+
+    @testset "Use reproject" begin
+        @testset "reciprocal reprojection of wkt" begin
+            wktpoint = GFT.WellKnownText(GFT.Geom(), "POINT (1120351.57 741921.42)")
+            result = GFT.WellKnownText(GFT.Geom(), "POINT (47.3488013802885 -122.598135130878)")
+            @test AG.reproject(wktpoint, GFT.EPSG(2927), GFT.EPSG(4326)) == result
+            @test convert(AG.Geometry, AG.reproject(result, GFT.EPSG(4326), GFT.EPSG(2927))) |> 
+                GeoInterface.coordinates ≈ [1120351.57, 741921.42]
+        end
+        @testset "reproject vector, vector of vector, or tuple" begin
+            coord = [1120351.57, 741921.42]
+            @test AG.reproject(coord, GFT.EPSG(2927), GFT.EPSG(4326)) ≈ [47.348801, -122.598135]
+            @test AG.reproject([coord], GFT.EPSG(2927), GFT.EPSG(4326)) ≈ [[47.348801, -122.598135]]
+            coord = (1120351.57, 741921.42)
+            @test AG.reproject(coord, GFT.EPSG(2927), GFT.EPSG(4326)) ≈ [47.348801, -122.598135]
+    end
 end
 
 @testset "Get Projection" begin

test/test_geometry.jl

@@ -1,5 +1,8 @@
 using Test
-import GeoInterface, GDAL, ArchGDAL; const AG = ArchGDAL
+import GeoInterface, GeoFormatTypes, GDAL, ArchGDAL;
+const AG = ArchGDAL
+const GFT = GeoFormatTypes
+
 
 @testset "Incomplete GeoInterface geometries" begin
     @test_logs (:warn, "unknown geometry type") GeoInterface.geotype(AG.creategeom(GDAL.wkbCircularString))
@@ -95,6 +98,25 @@ end
     @test AG.toJSON(point) == "{ \"type\": \"Point\", \"coordinates\": [ 100.0, 70.0 ] }"
 end
 
+
+
+@testset "convert point format" begin
+    point = AG.createpoint(100, 70)
+    json = convert(GFT.GeoJSON, point)
+    kml = convert(GFT.KML, point)
+    gml = convert(GFT.GML, point)
+    wkb = convert(GFT.WellKnownBinary, point) 
+    wkt = convert(GFT.WellKnownText, point) 
+    @test GFT.val(json) == AG.toJSON(point)
+    @test GFT.val(kml) == AG.toKML(point)
+    @test GFT.val(gml) == AG.toGML(point) # This isn't actually tested above
+    @test GFT.val(wkb) == AG.toWKB(point)
+    @test GFT.val(wkt) == AG.toWKT(point)
+    @test convert(GFT.GeoJSON, json) == json
+    @test convert(GFT.GeoJSON, wkb) == convert(GFT.GeoJSON, wkt) == convert(GFT.GeoJSON, gml)
+    @test convert(GFT.KML, gml) == convert(GFT.KML, wkt)
+end
+
 @testset "Testing construction of complex geometries" begin
     @test AG.toWKT(AG.createlinestring([1.,2.,3.], [4.,5.,6.])) == "LINESTRING (1 4,2 5,3 6)"
     AG.createlinestring([1.,2.,3.], [4.,5.,6.]) do geom
@@ -105,6 +127,8 @@ end
         @test AG.toWKT(geom) == "LINESTRING (1 4,2 5,3 6)"
         AG.setpoint!(geom, 1, 10, 10)
         @test AG.toWKT(geom) == "LINESTRING (1 4,10 10,3 6)"
+        # Test using convert
+        @test GFT.val(convert(GFT.WellKnownText, geom)) == AG.toWKT(geom)  
     end
     AG.createlinestring([1.,2.,3.], [4.,5.,6.], [7.,8.,9.]) do geom
         @test AG.toWKT(geom) == "LINESTRING (1 4 7,2 5 8,3 6 9)"