WIP: Geo format type wrappers and generic methods (#97)

* add convert methods and reproject for GeoFormatTypes

* bump geoformattypes version

* bugfix tests

* move conversions to convert.jl

* update convert

* tweaks to finish PR

* test importCRS

* bump GeoFormatTypes version to 0.3

* clean up docs

Project.toml

@@ -9,12 +9,14 @@ version = "0.3.1"
 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, 0.5"
+GeoFormatTypes = "0.3"
 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")
@@ -28,6 +32,7 @@ module ArchGDAL
     include("base/display.jl")
     include("datastreams.jl")
     include("geointerface.jl")
+    include("convert.jl")
 
     mutable struct DriverManager
         function DriverManager()

src/context.jl

@@ -202,7 +202,7 @@ for gdalfunc in (
         :gdalbuildvrt, :gdaldem, :gdalgrid, :gdalnearblack, :gdalrasterize,
         :gdaltranslate, :gdalvectortranslate, :gdalwarp, :getband,
         :getcolortable, :getfeature, :getgeom, :getlayer, :getmaskband,
-        :getoverview, :getpart, :getspatialref, :intersection, :importEPSG,
+        :getoverview, :getpart, :getspatialref, :importCRS, :intersection, :importEPSG,
         :importEPSGA, :importESRI, :importPROJ4, :importWKT, :importXML,
         :importURL, :lineargeom, :newspatialref, :nextfeature, :pointalongline,
         :pointonsurface, :polygonfromedges, :polygonize, :read, :sampleoverview,

src/convert.jl

@@ -0,0 +1,64 @@
+# This file contains friendly type-pyracy on `convert` for GeoFormatTypes.jl types
+
+"""
+    convert(::Type{<:AbstractGeometry}, mode::Geom source::GeoFormat)
+
+Convert a GeoFromat object to Geometry, then to the target format.
+The Geom trait is needed to separate out convert for CRS for WellKnownText
+and GML, which may contain both.
+"""
+Base.convert(target::Type{<:GFT.GeoFormat}, mode::Union{GFT.FormatMode,Type{GFT.FormatMode}}, 
+             source::GFT.GeoFormat) =
+    convert(target, convert(AbstractGeometry, source))
+
+"""
+    convert(::Type{<:AbstractGeometry}, source::GeoFormat)
+
+Convert `GeoFormat` geometry data to an ArchGDAL `Geometry` type
+"""
+Base.convert(::Type{<:AbstractGeometry}, source::GFT.AbstractWellKnownText) = 
+    fromWKT(GFT.val(source))
+Base.convert(::Type{<:AbstractGeometry}, source::GFT.WellKnownBinary) = 
+    fromWKB(GFT.val(source))
+Base.convert(::Type{<:AbstractGeometry}, source::GFT.GeoJSON) = 
+    fromJSON(GFT.val(source))
+Base.convert(::Type{<:AbstractGeometry}, source::GFT.GML) = 
+    fromGML(GFT.val(source))
+
+"""
+    convert(::Type{<:AbstractGeometry}, source::GeoFormat)
+
+Convert `AbstractGeometry` data to any gemoetry `GeoFormat` 
+"""
+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(target::Type{GeoFormat}, mode::CRS, source::GeoFormat)
+
+Convert `GeoFormat` crs data to another another `GeoFormat` crs type.
+"""
+Base.convert(target::Type{<:GFT.GeoFormat}, mode::Union{GFT.CRS,Type{GFT.CRS}}, 
+             source::GFT.GeoFormat) =
+    unsafe_convertcrs(target, importCRS(source))
+
+unsafe_convertcrs(::Type{<:GFT.CoordSys}, crsref) = 
+    GFT.CoordSys(toMICoordSys(crsref))
+unsafe_convertcrs(::Type{<:GFT.ProjString}, crsref) = 
+    GFT.ProjString(toPROJ4(crsref))
+unsafe_convertcrs(::Type{<:GFT.WellKnownText}, crsref) = 
+    GFT.WellKnownText(GFT.CRS(), toWKT(crsref))
+unsafe_convertcrs(::Type{<:GFT.ESRIWellKnownText}, crsref) =
+    GFT.ESRIWellKnownText(GFT.CRS(), toWKT(morphtoESRI!(crsref)))
+unsafe_convertcrs(::Type{<:GFT.GML}, crsref) = 
+    GFT.GML(toXML(crsref))
+

src/ogr/geometry.jl

@@ -1,3 +1,4 @@
+
 """
     fromWKB(data)
 

src/spatialref.jl

@@ -1,3 +1,99 @@
+"""
+    importCRS(x::GeoFormatTypes.GeoFormat)
+
+Import a coordinate reference system from a `GeoFormat` into GDAL,
+returning an [`AbstractSpatialRef`](@ref).
+"""
+importCRS(x::GFT.GeoFormat) = importCRS!(newspatialref(), x)
+
+unsafe_importCRS(x::GFT.GeoFormat) = importCRS!(unsafe_newspatialref(), x)
+
+"""
+    importCRS!(spref::AbstractSpatialRef, x::GeoFormatTypes.GeoFormat)
+
+Import a coordinate reference system from a `GeoFormat` into the spatial ref.
+"""
+function importCRS! end
+
+importCRS!(spref::AbstractSpatialRef, x::GFT.EPSG) =
+    importEPSG!(spref, GFT.val(x))
+importCRS!(spref::AbstractSpatialRef, x::GFT.AbstractWellKnownText) =
+    importWKT!(spref, GFT.val(x))
+importCRS!(spref::AbstractSpatialRef, x::GFT.ESRIWellKnownText) =
+    importESRI!(spref, GFT.val(x))
+importCRS!(spref::AbstractSpatialRef, x::GFT.ProjString) =
+    importPROJ4!(spref, GFT.val(x))
+importCRS!(spref::AbstractSpatialRef, x::GFT.GML) =
+    importXML!(spref, GFT.val(x))
+importCRS!(spref::AbstractSpatialRef, x::GFT.KML) =
+    importCRS!(spref, GFT.EPSG(4326))
+
+"""
+    reproject(points, sourceproj::GeoFormat, destproj::GeoFormat)
+
+Reproject points to a different coordinate reference system and/or format.
+
+## Arguments
+-`coord`: Vector of Geometry points
+-`sourcecrs`: The current coordinate reference system, as a `GeoFormat`
+-`targetcrs`: The coordinate reference system to transform to, using any CRS capable `GeoFormat`
+
+## Example
+
+```julia-repl
+julia> using ArchGDAL, GeoFormatTypes
+
+julia> ArchGDAL.reproject([[118, 34], [119, 35]], ProjString("+proj=longlat +datum=WGS84 +no_defs"), EPSG(2025))
+2-element Array{Array{Float64,1},1}:
+ [-2.60813482878655e6, 1.5770429674905164e7]
+ [-2.663928675953517e6, 1.56208905951487e7] 
+```
+"""
+reproject(coord, sourcecrs::GFT.GeoFormat, targetcrs::Nothing) = coord
+
+# These should be better integrated with geometry packages or follow some standard
+const ReprojectCoord = Union{<:NTuple{2,<:Number},<:NTuple{3,<:Number},AbstractVector{<:Number}}
+
+# Vector/Tuple coordinate(s)
+reproject(coord::ReprojectCoord, sourcecrs::GFT.GeoFormat, targetcrs::GFT.GeoFormat) =
+    GeoInterface.coordinates(reproject(createpoint(coord...), sourcecrs, targetcrs))
+reproject(coords::AbstractArray{<:ReprojectCoord}, sourcecrs::GFT.GeoFormat, targetcrs::GFT.GeoFormat) =
+    GeoInterface.coordinates.(reproject([createpoint(c...) for c in coords], sourcecrs, targetcrs))
+
+# GeoFormat
+reproject(geom::GFT.GeoFormat, sourcecrs::GFT.GeoFormat, targetcrs::GFT.GeoFormat) =
+    convert(typeof(geom), reproject(convert(Geometry, geom), sourcecrs, targetcrs))
+
+# Geometries
+function reproject(geom::AbstractGeometry, sourcecrs::GFT.GeoFormat, targetcrs::GFT.GeoFormat)
+    crs2transform(sourcecrs, targetcrs) do transform
+        transform!(geom, transform)
+    end
+end
+function reproject(geoms::AbstractArray{<:AbstractGeometry}, sourcecrs::GFT.GeoFormat,
+                   targetcrs::GFT.GeoFormat)
+    crs2transform(sourcecrs, targetcrs) do transform
+        transform!.(geoms, Ref(transform))
+    end
+end
+
+"""
+    crs2transform(f::Function, sourcecrs::GeoFormat, targetcrs::GeoFormat)
+
+Run the function `f` on a coord transform generated from the source and target 
+crs definitions. These can be any `GeoFormat` (from GeoFormatTypes) that holds 
+a coordinate reference system.
+"""
+function crs2transform(f::Function, sourcecrs::GFT.GeoFormat, targetcrs::GFT.GeoFormat)
+    importCRS(sourcecrs) do sourcecrs_ref
+        importCRS(targetcrs) do targetcrs_ref
+            createcoordtrans(sourcecrs_ref, targetcrs_ref) do transform
+                f(transform)
+            end
+        end
+    end
+end
+
 """
     newspatialref(wkt::AbstractString = "")
 

test/runtests.jl

@@ -8,6 +8,7 @@ include("remotefiles.jl")
 @testset "ArchGDAL" begin
     cd(dirname(@__FILE__)) do
         isdir("tmp") || mkpath("tmp")
+        include("test_convert.jl")
         include("test_datastreams.jl")
         include("test_gdal_tutorials.jl")
         include("test_geometry.jl")

test/test_convert.jl

@@ -0,0 +1,30 @@
+using Test
+import ArchGDAL; const AG = ArchGDAL
+import GeoFormatTypes; const GFT = GeoFormatTypes
+
+# Tests high level convert methods
+
+@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 json.val == AG.toJSON(point)
+    @test kml.val == AG.toKML(point)
+    @test gml.val == AG.toGML(point)
+    @test wkb.val == AG.toWKB(point)
+    @test wkt.val == AG.toWKT(point)
+    @test convert(GFT.GeoJSON, json) == convert(GFT.GeoJSON, wkb) == 
+          convert(GFT.GeoJSON, wkt) == convert(GFT.GeoJSON, gml) == json
+    @test convert(GFT.KML, gml) == convert(GFT.KML, wkt)
+end
+
+@testset "convert crs format" begin
+    proj4326 = GFT.ProjString("+proj=longlat +datum=WGS84 +no_defs")
+    @test convert(GFT.ProjString, GFT.CRS(), 
+                  convert(GFT.WellKnownText, GFT.CRS(), 
+                          convert(GFT.ESRIWellKnownText, GFT.CRS(), 
+                                  GFT.EPSG(4326)))) == proj4326
+end

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,23 @@ 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
 end
 
 @testset "Get Projection" begin

test/test_geometry.jl

@@ -1,5 +1,7 @@
 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))
@@ -105,6 +107,7 @@ 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 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)"