Enable Feature by using a geometry field instead of separate x y z.

Author: evetion

src/geointerface.jl

@@ -2,32 +2,38 @@ import GeoInterface as GI
 
 GeoInterface.isgeometry(::Type{<:Point}) = true
 GeoInterface.geomtrait(::Point) = PointTrait()
-# GeoInterface.trait(::Point) = FeatureTrait()
+GeoInterface.trait(::Point) = FeatureTrait()
 
 GeoInterface.ncoord(::PointTrait, p::Point) = 3
-GeoInterface.getcoord(::PointTrait, p::Point) = [p.x, p.y, p.z]
-GeoInterface.getcoord(::PointTrait, p::Point, i) = getfield(p, i)
+GeoInterface.getcoord(::PointTrait, p::Point) = p.geometry
+GeoInterface.getcoord(::PointTrait, p::Point, i) = p.geometry[i]
 
-GeoInterface.x(::PointTrait, p::Point) = p.x
-GeoInterface.y(::PointTrait, p::Point) = p.y
-GeoInterface.z(::PointTrait, p::Point) = p.z
+GeoInterface.x(::PointTrait, p::Point) = p.geometry[1]
+GeoInterface.y(::PointTrait, p::Point) = p.geometry[2]
+GeoInterface.z(::PointTrait, p::Point) = p.geometry[3]
 
 GeoInterface.isfeature(feat::Type{<:Point}) = true
 GeoInterface.properties(feat::T) where {T<:Point} = NamedTuple(zip(fieldnames(T), getfield(feat, n) for n in fieldnames(T)))
-GeoInterface.geometry(feat::Point) = feat
+GeoInterface.geometry(feat::Point) = feat.geom
 
 GeoInterface.isgeometry(::Type{<:Dataset}) = true
 GeoInterface.geomtrait(::Dataset) = MultiPointTrait()
-# GeoInterface.trait(::Dataset) = FeatureCollectionTrait()
+GeoInterface.trait(::Dataset) = FeatureCollectionTrait()
 GeoInterface.ngeom(::MultiPointTrait, ds::Dataset) = length(ds)
 GeoInterface.getgeom(::MultiPointTrait, ds::Dataset) = ds
 GeoInterface.getgeom(::MultiPointTrait, ds::Dataset, i) = ds[i]
 
-# GeoInterface.crs(geomtrait(geom), geom::customgeom)::GeoFormatTypes.GeoFormat}
+# GeoInterface.crs(geomtrait(geom), geom::customgeom)::GeoFormatTypes.GeoFormat
 GeoInterface.extent(::MultiPointTrait, ds::Dataset) = Extent(X=(ds.header.min_x, ds.header.max_x), Y=(ds.header.min_y, ds.header.max_y), Z=(ds.header.min_z, ds.header.max_z))
 
 GeoInterface.isfeaturecollection(::Type{Dataset}) = true
 GeoInterface.getfeature(::MultiPointTrait, ds::Dataset, i) = ds[i]
 GeoInterface.getfeature(::MultiPointTrait, ds::Dataset) = ds
 GeoInterface.nfeature(::MultiPointTrait, ds::Dataset) = length(ds)
-GeoInterface.geometrycolumns(::Dataset) = ()
+GeoInterface.geometrycolumns(::Dataset) = (:geometry,)
+
+GeoInterface.isgeometry(::Type{<:AbstractVector{<:Point}}) = true
+GeoInterface.geomtrait(::AbstractVector{<:Point}) = MultiPointTrait()
+GeoInterface.ngeom(::MultiPointTrait, ds::AbstractVector{<:Point}) = length(ds)
+GeoInterface.getgeom(::MultiPointTrait, ds::AbstractVector{<:Point}) = ds
+GeoInterface.getgeom(::MultiPointTrait, ds::AbstractVector{<:Point}, i) = ds[i]

src/point.jl

@@ -1,8 +1,6 @@
 abstract type Point end
 struct Point0 <: Point
-    x::Float64
-    y::Float64
-    z::Float64
+    geometry::SVector{3,Float64}
     intensity::UInt16
     return_number::UInt8
     number_of_returns::UInt8
@@ -17,9 +15,8 @@ struct Point0 <: Point
 end
 
 function Point0(rp, am)
-    x, y, z = am((rp.X, rp.Y, rp.Z))
     Point0(
-        x, y, z, rp.intensity,
+        am((rp.X, rp.Y, rp.Z)), rp.intensity,
         return_number(rp),
         number_of_returns(rp),
         scan_direction(rp),
@@ -33,9 +30,7 @@ function Point0(rp, am)
 end
 
 struct Point1 <: Point
-    x::Float64
-    y::Float64
-    z::Float64
+    geometry::SVector{3,Float64}
     intensity::UInt16
     return_number::UInt8
     number_of_returns::UInt8
@@ -51,9 +46,8 @@ struct Point1 <: Point
 end
 
 function Point1(rp, am)
-    x, y, z = am((rp.X, rp.Y, rp.Z))
     Point1(
-        x, y, z, rp.intensity,
+        am((rp.X, rp.Y, rp.Z)), rp.intensity,
         return_number(rp),
         number_of_returns(rp),
         scan_direction(rp),
@@ -68,9 +62,7 @@ function Point1(rp, am)
 end
 
 struct Point2 <: Point
-    x::Float64
-    y::Float64
-    z::Float64
+    geometry::SVector{3,Float64}
     intensity::UInt16
     return_number::UInt8
     number_of_returns::UInt8
@@ -88,9 +80,8 @@ struct Point2 <: Point
 end
 
 function Point2(rp, am)
-    x, y, z = am((rp.X, rp.Y, rp.Z))
     Point2(
-        x, y, z, rp.intensity,
+        am((rp.X, rp.Y, rp.Z)), rp.intensity,
         return_number(rp),
         number_of_returns(rp),
         scan_direction(rp),
@@ -108,9 +99,7 @@ end
 
 
 struct Point3 <: Point
-    x::Float64
-    y::Float64
-    z::Float64
+    geometry::SVector{3,Float64}
     intensity::UInt16
     return_number::UInt8
     number_of_returns::UInt8
@@ -129,9 +118,8 @@ struct Point3 <: Point
 end
 
 function Point3(rp, am)
-    x, y, z = am((rp.X, rp.Y, rp.Z))
     Point3(
-        x, y, z, rp.intensity,
+        am((rp.X, rp.Y, rp.Z)), rp.intensity,
         return_number(rp),
         number_of_returns(rp),
         scan_direction(rp),

src/table.jl

@@ -48,7 +48,11 @@ function Base.setproperty!(p::RawPoint, name, value, header)
         p.return_number = UInt8(value) << 7 | UInt8(LasIO.scan_direction(p)) << 6 | LasIO.number_of_returns(p) << 3 | LasIO.return_number(p)
     elseif name == :point_source_id
         p.point_source_ID = value
-    else
+    elseif name == :geometry
+        p.X = round(Int32, (value[1] - header.x_offset) / header.x_scale_factor)
+        p.Y = round(Int32, (value[2] - header.y_offset) / header.y_scale_factor)
+        p.Z = round(Int32, (value[3] - header.z_offset) / header.z_scale_factor)
+    elseif name in fieldnames(typeof(p))
         setproperty!(p, name, value)
     end
 end
@@ -58,7 +62,7 @@ function write(fn::AbstractString, table, bbox; scalex=0.01, scaley=0.01, scalez
 
     schema = Tables.schema(table)
     isnothing(schema) && error("A Schema is required")
-    all(name in column_names for name in schema.names) || error("Can't map all columns to RawPoint")
+    all(name in column_names for name in schema.names) || @warn("Can't map all columns to RawPoint")
 
     rows = Tables.rows(table)
 

test/testio.jl

@@ -92,7 +92,7 @@ end
 
     @test LazIO.readstring(ds.header.system_identifier) == "Laser shooter, pew pew!"
     @test length(ds) == 3
-    @test first(ds).x == 11000.01
+    @test first(ds).geometry[1] == 11000.01
     @test first(ds).classification == LazIO.classes.ground
 end
 
@@ -169,9 +169,9 @@ end
     p2 = first(ds_out)
     header1 = ds.header
     header2 = ds_out.header
-    @test p1.x === p2.x
-    @test p1.y === p2.y
-    @test p1.z === p2.z
+    @test p1.geometry[1] === p2.geometry[1]
+    @test p1.geometry[2] === p2.geometry[2]
+    @test p1.geometry[3] === p2.geometry[3]
     @test p1.intensity === p2.intensity
     @test header1.number_of_point_records === header2.number_of_point_records
 
@@ -197,9 +197,9 @@ end
     p2 = first(ds_out)
     header1 = ds.header
     header2 = ds_out.header
-    @test p1.x === p2.x
-    @test p1.y === p2.y
-    @test p1.z === p2.z
+    @test p1.geometry[1] === p2.geometry[1]
+    @test p1.geometry[2] === p2.geometry[2]
+    @test p1.geometry[3] === p2.geometry[3]
     @test p1.intensity === p2.intensity
     @test header1.number_of_point_records === header2.number_of_point_records
 
@@ -212,6 +212,9 @@ end
     ds = LazIO.open(testfile_str)
     GeoInterface.testgeometry(ds)
     GeoInterface.testgeometry(ds[1])
+    GeoInterface.testgeometry(collect(ds)[1:2])
     @test_broken GeoInterface.testfeature(ds[1])
     @test_broken GeoInterface.testfeaturecollection(ds)
+
+    @test GeoInterface.z(GeoInterface.getpoint(collect(ds)[1:2], 1)) == 846.66
 end