Table to IFeatureLayer

docs/src/tables.md

@@ -1,28 +1,70 @@
 # 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 must 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 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
+```
+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
+- 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/tables.jl

@@ -47,3 +47,218 @@ function schema_names(featuredefn::IFeatureDefnView)
     )
     return (geom_names, field_names, featuredefn, fielddefns)
 end
+
+"""
+    _convert_cleantype_to_AGtype(T)
+
+Converts type `T` into either:
+- a `OGRwkbGeometryType` or
+- a tuple of `OGRFieldType` and `OGRFieldSubType`
+
+"""
+function _convert_cleantype_to_AGtype end
+_convert_cleantype_to_AGtype(::Type{IGeometry}) = wkbUnknown
+@generated _convert_cleantype_to_AGtype(::Type{IGeometry{U}}) where U = :($U)
+@generated _convert_cleantype_to_AGtype(T::Type{U}) where U = :(convert(OGRFieldType, T), convert(OGRFieldSubType, T))
+
+
+"""
+    _convert_coltype_to_cleantype(T)
+
+Convert a table column type to a "clean" type:
+- Unions are flattened
+- Missing and Nothing are dropped
+- Resulting mixed types are approximated by their tightest common supertype
+
+"""
+function _convert_coltype_to_cleantype(T::Type)
+    flattened_T = Base.uniontypes(T)
+    clean_flattened_T = filter(t -> t ∉ [Missing, Nothing], flattened_T)
+    return promote_type(clean_flattened_T...)
+end
+
+"""
+    _fromtable(sch, rows; name)
+
+Converts a row table `rows` with schema `sch` to a layer (optionally named `name`) within a MEMORY dataset
+
+"""
+function _fromtable end
+
+"""
+    _fromtable(sch::Tables.Schema{names,types}, rows; name::String = "")
+
+Handles the case where names and types in `sch` are different from `nothing`
+
+# Implementation
+1. convert `rows`'s column types given in `sch` to either geometry types or field types and subtypes
+2. split `rows`'s columns into geometry typed columns and field typed columns
+3. create layer named `name` in a MEMORY dataset geomfields and fields types inferred from `rows`'s column types
+4. populate layer with `rows` values
+
+"""
+function _fromtable(
+    sch::Tables.Schema{names,types},
+    rows;
+    name::String = "",
+)::IFeatureLayer where {names,types}
+    # TODO maybe constrain `names`
+    strnames = string.(sch.names)
+
+    # Convert column types to either geometry types or field types and subtypes
+    AG_types = Vector{Union{OGRwkbGeometryType,Tuple{OGRFieldType,OGRFieldSubType}}}(undef, length(Tables.columnnames(rows)))
+    for (i, (coltype, colname)) in enumerate(zip(sch.types, strnames))
+        AG_types[i] = try
+            (_convert_cleantype_to_AGtype ∘ _convert_coltype_to_cleantype)(coltype)
+        catch e
+            if e isa MethodError
+                error("Cannot convert column \"$colname\" (type $coltype) to neither IGeometry{::OGRwkbGeometryType} or OGRFieldType and OGRFieldSubType",)
+            else
+                rethrow()
+            end
+        end
+    end
+
+    # Split names and types: between geometry type columns and field type columns
+    geomindices = isa.(AG_types, OGRwkbGeometryType)
+    !any(geomindices) && error("No column convertible to geometry")
+    geomtypes = AG_types[geomindices] # TODO consider to use a view
+    geomnames = strnames[geomindices]
+
+    fieldindices = isa.(AG_types, Tuple{OGRFieldType,OGRFieldSubType})
+    fieldtypes = AG_types[fieldindices] # TODO consider to use a view
+    fieldnames = strnames[fieldindices]
+
+    # Create layer
+    layer = createlayer(name = name, geom = first(geomtypes))
+    # TODO: create setname! for IGeomFieldDefnView. Probably needs first to fix issue #215
+    # TODO: "Model and handle relationships between GDAL objects systematically"
+    GDAL.ogr_gfld_setname(
+        getgeomdefn(layerdefn(layer), 0).ptr,
+        first(geomnames),
+    )
+
+    # Create FeatureDefn
+    if length(geomtypes) ≥ 2
+        for (j, geomtype) in enumerate(geomtypes[2:end])
+            creategeomdefn(geomnames[j+1], geomtype) do geomfielddefn
+                return addgeomdefn!(layer, geomfielddefn) # TODO check if necessary/interesting to set approx=true
+            end
+        end
+    end
+    for (j, (ft, fst)) in enumerate(fieldtypes)
+        createfielddefn(fieldnames[j], ft) do fielddefn
+            setsubtype!(fielddefn, fst)
+            return addfielddefn!(layer, fielddefn)
+        end
+    end
+
+    # Populate layer
+    for (i, row) in enumerate(rows)
+        rowvalues =
+            [Tables.getcolumn(row, col) for col in Tables.columnnames(row)]
+        rowgeoms = view(rowvalues, geomindices)
+        rowfields = view(rowvalues, fieldindices)
+        addfeature(layer) do feature
+            # For geometry fields both `missing` and `nothing` map to not geometry set
+            # since in GDAL <= v"3.3.2", special fields as geometry field cannot be NULL
+            # cf. `OGRFeature::IsFieldNull( int iField )` implemetation
+            for (j, val) in enumerate(rowgeoms)
+                val !== missing &&
+                    val !== nothing &&
+                    setgeom!(feature, j - 1, val)
+            end
+            for (j, val) in enumerate(rowfields)
+                if val === missing
+                    setfieldnull!(feature, j - 1)
+                elseif val !== nothing
+                    setfield!(feature, j - 1, val)
+                end
+            end
+        end
+    end
+
+    return layer
+end
+
+"""
+    _fromtable(::Tables.Schema{names,nothing}, rows; name::String = "")
+
+Handles the case where types in schema is `nothing`
+
+# Implementation
+Tables.Schema types are extracted from `rows`'s columns element types before calling `_fromtable(Tables.Schema(names, types), rows; name = name)`
+
+"""
+function _fromtable(
+    ::Tables.Schema{names,nothing},
+    rows;
+    name::String = "",
+)::IFeatureLayer where {names}
+    cols = Tables.columns(rows)
+    types = (eltype(collect(col)) for col in cols)
+    return _fromtable(Tables.Schema(names, types), rows; name = name)
+end
+
+"""
+    _fromtable(::Tables.Schema{names,nothing}, rows; name::String = "")
+
+Handles the case where schema is `nothing`
+
+# Implementation
+Tables.Schema names are extracted from `rows`'s columns names before calling `_fromtable(Tables.Schema(names, types), rows; name = name)`
+
+"""
+function _fromtable(::Nothing, rows; name::String = "")::IFeatureLayer
+    state = iterate(rows)
+    state === nothing && return IFeatureLayer()
+    row, _ = state
+    names = Tables.columnnames(row)
+    return _fromtable(Tables.Schema(names, nothing), rows; name = name)
+end
+
+"""
+    IFeatureLayer(table; name="")
+
+Construct an IFeatureLayer from a source implementing Tables.jl interface
+
+## Restrictions
+- Source must 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)
+
+## Returns
+An IFeatureLayer within a **MEMORY** driver dataset
+
+## Examples
+```jldoctest
+julia> using ArchGDAL; AG = ArchGDAL
+ArchGDAL
+
+julia> nt = NamedTuple([
+           :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"],
+       ])
+(point = Union{Missing, ArchGDAL.IGeometry{ArchGDAL.wkbPoint}}[Geometry: POINT (30 10), missing], mixedgeom = ArchGDAL.IGeometry[Geometry: POINT (5 10), Geometry: LINESTRING (30 10,10 30)], id = ["5.1", "5.2"], zoom = [1.0, 2.0], location = Union{Missing, String}[missing, "New Delhi"])
+
+julia> layer = AG.IFeatureLayer(nt; name="towns")
+Layer: towns
+  Geometry 0 (point): [wkbPoint]
+  Geometry 1 (mixedgeom): [wkbUnknown]
+     Field 0 (id): [OFTString], 5.1, 5.2
+     Field 1 (zoom): [OFTReal], 1.0, 2.0
+     Field 2 (location): [OFTString], missing, New Delhi
+```
+"""
+function IFeatureLayer(table; name::String = "")::IFeatureLayer
+    # Check tables interface's conformance
+    !Tables.istable(table) &&
+        throw(DomainError(table, "$table has not a Table interface"))
+    # Extract table data
+    rows = Tables.rows(table)
+    schema = Tables.schema(table)
+    return _fromtable(schema, rows; name = name)
+end