Implement a Tables interface for features

Project.toml

@@ -6,19 +6,19 @@ desc = "A high level API for GDAL - Geospatial Data Abstraction Library"
 version = "0.5.0"
 
 [deps]
-DataStreams = "9a8bc11e-79be-5b39-94d7-1ccc349a1a85"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DiskArrays = "3c3547ce-8d99-4f5e-a174-61eb10b00ae3"
 GDAL = "add2ef01-049f-52c4-9ee2-e494f65e021a"
 GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
+Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [compat]
-DataStreams = "0.4.2"
 DiskArrays = "0.2.4"
 GDAL = "1.1.3"
 GeoFormatTypes = "0.3"
 GeoInterface = "0.4, 0.5"
+Tables = "1"
 julia = "1.3"
 
 [extras]

docs/Project.toml

@@ -1,5 +1,6 @@
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DiskArrays = "3c3547ce-8d99-4f5e-a174-61eb10b00ae3"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 

docs/make.jl

@@ -16,6 +16,7 @@ makedocs(
         "GDAL Datasets" => "datasets.md",
         "Feature Data" => "features.md",
         "Raster Data" => "rasters.md",
+        "Tables Interface" => "tables.md",
         "Geometric Operations" => "geometries.md",
         "Spatial Projections" => "projections.md",
         # "Working with Spatialite" => "spatialite.md",

docs/src/datasets.md

@@ -101,3 +101,22 @@ dataset = ArchGDAL.<copy/create/read/update>(...)
 !!! note
 
     This pattern of using `do`-blocks to manage context plays a big way into the way we handle memory in this package. For details, see the section on Memory Management.
+
+The [`ArchGDAL.read`](@ref) method accepts keyword arguments(`kwargs`) viz. the GDAL [open-options](https://gdal.org/drivers/vector/csv.html#open-options) for reading `.csv` spatial datasets.
+
+Example: In a CSV the data is stored as `String`. 
+
+```@example datasets
+dataset1 = ArchGDAL.read("data/multi_geom.csv")
+layer1 = ArchGDAL.getlayer(dataset1, 0)
+```
+
+Well this is weird, the CSV driver recognises our point and linestring geometries as `String`. Now if you have a .csvt file of the same name with the geometry types as `WKT`, they types will be recognized, else, GDAL offers open-options to tweak the read parameters that are passed as `kwargs`.
+
+So for the above CSV, we want the driver to detect our geometries, so according to [open-options](https://gdal.org/drivers/vector/csv.html#open-options) we should use the `"GEOM_POSSIBLE_NAMES=point,linestring"` option. Also we want that the geometry columns should not be kept as regular `String` columns, so we add a `"KEEP_GEOM_COLUMNS=NO"` option too.
+
+```@example datasets
+dataset2 = ArchGDAL.read("data/multi_geom.csv", options = ["GEOM_POSSIBLE_NAMES=point,linestring", "KEEP_GEOM_COLUMNS=NO"])
+
+layer2 = ArchGDAL.getlayer(dataset2, 0)
+```

docs/src/reference.md

@@ -20,6 +20,12 @@ Pages   = ["dataset.jl"]
 Modules = [ArchGDAL]
 Pages   = ["feature.jl", "featuredefn.jl", "featurelayer.jl", "fielddefn.jl", "geometry.jl", "styletable.jl", "context.jl"]
 ```
+## [Tables Interface](@id API-Tables-Interface)
+
+```@autodocs
+Modules = [ArchGDAL]
+Pages   = ["tables.jl"]
+```
 
 ## [Raster Data](@id API-Raster-Data)
 

docs/src/tables.md

@@ -0,0 +1,72 @@
+# Tabular Interface
+
+```@setup tables
+using ArchGDAL
+using DataFrames
+```
+
+ArchGDAL now brings in greater flexibilty in terms of raster data handling via the
+[Tables.jl](https://github.com/JuliaData/Tables.jl) API, that aims to provide a fast and
+responsive tabular interface to data.
+
+In this section, we revisit the
+[`data/point.geojson`](https://github.com/yeesian/ArchGDALDatasets/blob/307f8f0e584a39a050c042849004e6a2bd674f99/data/point.geojson)
+dataset.
+
+```@example tables
+dataset = ArchGDAL.read("data/point.geojson")
+```
+
+Each layer can be represented as a separate Table.
+
+```@example tables
+layer = ArchGDAL.getlayer(dataset, 0)
+```
+
+The [`ArchGDAL.Table`](@ref) method accepts an `ArchGDAL.FeatureLayer`.
+```@example tables
+table = ArchGDAL.Table(layer)
+```
+
+Individual rows can be retrieved using the `Base.getindex(t::ArchGDAL.Table, idx::Int)` method or simply `table[idx]`.
+
+```@example tables
+row = table[1]
+```
+
+Layers are retrievable!
+One can get back the layer that a Table is made up of.
+```@example tables
+lyr = ArchGDAL.getlayer(table)
+```
+
+The Tables interface also support multiple geometries per layer.
+
+Here, we visit 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"])
+
+layer = ArchGDAL.getlayer(dataset, 0)
+table = ArchGDAL.Table(layer)
+```
+
+Exatracting a row from the table, we see that the row/feature is made up of two geometries
+viz. `point` and `linestring`.
+```@example tables
+row = table[1]
+```
+
+Finally layers can be converted to DataFrames to perform miscellaneous spatial operations.
+```@example tables
+df = DataFrame(table)
+```
+In some cases the `nextfeature` might become a bit tedious to use. In which case the `ArchGDAL.nextnamedtuple()` method comes in handy. Though built upon `nextfeature`, simply calling it, yields the `feature` as a `NamedTuple`. Though one might have to use `ArchGDAL.resetreading!(layer)` method to reset the layer reading to the start.
+
+```@example tables
+ArchGDAL.resetreading!(layer)
+feat1 = ArchGDAL.nextnamedtuple(layer)
+feat2 = ArchGDAL.nextnamedtuple(layer)
+```

src/ArchGDAL.jl

@@ -1,11 +1,10 @@
 module ArchGDAL
 
-    import GDAL, GeoInterface, GeoFormatTypes
-    import DataStreams: Data
-    import GeoInterface: coordinates, geotype
-    import Base: convert
-
     using Dates
+    using GDAL: GDAL
+    using GeoFormatTypes: GeoFormatTypes
+    using GeoInterface: GeoInterface
+    using Tables: Tables
 
     const GFT = GeoFormatTypes
 
@@ -30,7 +29,7 @@ module ArchGDAL
     include("context.jl")
     include("base/iterators.jl")
     include("base/display.jl")
-    include("datastreams.jl")
+    include("tables.jl")
     include("geointerface.jl")
     include("convert.jl")
 

src/dataset.jl

@@ -460,8 +460,10 @@ unsafe_getlayer(dataset::AbstractDataset, i::Integer) =
 
 """
     getlayer(dataset::AbstractDataset, name::AbstractString)
+    getlayer(table::Table)
 
-Fetch the feature layer corresponding to the given name.
+Fetch the feature layer corresponding to the given name. If it is called on a Table, which
+supports only one layer, a name is not needed.
 
 The returned layer remains owned by the GDALDataset and should not be deleted by
 the application.

src/ogr/feature.jl

@@ -95,7 +95,7 @@ getfielddefn(feature::Feature, i::Integer) =
     IFieldDefnView(GDAL.ogr_f_getfielddefnref(feature.ptr, i))
 
 """
-    findfieldindex(feature::Feature, name::AbstractString)
+    findfieldindex(feature::Feature, name::Union{AbstractString, Symbol})
 
 Fetch the field index given field name.
 
@@ -109,7 +109,7 @@ the field index, or -1 if no matching field is found.
 ### Remarks
 This is a cover for the `OGRFeatureDefn::GetFieldIndex()` method.
 """
-findfieldindex(feature::Feature, name::AbstractString) =
+findfieldindex(feature::Feature, name::Union{AbstractString, Symbol}) =
     GDAL.ogr_f_getfieldindex(feature.ptr, name)
 
 """
@@ -381,7 +381,7 @@ function getfield(feature::Feature, i::Integer)
     end
 end
 
-getfield(feature::Feature, name::AbstractString) =
+getfield(feature::Feature, name::Union{AbstractString, Symbol}) =
     getfield(feature, findfieldindex(feature, name))
 
 """

src/ogr/featuredefn.jl

@@ -91,7 +91,7 @@ getfielddefn(featuredefn::IFeatureDefnView, i::Integer) =
     IFieldDefnView(GDAL.ogr_fd_getfielddefn(featuredefn.ptr, i))
 
 """
-    findfieldindex(featuredefn::AbstractFeatureDefn, name::AbstractString)
+    findfieldindex(featuredefn::AbstractFeatureDefn, name::Union{AbstractString, Symbol})
 
 Find field by name.
 
@@ -101,7 +101,7 @@ the field index, or -1 if no match found.
 ### Remarks
 This uses the OGRFeatureDefn::GetFieldIndex() method.
 """
-findfieldindex(featuredefn::AbstractFeatureDefn, name::AbstractString) =
+findfieldindex(featuredefn::AbstractFeatureDefn, name::Union{AbstractString, Symbol}) =
     GDAL.ogr_fd_getfieldindex(featuredefn.ptr, name)
 
 """

src/ogr/featurelayer.jl

@@ -524,7 +524,7 @@ layerdefn(layer::AbstractFeatureLayer) =
     IFeatureDefnView(GDAL.ogr_l_getlayerdefn(layer.ptr))
 
 """
-    findfieldindex(layer::AbstractFeatureLayer, field::AbstractString, exactmatch::Bool)
+    findfieldindex(layer::AbstractFeatureLayer, field::Union{AbstractString, Symbol}, exactmatch::Bool)
 
 Find the index of the field in a layer, or -1 if the field doesn't exist.
 
@@ -534,7 +534,7 @@ the layer was created (eg. like `LAUNDER` in the OCI driver).
 """
 function findfieldindex(
         layer::AbstractFeatureLayer,
-        field::AbstractString,
+        field::Union{AbstractString, Symbol},
         exactmatch::Bool
     )
     return GDAL.ogr_l_findfieldindex(layer.ptr, field, exactmatch)

src/tables.jl

@@ -0,0 +1,73 @@
+"""
+Constructs `Table` out of `FeatureLayer`, where every row is a `Feature` consisting of Geometry and attributes.
+```
+ArchGDAL.Table(T::Union{IFeatureLayer, FeatureLayer})
+```
+"""
+struct Table{T<:Union{IFeatureLayer, FeatureLayer}}
+    layer::T
+end
+
+getlayer(t::Table) = Base.getfield(t, :layer)
+
+function Tables.schema(layer::AbstractFeatureLayer)
+    field_names, geom_names, featuredefn, fielddefns = schema_names(layer)
+    ngeom = ArchGDAL.ngeom(featuredefn)
+    geomdefns = (ArchGDAL.getgeomdefn(featuredefn, i) for i in 0:ngeom-1)
+    field_types = (_FIELDTYPE[gettype(fielddefn)] for fielddefn in fielddefns)
+    geom_types = (IGeometry for i in 1:ngeom)
+    Tables.Schema((field_names..., geom_names...), (field_types..., geom_types...))
+end
+
+Tables.istable(::Type{<:Table}) = true
+Tables.rowaccess(::Type{<:Table}) = true
+Tables.rows(t::Table) = t
+
+function Base.iterate(t::Table, st = 0)
+    layer = getlayer(t)
+    st >= nfeature(layer) && return nothing
+    if iszero(st)
+        resetreading!(layer)
+    end
+    return nextnamedtuple(layer), st + 1
+end
+
+function Base.getindex(t::Table, idx::Int)
+    layer = getlayer(t)
+    setnextbyindex!(layer, idx-1) 
+    return nextnamedtuple(layer)
+end
+
+Base.IteratorSize(::Type{<:Table}) = Base.HasLength()
+Base.size(t::Table) = nfeature(getlayer(t))
+Base.length(t::Table) = size(t)
+Base.IteratorEltype(::Type{<:Table}) = Base.HasEltype()
+Base.propertynames(t::Table) = Tables.schema(getlayer(t)).names
+Base.getproperty(t::Table, s::Symbol) = [getproperty(row, s) for row in t]
+
+"""
+Returns the feature row of a layer as a `NamedTuple`
+
+Calling it iteratively will work similar to `nextfeature` i.e. give the consecutive feature as `NamedTuple`
+"""
+function nextnamedtuple(layer::IFeatureLayer)
+    field_names, geom_names = schema_names(layer)
+    return nextfeature(layer) do feature
+        prop = (getfield(feature, name) for name in field_names)
+        geom = (getgeom(feature, idx-1) for idx in 1:length(geom_names))
+        NamedTuple{(field_names..., geom_names...)}((prop..., geom...))
+    end
+end
+
+function schema_names(layer::AbstractFeatureLayer)
+    featuredefn = layerdefn(layer)
+    fielddefns = (getfielddefn(featuredefn, i) for i in 0:nfield(layer)-1)
+    field_names = (Symbol(getname(fielddefn)) for fielddefn in fielddefns)
+    geom_names = (Symbol(getname(getgeomdefn(featuredefn, i-1))) for i in 1:ngeom(layer))
+    return (field_names, geom_names, featuredefn, fielddefns)
+end
+
+function Base.show(io::IO, t::Table)
+    println(io, "Table with $(nfeature(getlayer(t))) features")
+end
+Base.show(io::IO, ::MIME"text/plain", t::Table) = show(io, t)

test/remotefiles.jl

@@ -10,7 +10,18 @@ const testdatadir = @__DIR__
 REPO_URL = "https://github.com/yeesian/ArchGDALDatasets/blob/master/"
 
 # remote files with SHA-2 256 hash
+"""
+To add more files, follow the below steps to generate the SHA
+```
+julia> using SHA
+julia> open(filepath/filename) do f
+           bytes2hex(sha256(f))
+       end
+```
+"""
 remotefiles = [
+    ("data/multi_geom.csv", "00520017658b66ff21e40cbf553672fa8e280cddae6e7a5d1f8bd36bcd521770"),
+    ("data/missing_testcase.csv", "d49ba446aae9ef334350b64c876b4de652f28595fdecf78bea4e16af4033f7c6"),
     ("data/point.geojson", "8744593479054a67c784322e0c198bfa880c9388b39a2ddd4c56726944711bd9"),
     ("data/utmsmall.tif", "f40dae6e8b5e18f3648e9f095e22a0d7027014bb463418d32f732c3756d8c54f"),
     ("gdalworkshop/world.tif", "b376dc8af62f9894b5050a6a9273ac0763ae2990b556910d35d4a8f4753278bb"),