Added Table write interface for XYZ fields, includes conversion. (#26)

* Added Table write interface for XYZ fields, includes conversion.

* Added conversion of other fields as well.

Project.toml

@@ -1,9 +1,10 @@
 name = "LazIO"
 uuid = "c3605908-9f0f-11e8-0a72-0d361c15a277"
 authors = ["Maarten Pronk <[email protected]>"]
-version = "0.2.1"
+version = "0.3.0"
 
 [deps]
+Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 FixedPointNumbers = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
 LASzip_jll = "8372b9c3-1e34-5cc3-bfab-1a98e101de11"

src/dataset.jl

@@ -10,6 +10,10 @@ function Base.show(io::IO, ds::LazDataset)
     println(io, "LazDataset of $(ds.filename) with $n points.")
 end
 
+function bounds(ds::LazDataset)
+    bounds(ds.header)
+end
+
 function open(f::AbstractString)
     # Setup laszip reader
     laszip_reader = Ref{Ptr{Cvoid}}()

src/laszip_h.jl

@@ -1,6 +1,5 @@
-# Automatically generated using Clang.jl wrap_c, version 0.0.0
-
 using Parameters
+using Dates
 
 @with_kw mutable struct LazGeoKey
     key_id::UInt16 = UInt16(0)
@@ -11,10 +10,10 @@ end
 
 @with_kw mutable struct LazVLR
     reserved::UInt16 = UInt16(0)
-    user_id::NTuple{16,UInt8} = ntuple(i->UInt8(0x0), 16)
+    user_id::NTuple{16,UInt8} = ntuple(i -> UInt8(0x0), 16)
     record_id::UInt16 = UInt16(0)
     record_length_after_header::UInt16 = UInt16(0)
-    description::NTuple{32,UInt8} = ntuple(i->UInt8(0x0), 32)
+    description::NTuple{32,UInt8} = ntuple(i -> UInt8(0x0), 32)
     data::Ptr{UInt8} = pointer("")
 end
 
@@ -33,27 +32,27 @@ Base.convert(::Type{LasIO.LasVariableLengthRecord}, vlr::LazVLR) =
     project_ID_GUID_data_1::UInt32 = UInt32(0)
     project_ID_GUID_data_2::UInt16 = UInt16(0)
     project_ID_GUID_data_3::UInt16 = UInt16(0)
-    project_ID_GUID_data_4::NTuple{8,UInt8} = ntuple(i->UInt8(20), 8)
+    project_ID_GUID_data_4::NTuple{8,UInt8} = ntuple(i -> UInt8(20), 8)
     # project_ID_GUID_data_4::Array{UInt8, 1} = Array{UInt8, 1}(zeros(0, 8))
-    version_major::UInt8 = UInt8(0)
-    version_minor::UInt8 = UInt8(0)
-    system_identifier::NTuple{32,UInt8} = ntuple(i->UInt8(20), 32)
+    version_major::UInt8 = UInt8(1)
+    version_minor::UInt8 = UInt8(2)
+    system_identifier::NTuple{32,UInt8} = ntuple(i -> UInt8(20), 32)
     # system_identifier::Array{UInt8} = Array{UInt8}(32)
-    generating_software::NTuple{32,UInt8} = ntuple(i->UInt8(20), 32)
+    generating_software::NTuple{32,UInt8} = ntuple(i -> UInt8(20), 32)
     # generating_software::Array{UInt8, 1} = Array{UInt8, 1}(zeros(0, 32))
-    file_creation_day::UInt16 = UInt16(0)
-    file_creation_year::UInt16 = UInt16(0)
-    header_size::UInt16 = UInt16(0)
-    offset_to_point_data::UInt32 = UInt32(0)
+    file_creation_day::UInt16 = UInt16((today() - Date(year(today()))).value)
+    file_creation_year::UInt16 = UInt16(year(today()))
+    header_size::UInt16 = UInt16(227)
+    offset_to_point_data::UInt32 = UInt32(227)
     number_of_variable_length_records::UInt32 = UInt32(0)
     point_data_format::UInt8 = UInt8(0)
-    point_data_record_length::UInt16 = UInt16(0)
+    point_data_record_length::UInt16 = UInt16(20)
     number_of_point_records::UInt32 = UInt32(0)
-    number_of_points_by_return::NTuple{5,UInt32} = ntuple(i->UInt32(0), 5)
+    number_of_points_by_return::NTuple{5,UInt32} = ntuple(i -> UInt32(0), 5)
     # number_of_points_by_return::Array{UInt32, 1} = Array{UInt32, 1}(zeros(0, 5))
-    x_scale_factor::Float64 = Float64(0.0)
-    y_scale_factor::Float64 = Float64(0.0)
-    z_scale_factor::Float64 = Float64(0.0)
+    x_scale_factor::Float64 = Float64(1.0)
+    y_scale_factor::Float64 = Float64(1.0)
+    z_scale_factor::Float64 = Float64(1.0)
     x_offset::Float64 = Float64(0.0)
     y_offset::Float64 = Float64(0.0)
     z_offset::Float64 = Float64(0.0)
@@ -67,7 +66,7 @@ Base.convert(::Type{LasIO.LasVariableLengthRecord}, vlr::LazVLR) =
     start_of_first_extended_variable_length_record::UInt64 = UInt64(0)
     number_of_extended_variable_length_records::UInt32 = UInt32(0)
     extended_number_of_point_records::UInt64 = UInt64(0)
-    extended_number_of_points_by_return::NTuple{15,UInt64} = ntuple(i->UInt64(0), 15)
+    extended_number_of_points_by_return::NTuple{15,UInt64} = ntuple(i -> UInt64(0), 15)
     # extended_number_of_points_by_return::Array{UInt64, 1} = Array{UInt64, 1}(zeros(0, 15))
     user_data_in_header_size::UInt32 = UInt32(0)
     user_data_in_header::Ptr{UInt8} = pointer("")
@@ -76,6 +75,10 @@ Base.convert(::Type{LasIO.LasVariableLengthRecord}, vlr::LazVLR) =
     user_data_after_header::Ptr{UInt8} = pointer("")
 end
 
+function bounds(h::LazHeader)
+    (min_x = h.min_x, max_x = h.max_x, min_y = h.min_y, max_y = h.max_y, min_z = h.min_z, max_z = h.max_z)
+end
+
 @with_kw mutable struct LazPoint
     X::Int32 = Int32(0)
     Y::Int32 = Int32(0)
@@ -119,14 +122,47 @@ end
     extended_classification::UInt8 = UInt8(0)
     extended_return_number::UInt8 = UInt8(0)
     # extended_number_of_returns::UInt8 = UInt8(0)
-    dummy::NTuple{7,UInt8} = ntuple(i->UInt8(0), 7)
+    dummy::NTuple{7,UInt8} = ntuple(i -> UInt8(0), 7)
     gps_time::Float64 = Float64(0.0)
-    rgb::NTuple{4,UInt16} = ntuple(i->UInt16(0), 4)
-    wave_packet::NTuple{29,UInt8} = ntuple(i->UInt8(0), 29)
+    rgb::NTuple{4,UInt16} = ntuple(i -> UInt16(0), 4)
+    wave_packet::NTuple{29,UInt8} = ntuple(i -> UInt8(0), 29)
     num_extra_bytes::Int32 = Int32(0)
     extra_bytes::Ptr{UInt8} = pointer("")
 end
 
+const classes = (created = 0,
+    unclassified = 1,
+    ground = 2,
+    low_vegation = 3,
+    medium_vegation = 4,
+    high_vegetation = 5,
+    building = 6,
+    noise = 7,
+    key_point = 8,
+    water = 9,
+    overlap = 12,)
+const user_defined_class = 31 # actually reserved
+const classes_extended = (created = 0,
+    unclassified = 1,
+    ground = 2,
+    low_vegation = 3,
+    medium_vegation = 4,
+    high_vegetation = 5,
+    building = 6,
+    noise = 7,
+    reserved = 8,
+    water = 9,
+    rail = 10,
+    road = 11,
+    reserved_ = 12,
+    wire_guard = 13,
+    wire_conductor = 14,
+    tower_transmission = 15,
+    wire_construct = 16,
+    bridge = 17,
+    noise_high = 18,)
+const user_defined_class_extended = 64
+
 LasIO.return_number(p::LazIO.LazPoint) = (p.return_number & 0b00000111)
 LasIO.number_of_returns(p::LazIO.LazPoint) = (p.return_number & 0b00111000) >> 3
 LasIO.scan_direction(p::LazIO.LazPoint) = Bool((p.return_number & 0b01000000) >> 6)

src/table.jl

@@ -3,4 +3,69 @@ using Tables
 Tables.istable(::Type{<:LazDataset}) = true
 Tables.rowaccess(::Type{<:LazDataset}) = true
 Tables.rows(ds::LazDataset) = ds
+# TODO Only support fields that are there based on las version
 Tables.schema(ds::LazDataset) = Tables.Schema(map(Symbol, fieldnames(LazIO.LazPoint)), Tuple{fieldtypes(LazIO.LazPoint)...})
+
+const column_names = (fieldnames(LazPoint)..., :number_of_returns)
+
+
+"""Determine offset as implemented by LasTools."""
+function determine_offset(min_value, max_value, scale; threshold=10^7)
+    s = round(Int64, ((min_value + max_value) / 2) / scale / threshold)
+    s *= threshold * scale
+
+    # Try to convert back and forth and check overflow
+    (muladd(round(Int32, (min_value - s) / scale), scale, s) > 0) == (min_value > 0) || error("Can't fit offset in this scale, try to coarsen it.")
+    (muladd(round(Int32, (max_value - s) / scale), scale, s) > 0) == (max_value > 0) || error("Can't fit offset in this scale, try to coarsen it.")
+end
+
+"""Correctly set fields that require conversion or packing."""
+function Base.setproperty!(p::LazPoint, name, value, header)
+    if name == :X && typeof(value) != Int32
+        p.X = round(Int32, (value - header.x_offset) / header.x_scale_factor)
+    elseif name == :Y && typeof(value) != Int32
+        p.Y = round(Int32, (value - header.y_offset) / header.y_scale_factor)
+    elseif name == :Z && typeof(value) != Int32
+        p.Z = round(Int32, (value - header.z_offset) / header.z_scale_factor)
+    elseif name == :classification && typeof(value) != UInt8
+        classid = get(classes, Symbol(value), user_defined_class)
+        p.classification = UInt8(LasIO.withheld(p)) << 7 | UInt8(LasIO.key_point(p)) << 6 | UInt8(LasIO.synthetic(p)) << 5 | UInt8(classid)
+    elseif name == :gps_time && typeof(value) != Float64
+        p.gps_time = LasIO.gps_time(value)
+    elseif name == :return_number
+        p.return_number = UInt8(LasIO.edge_of_flight_line(p)) << 7 | UInt8(LasIO.scan_direction(p)) << 6 | LasIO.number_of_returns(p) << 3 | UInt8(value)
+    elseif name == :number_of_returns
+        p.return_number = UInt8(LasIO.edge_of_flight_line(p)) << 7 | UInt8(LasIO.scan_direction(p)) << 6 | UInt8(value) << 3 | LasIO.return_number(p)
+    else
+        setproperty!(p, name, value)
+    end
+end
+
+function write(fn::AbstractString, table, bbox; scale=0.01)
+
+    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 LazPoint")
+
+    rows = Tables.rows(table)
+
+    header = LazHeader()
+    header.x_offset = determine_offset(bbox.min_x, bbox.max_x, scale)
+    header.y_offset = determine_offset(bbox.min_y, bbox.max_y, scale)
+    header.z_offset = determine_offset(bbox.min_z, bbox.max_z, scale)
+    header.x_scale_factor = scale
+    header.y_scale_factor = scale
+    header.z_scale_factor = scale
+
+    p = LazPoint()
+
+    LazIO.write(fn, header) do io
+        for row in rows
+            Tables.eachcolumn(schema, row) do value, i, name
+                setproperty!(p, name, value, header)
+            end
+            LazIO.writepoint(io, p)
+        end
+    end
+    fn
+end

test/testio.jl

@@ -3,6 +3,7 @@ using LazIO
 using LasIO
 using Test
 using Tables
+using Dates
 
 workdir = @__DIR__
 lasio_testdir = joinpath(dirname(pathof(LasIO)), "..", "test")
@@ -26,9 +27,9 @@ header, pointdata_all = load(testfile)
 end
 
 @testset "Range indexing" begin
-    _, pointdata_276 = load(testfile, range = 276)
-    _, pointdata_array = load(testfile, range = [1,276,277,497536])
-    _, pointdata_colon = load(testfile, range = :)
+    _, pointdata_276 = load(testfile, range=276)
+    _, pointdata_array = load(testfile, range=[1,276,277,497536])
+    _, pointdata_colon = load(testfile, range=:)
 
     @test pointdata_all[1] == pointdata_colon[1] == pointdata_array[1]
     @test pointdata_276[1] == pointdata_colon[276] == pointdata_array[2]
@@ -69,7 +70,27 @@ end
     @test Tables.rowaccess(LazIO.LazDataset)
     @test first(ds) isa LazIO.LazPoint
     @inferred first(Tables.rows(ds))
+
+    LazIO.write("test_table.laz", ds, LazIO.bounds(ds), scale=0.1)
     close(ds)
+
+    manual_fn = "test_table_manual.laz"
+    table = (X = [11000.01, 12000., 32000.],
+             Y = [11000., 12000., 32000.],
+             Z = [11000., 12000., 32000.01],
+             classification = ["ground", "test", "unclassified"],
+             gps_time = fill(now(), 3),
+             return_number = [1,2,1],
+             number_of_returns = [1,2,3])
+    bounds = (min_x = 11000.01, max_x = 32000., min_y = 11000., max_y = 32000., min_z = 11000., max_z = 32000.01)
+    LazIO.write(manual_fn, table, bounds, scale=0.01)
+
+    ds = LazIO.open(manual_fn)
+
+    @test length(ds) == 3
+    @test first(ds).X == 1099901
+    @test muladd(first(ds).X, ds.header.x_scale_factor, ds.header.x_offset) ≈ 11000.01
+    @test LasIO.classification(first(ds)) == LazIO.classes.ground
 end
 
 @testset "Writing" begin