diff --git a/benchmark/remotefiles.jl b/benchmark/remotefiles.jl index a5af3477..d0e7d40b 100644 --- a/benchmark/remotefiles.jl +++ b/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) diff --git a/docs/src/tables.md b/docs/src/tables.md index bd80000b..cd67934a 100644 --- a/docs/src/tables.md +++ b/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 diff --git a/src/base/display.jl b/src/base/display.jl index 3f7cbe11..ea0c1f9e 100644 --- a/src/base/display.jl +++ b/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,15 +252,15 @@ 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])") @@ -267,10 +268,12 @@ function Base.show(io::IO, geom::AbstractGeometry)::Nothing 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 diff --git a/src/ogr/geometry.jl b/src/ogr/geometry.jl index 6c0c65c9..bb12473e 100644 --- a/src/ogr/geometry.jl +++ b/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 diff --git a/src/tables.jl b/src/tables.jl index 7d42f0fa..0a5c6c17 100644 --- a/src/tables.jl +++ b/src/tables.jl @@ -47,3 +47,266 @@ 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 +@generated function _convert_cleantype_to_AGtype( + T::Type{U}, +) where {U<:GeoInterface.AbstractGeometry} + return :(convert(OGRwkbGeometryType, T)) +end +@generated function _convert_cleantype_to_AGtype(T::Type{U}) where {U} + return :(convert(OGRFieldType, T), convert(OGRFieldSubType, T)) +end + +""" + _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 + +function _create_empty_layer_from_AGtypes( + colnames::NTuple{N,String}, + AGtypes::Vector{ + Union{OGRwkbGeometryType,Tuple{OGRFieldType,OGRFieldSubType}}, + }, + name::String, +) where {N} + # Split names and types: between geometry type columns and field type columns + geomindices = isa.(AGtypes, OGRwkbGeometryType) + !any(geomindices) && error("No column convertible to geometry") + geomtypes = AGtypes[geomindices] # TODO consider to use a view + geomnames = colnames[geomindices] + + fieldindices = isa.(AGtypes, Tuple{OGRFieldType,OGRFieldSubType}) + fieldtypes = AGtypes[fieldindices] # TODO consider to use a view + fieldnames = colnames[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 (geomtype, geomname) in zip(geomtypes[2:end], geomnames[2:end]) + creategeomdefn(geomname, geomtype) do geomfielddefn + return addgeomdefn!(layer, geomfielddefn) # TODO check if necessary/interesting to set approx=true + end + end + end + for (fieldname, (fieldtype, fieldsubtype)) in zip(fieldnames, fieldtypes) + createfielddefn(fieldname, fieldtype) do fielddefn + setsubtype!(fielddefn, fieldsubtype) + return addfielddefn!(layer, fielddefn) + end + end + + return layer, geomindices, fieldindices +end + +""" + _infergeometryorfieldtypes(sch, rows) + +Infer ArchGDAL field and geometry types from schema + +""" +function _infergeometryorfieldtypes( + sch::Tables.Schema{names,types}, + rows, +) where {names,types} + colnames = string.(sch.names) + + # Convert column types to either geometry types or field types and subtypes + AGtypes = + Vector{Union{OGRwkbGeometryType,Tuple{OGRFieldType,OGRFieldSubType}}}( + undef, + length(Tables.columnnames(rows)), + ) + for (j, (coltype, colname)) in enumerate(zip(sch.types, colnames)) + # we wrap the following in a try-catch block to surface the original column type (rather than clean/converted type) in the error message + AGtypes[j] = 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 + throw(e) + end + end + end + + return AGtypes +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; + layer_name::String, +)::IFeatureLayer where {names,types} + AGtypes = _infergeometryorfieldtypes(sch, rows) + + (layer, geomindices, fieldindices) = _create_empty_layer_from_AGtypes( + string.(sch.names), + AGtypes, + layer_name, + ) + + # Populate layer + for row in 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) + if val !== missing && val !== nothing + setgeom!(feature, j - 1, convert(IGeometry, val)) + end + 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; + kwargs..., +)::IFeatureLayer where {names} + cols = Tables.columns(rows) + types = (eltype(collect(col)) for col in cols) + return _fromtable(Tables.Schema(names, types), rows; kwargs...) +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; kwargs...)::IFeatureLayer + state = iterate(rows) + state === nothing && return IFeatureLayer() + row, _ = state + names = Tables.columnnames(row) + return _fromtable(Tables.Schema(names, nothing), rows; kwargs...) +end + +""" + IFeatureLayer(table; kwargs...) + +Construct an IFeatureLayer from a source implementing Tables.jl interface + +## Keyword arguments +- `layer_name::String = ""`: name of the layer + +## 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}` or + - `Union{GeoInterface.AbstractGeometry, Nothing, Missing}` or + - `Union{String, Nothing, Missing}` provided that String values can be decoded as WKT or + - `Union{Vector{UInt8}, Nothing, Missing}` provided that Vector{UInt8} values can be decoded as WKB +- 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 + +## Example +```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}}[IGeometry: POINT (30 10), missing], mixedgeom = ArchGDAL.IGeometry[IGeometry: POINT (5 10), IGeometry: 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; layer_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; layer_name::String = "layer") + # 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; layer_name = layer_name) +end diff --git a/src/types.jl b/src/types.jl index a2936a2b..e994dd88 100644 --- a/src/types.jl +++ b/src/types.jl @@ -291,15 +291,45 @@ end OFTInteger64List::GDAL.OFTInteger64List, ) +const OGRFieldcompatibleDataTypes = Dict( + Bool => (OFTInteger, OFSTBoolean), + Int8 => (OFTInteger, OFSTNone), + Int16 => (OFTInteger, OFSTInt16), + Int32 => (OFTInteger, OFSTNone), + Vector{Bool} => (OFTIntegerList, OFSTBoolean), + Vector{Int16} => (OFTIntegerList, OFSTInt16), + Vector{Int32} => (OFTIntegerList, OFSTNone), + Float16 => (OFTReal, OFSTNone), + Float32 => (OFTReal, OFSTFloat32), + Float64 => (OFTReal, OFSTNone), + Vector{Float16} => (OFTRealList, OFSTNone), + Vector{Float32} => (OFTRealList, OFSTFloat32), + Vector{Float64} => (OFTRealList, OFSTNone), + String => (OFTString, OFSTNone), + Vector{String} => (OFTStringList, OFSTNone), + Vector{UInt8} => (OFTBinary, OFSTNone), + Dates.Date => (OFTDate, OFSTNone), + Dates.Time => (OFTTime, OFSTNone), + Dates.DateTime => (OFTDateTime, OFSTNone), + Int64 => (OFTInteger64, OFSTNone), + Vector{Int64} => (OFTInteger64List, OFSTNone), +) + @convert( OGRFieldType::DataType, OFTInteger::Bool, + OFTInteger::Int8, OFTInteger::Int16, OFTInteger::Int32, # default type comes last - OFTIntegerList::Vector{Int32}, + OFTIntegerList::Vector{Bool}, + OFTIntegerList::Vector{Int16}, + OFTIntegerList::Vector{Int32}, # default type comes last + OFTReal::Float16, OFTReal::Float32, OFTReal::Float64, # default type comes last - OFTRealList::Vector{Float64}, + OFTRealList::Vector{Float16}, + OFTRealList::Vector{Float32}, + OFTRealList::Vector{Float64}, # default type comes last OFTString::String, OFTStringList::Vector{String}, OFTBinary::Vector{UInt8}, @@ -321,11 +351,29 @@ end @convert( OGRFieldSubType::DataType, - OFSTNone::Nothing, - OFSTBoolean::Bool, - OFSTInt16::Int16, - OFSTFloat32::Float32, - OFSTJSON::String, + OFSTNone::Int8, + OFSTNone::Int32, + OFSTBoolean::Vector{Bool}, + OFSTBoolean::Bool, # default type comes last + OFSTInt16::Vector{Int16}, + OFSTInt16::Int16, # default type comes last + OFSTNone::Vector{Int32}, + OFSTInt16::Float16, + OFSTNone::Float64, + OFSTInt16::Vector{Float16}, + OFSTFloat32::Vector{Float32}, + OFSTFloat32::Float32, # default type comes last + OFSTNone::Vector{Float64}, + OFSTNone::String, + OFSTNone::Vector{String}, + OFSTNone::Vector{UInt8}, + OFSTNone::Dates.Date, + OFSTNone::Dates.Time, + OFSTNone::Dates.DateTime, + OFSTNone::Int64, + OFSTNone::Vector{Int64}, + # Lacking OFSTUUID and OFSTJSON defined in GDAL ≥ v"3.3" + OFSTNone::Nothing, # default type comes last ) @convert( @@ -510,6 +558,96 @@ end wkbGeometryCollection25D::GDAL.wkbGeometryCollection25D, ) +@generated function convert( + T1::Type{OGRwkbGeometryType}, + T2::Type{U}, +) where {U<:GeoInterface.AbstractGeometry} + U <: GeoInterface.AbstractPoint && return :(wkbPoint) + U <: GeoInterface.AbstractMultiPoint && return :(wkbMultiPoint) + U <: GeoInterface.AbstractLineString && return :(wkbLineString) + U <: GeoInterface.AbstractMultiLineString && return :(wkbMultiLineString) + U <: GeoInterface.AbstractPolygon && return :(wkbPolygon) + U <: GeoInterface.AbstractMultiPolygon && return :(wkbMultiPolygon) + U == GeoInterface.AbstractGeometry && return :(wkbUnknown) + return :(error("No convert method to convert $T2 to $T1")) +end + +@convert( + OGRwkbGeometryType::IGeometry, + wkbUnknown::IGeometry{wkbUnknown}, + wkbPoint::IGeometry{wkbPoint}, + wkbLineString::IGeometry{wkbLineString}, + wkbPolygon::IGeometry{wkbPolygon}, + wkbMultiPoint::IGeometry{wkbMultiPoint}, + wkbMultiLineString::IGeometry{wkbMultiLineString}, + wkbMultiPolygon::IGeometry{wkbMultiPolygon}, + wkbGeometryCollection::IGeometry{wkbGeometryCollection}, + wkbCircularString::IGeometry{wkbCircularString}, + wkbCompoundCurve::IGeometry{wkbCompoundCurve}, + wkbCurvePolygon::IGeometry{wkbCurvePolygon}, + wkbMultiCurve::IGeometry{wkbMultiCurve}, + wkbMultiSurface::IGeometry{wkbMultiSurface}, + wkbCurve::IGeometry{wkbCurve}, + wkbSurface::IGeometry{wkbSurface}, + wkbPolyhedralSurface::IGeometry{wkbPolyhedralSurface}, + wkbTIN::IGeometry{wkbTIN}, + wkbTriangle::IGeometry{wkbTriangle}, + wkbNone::IGeometry{wkbNone}, + wkbLinearRing::IGeometry{wkbLinearRing}, + wkbCircularStringZ::IGeometry{wkbCircularStringZ}, + wkbCompoundCurveZ::IGeometry{wkbCompoundCurveZ}, + wkbCurvePolygonZ::IGeometry{wkbCurvePolygonZ}, + wkbMultiCurveZ::IGeometry{wkbMultiCurveZ}, + wkbMultiSurfaceZ::IGeometry{wkbMultiSurfaceZ}, + wkbCurveZ::IGeometry{wkbCurveZ}, + wkbSurfaceZ::IGeometry{wkbSurfaceZ}, + wkbPolyhedralSurfaceZ::IGeometry{wkbPolyhedralSurfaceZ}, + wkbTINZ::IGeometry{wkbTINZ}, + wkbTriangleZ::IGeometry{wkbTriangleZ}, + wkbPointM::IGeometry{wkbPointM}, + wkbLineStringM::IGeometry{wkbLineStringM}, + wkbPolygonM::IGeometry{wkbPolygonM}, + wkbMultiPointM::IGeometry{wkbMultiPointM}, + wkbMultiLineStringM::IGeometry{wkbMultiLineStringM}, + wkbMultiPolygonM::IGeometry{wkbMultiPolygonM}, + wkbGeometryCollectionM::IGeometry{wkbGeometryCollectionM}, + wkbCircularStringM::IGeometry{wkbCircularStringM}, + wkbCompoundCurveM::IGeometry{wkbCompoundCurveM}, + wkbCurvePolygonM::IGeometry{wkbCurvePolygonM}, + wkbMultiCurveM::IGeometry{wkbMultiCurveM}, + wkbMultiSurfaceM::IGeometry{wkbMultiSurfaceM}, + wkbCurveM::IGeometry{wkbCurveM}, + wkbSurfaceM::IGeometry{wkbSurfaceM}, + wkbPolyhedralSurfaceM::IGeometry{wkbPolyhedralSurfaceM}, + wkbTINM::IGeometry{wkbTINM}, + wkbTriangleM::IGeometry{wkbTriangleM}, + wkbPointZM::IGeometry{wkbPointZM}, + wkbLineStringZM::IGeometry{wkbLineStringZM}, + wkbPolygonZM::IGeometry{wkbPolygonZM}, + wkbMultiPointZM::IGeometry{wkbMultiPointZM}, + wkbMultiLineStringZM::IGeometry{wkbMultiLineStringZM}, + wkbMultiPolygonZM::IGeometry{wkbMultiPolygonZM}, + wkbGeometryCollectionZM::IGeometry{wkbGeometryCollectionZM}, + wkbCircularStringZM::IGeometry{wkbCircularStringZM}, + wkbCompoundCurveZM::IGeometry{wkbCompoundCurveZM}, + wkbCurvePolygonZM::IGeometry{wkbCurvePolygonZM}, + wkbMultiCurveZM::IGeometry{wkbMultiCurveZM}, + wkbMultiSurfaceZM::IGeometry{wkbMultiSurfaceZM}, + wkbCurveZM::IGeometry{wkbCurveZM}, + wkbSurfaceZM::IGeometry{wkbSurfaceZM}, + wkbPolyhedralSurfaceZM::IGeometry{wkbPolyhedralSurfaceZM}, + wkbTINZM::IGeometry{wkbTINZM}, + wkbTriangleZM::IGeometry{wkbTriangleZM}, + wkbPoint25D::IGeometry{wkbPoint25D}, + wkbLineString25D::IGeometry{wkbLineString25D}, + wkbPolygon25D::IGeometry{wkbPolygon25D}, + wkbMultiPoint25D::IGeometry{wkbMultiPoint25D}, + wkbMultiLineString25D::IGeometry{wkbMultiLineString25D}, + wkbMultiPolygon25D::IGeometry{wkbMultiPolygon25D}, + wkbGeometryCollection25D::IGeometry{wkbGeometryCollection25D}, + wkbUnknown::IGeometry +) + function basetype(gt::OGRwkbGeometryType)::OGRwkbGeometryType wkbGeomType = convert(GDAL.OGRwkbGeometryType, gt) wkbGeomType &= (~0x80000000) # Remove 2.5D flag. diff --git a/src/utils.jl b/src/utils.jl index c27efac6..97222a1b 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -207,7 +207,7 @@ macro ogrerr(code, message) "Unknown error." end - error($message * " ($detailmsg)") + error($(esc(message)) * " ($detailmsg)") end end end diff --git a/test/Project.toml b/test/Project.toml index 134631c2..b1989f37 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -8,6 +8,7 @@ GDAL = "add2ef01-049f-52c4-9ee2-e494f65e021a" GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f" GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f" ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534" +LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" diff --git a/test/test_convert.jl b/test/test_convert.jl index 4f5d209c..90e66334 100644 --- a/test/test_convert.jl +++ b/test/test_convert.jl @@ -11,7 +11,7 @@ const GFT = GeoFormatTypes; point = AG.createpoint(100, 70) json = convert(GFT.GeoJSON, point) @test sprint(print, convert(AG.IGeometry, json)) == - "Geometry: POINT (100 70)" + "IGeometry: POINT (100 70)" kml = convert(GFT.KML, point) gml = convert(GFT.GML, point) wkb = convert(GFT.WellKnownBinary, point) diff --git a/test/test_display.jl b/test/test_display.jl index 3a995fea..a583bb48 100644 --- a/test/test_display.jl +++ b/test/test_display.jl @@ -51,6 +51,22 @@ const AG = ArchGDAL; "pointname (OFTString)" @test sprint(print, AG.getgeomdefn(feature, 0)) == " (wkbPoint)" end + + # Test display for Geometries and IGeometries + @test sprint(print, AG.createpoint(1, 2)) == + "IGeometry: POINT (1 2)" + @test sprint( + print, + AG.createpoint(1, 2), + context = :compact => true, + ) == "IGeometry: wkbPoint" + AG.createpoint(1, 2) do point + @test sprint(print, point) == "Geometry: POINT (1 2)" + end + AG.createpoint(1, 2) do point + @test sprint(print, point, context = :compact => true) == + "Geometry: wkbPoint" + end end AG.read("gdalworkshop/world.tif") do dataset diff --git a/test/test_feature.jl b/test/test_feature.jl index 4c9a2553..ebef9a2f 100644 --- a/test/test_feature.jl +++ b/test/test_feature.jl @@ -241,11 +241,11 @@ const AG = ArchGDAL; @test !AG.isfieldnull(feature, i - 1) @test AG.isfieldsetandnotnull(feature, i - 1) end - @test sprint(print, AG.getgeom(feature)) == "NULL Geometry" + @test sprint(print, AG.getgeom(feature)) == "NULL IGeometry" AG.getgeom(feature) do geom @test sprint(print, geom) == "NULL Geometry" end - @test sprint(print, AG.getgeom(feature, 0)) == "NULL Geometry" + @test sprint(print, AG.getgeom(feature, 0)) == "NULL IGeometry" AG.getgeom(feature, 0) do geom @test sprint(print, geom) == "NULL Geometry" end diff --git a/test/test_featurelayer.jl b/test/test_featurelayer.jl index 15a4ec76..b23f6dca 100644 --- a/test/test_featurelayer.jl +++ b/test/test_featurelayer.jl @@ -95,13 +95,13 @@ const AG = ArchGDAL; ) @test sprint(print, AG.getgeom(newfeature)) == - "Geometry: POINT EMPTY" + "IGeometry: POINT EMPTY" @test sprint(print, AG.getgeom(newfeature, 0)) == - "Geometry: POINT EMPTY" + "IGeometry: POINT EMPTY" @test sprint(print, AG.getgeom(newfeature, 1)) == - "Geometry: LINESTRING EMPTY" + "IGeometry: LINESTRING EMPTY" @test sprint(print, AG.getgeom(newfeature, 2)) == - "Geometry: POLYGON ((0 0,1 1,0 1))" + "IGeometry: POLYGON ((0 0,1 1,0 1))" AG.getgeom(newfeature) do g @test sprint(print, g) == "Geometry: POINT EMPTY" end @@ -152,7 +152,7 @@ const AG = ArchGDAL; @test n == 2 AG.clearspatialfilter!(layer) @test sprint(print, AG.getspatialfilter(layer)) == - "NULL Geometry" + "NULL IGeometry" n = 0 for feature in layer n += 1 diff --git a/test/test_geometry.jl b/test/test_geometry.jl index 417e5ec3..88e1b30b 100644 --- a/test/test_geometry.jl +++ b/test/test_geometry.jl @@ -2,6 +2,8 @@ using Test import GeoInterface, GeoFormatTypes, ArchGDAL; const AG = ArchGDAL const GFT = GeoFormatTypes +using LibGEOS; +const LG = LibGEOS; @testset "test_geometry.jl" begin @testset "Incomplete GeoInterface geometries" begin @@ -653,12 +655,12 @@ const GFT = GeoFormatTypes ) AG.clone(geom3) do geom4 @test sprint(print, AG.clone(geom3)) == - "Geometry: GEOMETRYCOLLECTION (" * + "IGeometry: GEOMETRYCOLLECTION (" * "POINT (2 5 8)," * "POLYGON ((0 0 8," * " ... MPTY)" @test sprint(print, AG.clone(geom4)) == - "Geometry: GEOMETRYCOLLECTION (" * + "IGeometry: GEOMETRYCOLLECTION (" * "POINT (2 5 8)," * "POLYGON ((0 0 8," * " ... MPTY)" @@ -727,8 +729,8 @@ const GFT = GeoFormatTypes @test AG.getgeomtype(AG.getgeom(geom3, 0)) == AG.wkbPoint25D @test AG.getgeomtype(AG.getgeom(geom3, 1)) == AG.wkbPolygon25D @test AG.getgeomtype(AG.getgeom(geom3, 2)) == AG.wkbPolygon25D - @test sprint(print, AG.getgeom(geom3, 3)) == "NULL Geometry" - @test sprint(print, AG.getgeom(AG.IGeometry(), 3)) == "NULL Geometry" + @test sprint(print, AG.getgeom(geom3, 3)) == "NULL IGeometry" + @test sprint(print, AG.getgeom(AG.IGeometry(), 3)) == "NULL IGeometry" AG.getgeom(geom3, 0) do geom4 @test AG.getgeomtype(geom4) == AG.wkbPoint25D end @@ -800,9 +802,25 @@ const GFT = GeoFormatTypes @testset "Cloning NULL geometries" begin geom = AG.IGeometry() @test AG.geomname(geom) === missing - @test sprint(print, AG.clone(geom)) == "NULL Geometry" + @test sprint(print, AG.clone(geom)) == "NULL IGeometry" AG.clone(geom) do g @test sprint(print, g) == "NULL Geometry" end end + + @testset "GeoInterface to IGeometry conversions" begin + wktgeoms = [ + "POINT(0.12345 2.000 0.1)", + "MULTIPOINT (130 240, 130 240, 130 240, 570 240, 570 240, 570 240, 650 240)", + "LINESTRING (130 240, 650 240)", + "MULTILINESTRING ((0 0, 10 10), (0 0, 10 0), (10 0, 10 10))", + "POLYGON((1 1,1 5,5 5,5 1,1 1))", + "MULTIPOLYGON(((0 0,5 10,10 0,0 0),(1 1,1 2,2 2,2 1,1 1),(100 100,100 102,102 102,102 100,100 100)))", + "GEOMETRYCOLLECTION(MULTIPOINT(0 0, 0 0, 1 1),LINESTRING(1 1, 2 2, 2 2, 0 0),POLYGON((5 5, 0 0, 0 2, 2 2, 5 5)))", + ] + for wktgeom in wktgeoms + @test (AG.toWKT ∘ convert)(AG.IGeometry, LG.readgeom(wktgeom)) == + (AG.toWKT ∘ AG.fromWKT)(wktgeom) + end + end end diff --git a/test/test_tables.jl b/test/test_tables.jl index 2e463739..3c5acfbc 100644 --- a/test/test_tables.jl +++ b/test/test_tables.jl @@ -2,6 +2,7 @@ using Test import ArchGDAL; const AG = ArchGDAL; using Tables +using LibGEOS @testset "test_tables.jl" begin @testset "Tables Support" begin @@ -795,5 +796,292 @@ using Tables clean_test_dataset_files() end + + @testset "Table to layer conversion" begin + # Helper functions + toWKT_withmissings(x::Missing) = missing + toWKT_withmissings(x::AG.AbstractGeometry) = AG.toWKT(x) + toWKT_withmissings(x::Any) = x + + function ctv_toWKT( + x::T, + ) where {T<:NTuple{N,AbstractArray{S,D} where S}} where {N,D} + return Tuple(toWKT_withmissings.(x[i]) for i in 1:N) + end + + """ + equals_for_columntables_with_IGeometries(ct1, ct2) + + Compares two `NamedTuple` containing values `<: IGeometry` in the first `ngeom` columns of `ct1`, regarless of key order + + """ + function equals_for_columntables_with_IGeometries(ct1, ct2, ngeom) + # we convert IGeometry values to WKT + (ctv1, ctv2) = ctv_toWKT.(values.((ct1, ct2))) + # we use two index functions to map ctv1 and ctv2 indices to the + # sorted key list indices + (spidx_in, spidx_out) = + sortperm.(([keys(ct1)...], [keys(ct2)...])) + return all([ + sort([keys(ct1)...]) == sort([keys(ct2)...]), + all( + all.([ + isequal.( + ( + # if we are comparing two geometry columns values, we + # convert `nothing` values to `missing`, see note #2 + spidx_out[i] <= ngeom ? + map( + val -> + ( + val === nothing || + val === missing + ) ? missing : val, + ctv1[spidx_in[i]], + ) : ctv1[spidx_in[i]] + ), + ctv2[spidx_out[i]], + ) for i in 1:length(ctv2) + ]), + ), + ]) + end + + """ + nt2layer2nt_equals_nt(nt; force_no_schema=true) + + Takes a NamedTuple, converts it to an IFeatureLayer and compares the NamedTuple + to the one obtained from the IFeatureLayer conversion to table + + _Notes:_ + 1. _Table columns have geometry column first and then field columns as + enforced by `Tables.columnnames`_ + 2. _`nothing` values in geometry column are returned as `missing` from + the NamedTuple roundtrip conversion, since geometry fields do not have the + same distinction between NULL and UNSET values the fields have_ + + """ + function nt2layer2nt_equals_nt( + nt::NamedTuple; + force_no_schema::Bool = false, + )::Bool + force_no_schema ? + layer = AG._fromtable( + nothing, + Tables.rows(nt); + layer_name = "layer", + ) : layer = AG.IFeatureLayer(nt) + (ct_in, ct_out) = Tables.columntable.((nt, layer)) + ngeom = AG.ngeom(layer) + return equals_for_columntables_with_IGeometries( + ct_in, + ct_out, + ngeom, + ) + end + + @testset "Tables with IGeometry" begin + # Test with mixed IGeometry and Float + nt = NamedTuple([ + :point => [AG.createpoint(30, 10), 1.0], + :name => ["point1", "point2"], + ]) + @test_throws ErrorException nt2layer2nt_equals_nt(nt) + + # Test with mixed String and Float64 + nt = NamedTuple([ + :point => [ + AG.createpoint(30, 10), + AG.createlinestring([ + (30.0, 10.0), + (10.0, 30.0), + (40.0, 40.0), + ]), + ], + :name => ["point1", 2.0], + ]) + @test_throws ErrorException nt2layer2nt_equals_nt(nt) + + # Test with Int128 not convertible to OGRFieldType + nt = NamedTuple([ + :point => [ + AG.createpoint(30, 10), + AG.createlinestring([ + (30.0, 10.0), + (10.0, 30.0), + (40.0, 40.0), + ]), + ], + :id => Int128[1, 2], + ]) + @test_throws ErrorException nt2layer2nt_equals_nt(nt) + + # Test with `missing` and `nothing`values + nt = NamedTuple([ + :point => [ + AG.createpoint(30, 10), + nothing, + AG.createpoint(35, 15), + ], + :linestring => [ + AG.createlinestring([ + (30.0, 10.0), + (10.0, 30.0), + (40.0, 40.0), + ]), + AG.createlinestring([ + (35.0, 15.0), + (15.0, 35.0), + (45.0, 45.0), + ]), + missing, + ], + :id => [nothing, "5.1", "5.2"], + :zoom => [1.0, 2.0, 3], + :location => ["Mumbai", missing, "New Delhi"], + :mixedgeom1 => [ + AG.createpoint(30, 10), + AG.createlinestring([ + (30.0, 10.0), + (10.0, 30.0), + (40.0, 40.0), + ]), + AG.createpoint(35, 15), + ], + :mixedgeom2 => [ + AG.createpoint(30, 10), + AG.createlinestring([ + (30.0, 10.0), + (10.0, 30.0), + (40.0, 40.0), + ]), + AG.createmultilinestring([ + [(25.0, 5.0), (5.0, 25.0), (35.0, 35.0)], + [(35.0, 15.0), (15.0, 35.0), (45.0, 45.0)], + ]), + ], + ]) + @test nt2layer2nt_equals_nt(nt; force_no_schema = true) + @test nt2layer2nt_equals_nt(nt) + + # Test with `missing` values + nt = NamedTuple([ + :point => [ + AG.createpoint(30, 10), + missing, + AG.createpoint(35, 15), + ], + :linestring => [ + AG.createlinestring([ + (30.0, 10.0), + (10.0, 30.0), + (40.0, 40.0), + ]), + AG.createlinestring([ + (35.0, 15.0), + (15.0, 35.0), + (45.0, 45.0), + ]), + missing, + ], + :id => [missing, "5.1", "5.2"], + :zoom => [1.0, 2.0, 3], + :location => ["Mumbai", missing, "New Delhi"], + :mixedgeom1 => [ + AG.createpoint(30, 10), + AG.createlinestring([ + (30.0, 10.0), + (10.0, 30.0), + (40.0, 40.0), + ]), + AG.createpoint(35, 15), + ], + :mixedgeom2 => [ + AG.createpoint(30, 10), + AG.createlinestring([ + (30.0, 10.0), + (10.0, 30.0), + (40.0, 40.0), + ]), + AG.createmultilinestring([ + [(25.0, 5.0), (5.0, 25.0), (35.0, 35.0)], + [(35.0, 15.0), (15.0, 35.0), (45.0, 45.0)], + ]), + ], + ]) + @test nt2layer2nt_equals_nt(nt; force_no_schema = true) + @test nt2layer2nt_equals_nt(nt) + end + + @testset "Tables with mixed IGeometry, GeoInterface, WKT/WKB" begin + nt = NamedTuple([ + :point => [ + AG.createpoint(30, 10), + nothing, + AG.createpoint(35, 15), + ], + :linestring => [ + AG.createlinestring([ + (30.0, 10.0), + (10.0, 30.0), + (40.0, 40.0), + ]), + AG.createlinestring([ + (35.0, 15.0), + (15.0, 35.0), + (45.0, 45.0), + ]), + missing, + ], + :id => [nothing, "5.1", "5.2"], + :zoom => [1.0, 2.0, 3], + :location => ["Mumbai", missing, "New Delhi"], + :mixedgeom1 => [ + AG.createpoint(5, 15), + AG.createlinestring([ + (30.0, 10.0), + (10.0, 30.0), + (40.0, 40.0), + ]), + AG.createpoint(35, 15), + ], + :mixedgeom2 => [ + AG.createpoint(10, 20), + AG.createlinestring([ + (30.0, 10.0), + (10.0, 30.0), + (40.0, 40.0), + ]), + AG.createmultilinestring([ + [(25.0, 5.0), (5.0, 25.0), (35.0, 35.0)], + [(35.0, 15.0), (15.0, 35.0), (45.0, 45.0)], + ]), + ], + ]) + + # Test a table conversion with geometries as `GeoInterface.AbstractGeometry` only + nt_native = (; + zip(Symbol.((.*)(String.(keys(nt)), "_GI")), values(nt))..., + ) + nt_GI = (; + zip( + Symbol.((.*)(String.(keys(nt)), "_GI")), + map.( + x -> + typeof(x) <: AG.IGeometry ? + LibGEOS.readgeom(AG.toWKT(x)) : x, + values(nt), + ), + )..., + ) + @test all([ + string( + Tables.columntable(AG.IFeatureLayer(nt_native))[colname], + ) == string( + Tables.columntable(AG.IFeatureLayer(nt_GI))[colname], + ) for colname in keys(nt_native) + ]) + end + end end end diff --git a/test/test_types.jl b/test/test_types.jl index 0489f73a..9bec3417 100644 --- a/test/test_types.jl +++ b/test/test_types.jl @@ -30,6 +30,26 @@ import ImageCore end end + @testset "Convert GeoInterface.AbstractGeometry types to OGRwkbGeometryType" begin + @test convert(AG.OGRwkbGeometryType, GeoInterface.Point) == AG.wkbPoint + @test convert(AG.OGRwkbGeometryType, GeoInterface.MultiPoint) == + AG.wkbMultiPoint + @test convert(AG.OGRwkbGeometryType, GeoInterface.LineString) == + AG.wkbLineString + @test convert(AG.OGRwkbGeometryType, GeoInterface.MultiLineString) == + AG.wkbMultiLineString + @test convert(AG.OGRwkbGeometryType, GeoInterface.Polygon) == + AG.wkbPolygon + @test convert(AG.OGRwkbGeometryType, GeoInterface.MultiPolygon) == + AG.wkbMultiPolygon + @test convert(AG.OGRwkbGeometryType, GeoInterface.AbstractGeometry) == + AG.wkbUnknown + @test_throws ErrorException convert( + AG.OGRwkbGeometryType, + GeoInterface.GeometryCollection, + ) + end + @testset "Testing GDAL Type Methods" begin @testset "GDAL Open Flags" begin @test AG.OF_READONLY | 0x04 == 0x04 diff --git a/test/test_utils.jl b/test/test_utils.jl index ef2e01d4..2e5367fa 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -4,7 +4,10 @@ import ArchGDAL; const AG = ArchGDAL; "Test both that an ErrorException is thrown and that the message is as expected" -eval_ogrerr(err, expected_message) = @test (@test_throws ErrorException AG.@ogrerr err "e:").value.msg == "e: ($expected_message)" +function eval_ogrerr(err, expected_message) + @test (@test_throws ErrorException AG.@ogrerr err "e:").value.msg == + "e: ($expected_message)" +end @testset "test_utils.jl" begin @testset "metadataitem" begin @@ -14,15 +17,29 @@ eval_ogrerr(err, expected_message) = @test (@test_throws ErrorException AG.@ogre @test AG.metadataitem(driver, "DMD_EXTENSIONS") == "tif tiff" end + @testset "gdal error macros" begin + @test_throws ErrorException AG.createlayer() do layer + AG.addfeature(layer) do feature + return AG.setgeom!(feature, 1, AG.createpoint(1, 1)) + end + end + end + @testset "OGR Errors" begin @test isnothing(AG.@ogrerr GDAL.OGRERR_NONE "not an error") eval_ogrerr(GDAL.OGRERR_NOT_ENOUGH_DATA, "Not enough data.") eval_ogrerr(GDAL.OGRERR_NOT_ENOUGH_MEMORY, "Not enough memory.") - eval_ogrerr(GDAL.OGRERR_UNSUPPORTED_GEOMETRY_TYPE, "Unsupported geometry type.") + eval_ogrerr( + GDAL.OGRERR_UNSUPPORTED_GEOMETRY_TYPE, + "Unsupported geometry type.", + ) eval_ogrerr(GDAL.OGRERR_UNSUPPORTED_OPERATION, "Unsupported operation.") eval_ogrerr(GDAL.OGRERR_CORRUPT_DATA, "Corrupt data.") eval_ogrerr(GDAL.OGRERR_FAILURE, "Failure.") - eval_ogrerr(GDAL.OGRERR_UNSUPPORTED_SRS, "Unsupported spatial reference system.") + eval_ogrerr( + GDAL.OGRERR_UNSUPPORTED_SRS, + "Unsupported spatial reference system.", + ) eval_ogrerr(GDAL.OGRERR_INVALID_HANDLE, "Invalid handle.") eval_ogrerr(GDAL.OGRERR_NON_EXISTING_FEATURE, "Non-existing feature.") # OGRERR_NON_EXISTING_FEATURE is the highest error code currently in GDAL. If another one is