Table to IFeatureLayer

benchmark/remotefiles.jl

@@ -16,12 +16,10 @@ julia> open(filepath/filename) do f
        end
 ```
 """
-remotefiles = [
-    (
-        "data/road.zip",
-        "058bdc549d0fc5bfb6deaef138e48758ca79ae20df79c2fb4c40cb878f48bfd8",
-    ),
-]
+remotefiles = [(
+    "data/road.zip",
+    "058bdc549d0fc5bfb6deaef138e48758ca79ae20df79c2fb4c40cb878f48bfd8",
+)]
 
 function verify(path::AbstractString, hash::AbstractString)
     @assert occursin(r"^[0-9a-f]{64}$", hash)

docs/src/tables.md

@@ -1,28 +1,71 @@
 # Tabular Interface
 
 ```@setup tables
-using ArchGDAL, DataFrames
+using ArchGDAL; AG = ArchGDAL
+using DataFrames
 ```
 
 ArchGDAL now brings in greater flexibilty in terms of vector data handling via the
 [Tables.jl](https://github.com/JuliaData/Tables.jl) API. In general, tables are modelled based on feature layers and support multiple geometries per layer. Namely, the layer(s) of a dataset can be converted to DataFrame(s) to perform miscellaneous spatial operations.
 
+## Conversion to table
+
 Here is a quick example based on the
 [`data/point.geojson`](https://github.com/yeesian/ArchGDALDatasets/blob/307f8f0e584a39a050c042849004e6a2bd674f99/data/point.geojson)
 dataset:
 
-```@example tables
-dataset = ArchGDAL.read("data/point.geojson")
-
-DataFrames.DataFrame(ArchGDAL.getlayer(dataset, 0))
+```@repl tables
+ds = AG.read("data/point.geojson")
+DataFrame(AG.getlayer(ds, 0))
 ```
 
 To illustrate multiple geometries, here is a second example based on the
 [`data/multi_geom.csv`](https://github.com/yeesian/ArchGDALDatasets/blob/master/data/multi_geom.csv)
 dataset:
 
-```@example tables
-dataset1 = ArchGDAL.read("data/multi_geom.csv", options = ["GEOM_POSSIBLE_NAMES=point,linestring", "KEEP_GEOM_COLUMNS=NO"])
+```@repl tables
+ds = AG.read("data/multi_geom.csv", options = ["GEOM_POSSIBLE_NAMES=point,linestring", "KEEP_GEOM_COLUMNS=NO"])
+DataFrame(AG.getlayer(ds, 0))
+```
+## Conversion to layer
+A table-like source implementing Tables.jl interface can be converted to a layer, provided that:
+- Source contains at least one geometry column
+- Geometry columns are recognized by their element type being a subtype of `Union{IGeometry, Nothing,  Missing}`
+- Non geometry columns must contain types handled by GDAL/OGR (e.g. not `Int128` nor composite type)
+
+**Note**: As geometries and fields are stored separately in GDAL features, the backward conversion of the layer won't have the same column ordering. Geometry columns will be the first columns.
 
-DataFrames.DataFrame(ArchGDAL.getlayer(dataset1, 0))
+```@repl tables
+df = DataFrame([
+    :point => [AG.createpoint(30, 10), missing],
+    :mixedgeom => [AG.createpoint(5, 10), AG.createlinestring([(30.0, 10.0), (10.0, 30.0)])],
+    :id => ["5.1", "5.2"],
+    :zoom => [1.0, 2],
+    :location => [missing, "New Delhi"],
+])
+layer = AG.IFeatureLayer(df)
+```
+
+The layer, converted from a source implementing the Tables.jl interface, will be in a [memory](https://gdal.org/drivers/vector/memory.html) dataset.  
+Hence you can:
+- Add other layers to it
+- Copy it to a dataset with another driver
+- Write it to a file
+### Example of writing with ESRI Shapefile driver
+```@repl tables
+ds = AG.write(layer.ownedby, "test.shp", driver=AG.getdriver("ESRI Shapefile"))
+DataFrame(AG.getlayer(AG.read("test.shp"), 0))
+rm.(["test.shp", "test.shx", "test.dbf"]) # hide
+```
+!!! warning
+    As OGR ESRI Shapefile driver [does not support multi geometries](https://gdal.org/development/rfc/rfc41_multiple_geometry_fields.html#drivers), the second geometry has been dropped
+!!! warning
+    As OGR ESRI Shapefile driver does not support nullable fields, the `missing` location has been replaced by `""`
+### Example of writing with GML driver
+Using the GML 3.2.1 more capable driver/format, you can write more information to the file
+```@repl tables
+ds = AG.write(layer.ownedby, "test.gml", driver=AG.getdriver("GML"), options=["FORMAT=GML3.2"])
+DataFrame(AG.getlayer(AG.read("test.gml", options=["EXPOSE_GML_ID=NO"]), 0))
+rm.(["test.gml", "test.xsd"]) # hide
 ```
+**Note:** [OGR GML driver](https://gdal.org/drivers/vector/gml.html#open-options) option **EXPOSE\_GML\_ID=NO** avoids to read the `gml_id` field, mandatory in GML 3.x format and automatically created by the OGR GML driver

src/base/display.jl

@@ -128,7 +128,8 @@ function Base.show(io::IO, layer::AbstractFeatureLayer)::Nothing
         end
         if ngeomdisplay == 1 # only support printing of a single geom column
             for f in layer
-                geomwkt = toWKT(getgeom(f))
+                geom = getgeom(f)
+                geomwkt = geom.ptr != C_NULL ? toWKT(geom) : "NULL Geometry"
                 length(geomwkt) > 25 && (geomwkt = "$(geomwkt[1:20])...)")
                 newdisplay = "$display, $geomwkt"
                 if length(newdisplay) > 75
@@ -251,26 +252,28 @@ function Base.show(io::IO, spref::AbstractSpatialRef)::Nothing
     return nothing
 end
 
-function Base.show(io::IO, geom::AbstractGeometry)::Nothing
+function Base.show(io::IO, geom::AbstractGeometry, prefix::String)::Nothing
     if geom.ptr == C_NULL
-        print(io, "NULL Geometry")
+        print(io, "NULL $(prefix)Geometry")
         return nothing
     end
     compact = get(io, :compact, false)
 
     if !compact
-        print(io, "Geometry: ")
+        print(io, "$(prefix)Geometry: ")
         geomwkt = toWKT(geom)
         if length(geomwkt) > 60
             print(io, "$(geomwkt[1:50]) ... $(geomwkt[(end - 4):end])")
         else
             print(io, "$geomwkt")
         end
     else
-        print(io, "Geometry: $(getgeomtype(geom))")
+        print(io, "$(prefix)Geometry: $(getgeomtype(geom))")
     end
     return nothing
 end
+Base.show(io::IO, geom::IGeometry) = show(io, geom, "I")
+Base.show(io::IO, geom::Geometry) = show(io, geom, "")
 
 function Base.show(io::IO, ct::ColorTable)::Nothing
     if ct.ptr == C_NULL

src/ogr/geometry.jl

@@ -34,6 +34,8 @@ function unsafe_fromWKB(data)::Geometry
     return Geometry(geom[])
 end
 
+convert(::Type{IGeometry}, data::Vector{UInt8}) = fromWKB(data)
+
 """
     fromWKT(data::Vector{String})
 
@@ -74,6 +76,8 @@ fromWKT(data::String, args...)::IGeometry = fromWKT([data], args...)
 unsafe_fromWKT(data::String, args...)::Geometry =
     unsafe_fromWKT([data], args...)
 
+convert(::Type{IGeometry}, s::String) = fromWKT(s)
+
 """
 Destroy geometry object.
 
@@ -1636,3 +1640,32 @@ for (f, rt) in ((:create, :IGeometry), (:unsafe_create, :Geometry))
         end
     end
 end
+
+# Conversion from GeoInterface geometry
+function convert(::Type{IGeometry}, g::GeoInterface.AbstractPoint)
+    return createpoint(GeoInterface.coordinates(g))
+end
+function convert(::Type{IGeometry}, g::GeoInterface.AbstractMultiPoint)
+    return createmultipoint(GeoInterface.coordinates(g))
+end
+function convert(::Type{IGeometry}, g::GeoInterface.AbstractLineString)
+    return createlinestring(GeoInterface.coordinates(g))
+end
+function convert(::Type{IGeometry}, g::GeoInterface.AbstractMultiLineString)
+    return createmultilinestring(GeoInterface.coordinates(g))
+end
+function convert(::Type{IGeometry}, g::GeoInterface.AbstractPolygon)
+    return createpolygon(GeoInterface.coordinates(g))
+end
+function convert(::Type{IGeometry}, g::GeoInterface.AbstractMultiPolygon)
+    return createmultipolygon(GeoInterface.coordinates(g))
+end
+function convert(::Type{IGeometry}, g::GeoInterface.AbstractGeometryCollection)
+    ag_geom = creategeom(wkbGeometryCollection)
+    for gi_subgeom in GeoInterface.geometries(g)
+        ag_subgeom = convert(IGeometry, gi_subgeom)
+        result = GDAL.ogr_g_addgeometry(ag_geom.ptr, ag_subgeom.ptr)
+        @ogrerr result "Failed to add $ag_subgeom"
+    end
+    return ag_geom
+end