diff --git a/Project.toml b/Project.toml index 023b724..4106268 100644 --- a/Project.toml +++ b/Project.toml @@ -1,9 +1,10 @@ name = "LazIO" uuid = "c3605908-9f0f-11e8-0a72-0d361c15a277" authors = ["Maarten Pronk "] -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" diff --git a/src/dataset.jl b/src/dataset.jl index e473799..ad64b06 100644 --- a/src/dataset.jl +++ b/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}}() diff --git a/src/laszip_h.jl b/src/laszip_h.jl index 260768d..8821d38 100644 --- a/src/laszip_h.jl +++ b/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) diff --git a/src/table.jl b/src/table.jl index 281cbca..3756fc4 100644 --- a/src/table.jl +++ b/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 diff --git a/test/testio.jl b/test/testio.jl index 20fe572..9933f81 100644 --- a/test/testio.jl +++ b/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