diff --git a/Project.toml b/Project.toml index 307ea24..bf77b3c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "LASDatasets" uuid = "cc498e2a-d443-4943-8f26-2a8a0f3c7cdb" authors = ["BenCurran98 "] -version = "0.3.3" +version = "0.4.0" [deps] BufferedStreams = "e1450e63-4bb3-523b-b2a4-4ffa8c0fd77d" diff --git a/docs/src/header.md b/docs/src/header.md index 3bdace3..b841a41 100644 --- a/docs/src/header.md +++ b/docs/src/header.md @@ -38,6 +38,7 @@ You can also modify certain fields in the header, but one should note that for s ```@docs; canonical = false set_las_version! +set_point_format! set_spatial_info! set_point_data_offset! set_point_record_length! diff --git a/docs/src/interface.md b/docs/src/interface.md index 95bd78f..c0980f5 100644 --- a/docs/src/interface.md +++ b/docs/src/interface.md @@ -70,18 +70,47 @@ Alternatively, if you just have the point cloud data as a Table: using StaticArrays using TypedTables -pc = Table(position = rand(SVector{3, Float64}, 10), classification = rand(UIn8, 10)) +pc = Table(position = rand(SVector{3, Float64}, 10), classification = rand(UInt8, 10)) save_las("my_las.las", pc) ``` Note that when you supply just the point cloud outside of a `LASDataset`, *LASDatasets.jl* will automatically construct the appropriate header for you so you don't need to worry about the specifics of appropriate point formats etc. ## Modifying LAS Contents -You can modify point fields in your `LASDataset` by adding new columns or merging in values from an existing vector. +You can modify point fields in your `LASDataset` by adding new columns or merging in values from an existing vector. Additionally, you can add and remove points from a dataset. When adding points, the user is responsible for correctly setting the appropriate fields (e.g. synthetic flags). ```@docs; canonical = false add_column! merge_column! +add_points! +remove_points! +``` + +For example, if you want to add a set of synthetic points to your dataset, you can run: +```julia +las = load_las("my_las.las") +# note - we need to set a synthetic column here for the existing points before we append points with this field +add_column!(las, :synthetic, falses(number_of_points(las))) +synthetic_points = Table(position = rand(SVector{3, Float64}, 5), classification = rand(UInt8, 5), synthetic = trues(5)) +add_points!(las, synthetic_points) +``` + +You can remove points from your data using the `remove_points!` function and specifying the indices of the points you wish to delete (these will be indexing into the list of points in order). E.g. +```julia +remove_points!(las, 11:15) +``` + +Note that you can also modify the contents of your points by acting directly on the tabular pointcloud data. **Note:** this should **not** be used to add/remove points or point fields, since this will cause a conflict between the data in your points and the file header. Intended use is for operations that preserve the number of points and the existing fields. +For example: + +```julia +pc = get_pointcloud(las) + +# shuffle the order of the points based on point positions +pc = pc[sortperm(pc.position)] + +# set all classifications to 0 +pc.classification .= 0 ``` You can also add or remove *(E)VLRs* using the following functions, and set an existing *(E)VLR* as *superseded* if it's an old copy of a record. diff --git a/docs/src/vlrs.md b/docs/src/vlrs.md index 8ece1bf..e050897 100644 --- a/docs/src/vlrs.md +++ b/docs/src/vlrs.md @@ -62,7 +62,7 @@ You can then read the *LAS* data and extract the classification lookup: las = load_las("pc.las") # look for the classification lookup VLR by checking for its user and record IDs lookup_vlr = extract_vlr_type(get_vlrs(las), "LASF_Spec", 0) -lookup = get_data(lookup_vlr[1]) +lookup = get_data(lookup_vlr) ``` ### Text Area Descriptions @@ -150,7 +150,7 @@ So in our example, we can tell the system that records containing data of type ` @register_vlr_type MyType "My Custom Records" collect(1:100) ``` -And now we can save our `MyType` *VLRs* into a *LAS* file in the same way as we did above for the register *VLR* types. Note that you can use the function `extract_vlr_type` on your collection of *VLRs* to pull out any *VLRs* with a specific user ID and record ID. +And now we can save our `MyType` *VLRs* into a *LAS* file in the same way as we did above for the register *VLR* types. Note that you can use the function `extract_vlr_type` on your collection of *VLRs* to pull out the *VLR* with a specific user ID and record ID. ```@docs; canonical = false extract_vlr_type diff --git a/src/LASDatasets.jl b/src/LASDatasets.jl index 6c6babc..3a28805 100644 --- a/src/LASDatasets.jl +++ b/src/LASDatasets.jl @@ -39,7 +39,7 @@ export scan_direction, user_data # header interface export LasHeader, is_standard_gps, is_wkt, file_source_id, global_encoding, las_version, system_id, software_id, creation_day_of_year, creation_year export header_size, point_data_offset, point_record_length, record_format, point_format, number_of_points, number_of_vlrs, number_of_evlrs, evlr_start, spatial_info, scale, num_return_channels -export set_las_version!, set_spatial_info!, set_point_data_offset!, set_point_format!, set_point_record_length!, set_point_record_count!, set_num_vlr!, set_num_evlr!, set_gps_week_time_bit! +export set_las_version!, set_point_format!, set_spatial_info!, set_point_data_offset!, set_point_format!, set_point_record_length!, set_point_record_count!, set_num_vlr!, set_num_evlr!, set_gps_week_time_bit! export set_gps_standard_time_bit!, is_internal_waveform, is_external_waveform, unset_wkt_bit! export set_waveform_external_bit!, set_waveform_internal_bit!, set_synthetic_return_numbers_bit!, unset_synthetic_return_numbers_bit! export set_wkt_bit!, get_number_of_points_by_return, set_number_of_points_by_return!, waveform_record_start @@ -55,7 +55,7 @@ export get_classes, get_description, set_description! export @register_vlr_type, read_vlr_data, extract_vlr_type export LASDataset, get_header, get_pointcloud, get_vlrs, get_evlrs, get_user_defined_bytes, get_unit_conversion -export add_column!, merge_column!, add_vlr!, remove_vlr!, set_superseded! +export add_column!, add_points!, remove_points!, merge_column!, add_vlr!, remove_vlr!, set_superseded! # I/O methods export load_las, load_pointcloud, save_las, load_header, load_vlrs diff --git a/src/dataset.jl b/src/dataset.jl index 1fe8f34..75b1433 100644 --- a/src/dataset.jl +++ b/src/dataset.jl @@ -5,35 +5,32 @@ A wrapper around a LAS dataset. Contains point cloud data in tabular format as w $(TYPEDFIELDS) """ -mutable struct LASDataset +struct LASDataset """The header from the LAS file the points were extracted from""" - const header::LasHeader + header::LasHeader - """LAS point cloud data stored in a Tabular format for convenience""" - const pointcloud::Table - - """Additional user data assigned to each point that aren't standard LAS fields""" - _user_data::Union{Missing, FlexTable} + """Point cloud data stored in a Tabular format for convenience""" + pointcloud::FlexTable """Collection of Variable Length Records from the LAS file""" - const vlrs::Vector{LasVariableLengthRecord} + vlrs::Vector{LasVariableLengthRecord} """Collection of Extended Variable Length Records from the LAS file""" - const evlrs::Vector{LasVariableLengthRecord} + evlrs::Vector{LasVariableLengthRecord} """Extra user bytes packed between the Header block and the first VLR of the source LAS file""" - const user_defined_bytes::Vector{UInt8} + user_defined_bytes::Vector{UInt8} """Unit conversion factors applied to each axis when the dataset is ingested. This is reversed when you save the dataset to keep header/coordinate system information consistent""" - const unit_conversion::SVector{3, Float64} + unit_conversion::SVector{3, Float64} function LASDataset(header::LasHeader, - pointcloud::Table, + pointcloud::AbstractVector{<:NamedTuple}, vlrs::Vector{<:LasVariableLengthRecord}, evlrs::Vector{<:LasVariableLengthRecord}, user_defined_bytes::Vector{UInt8}, unit_conversion::SVector{3, Float64} = NO_CONVERSION) - + pointcloud = FlexTable(pointcloud) # do a few checks to make sure everything is consistent between the header and other data point_format_from_table = get_point_format(pointcloud) point_format_from_header = point_format(header) @@ -59,19 +56,16 @@ mutable struct LASDataset # if points don't have an ID column assigned to them, add it in if :id ∉ columnnames(pointcloud) - pointcloud = Table(pointcloud, id = collect(1:length(pointcloud))) + pointcloud.id = collect(1:length(pointcloud)) end # need to split out our "standard" las columns from user-specific ones cols = collect(columnnames(pointcloud)) these_are_las_cols = cols .∈ Ref(RECOGNISED_LAS_COLUMNS) - las_cols = cols[these_are_las_cols] other_cols = cols[.!these_are_las_cols] - las_pc = Table(NamedTuple{ (las_cols...,) }( (map(col -> getproperty(pointcloud, col), las_cols)...,) )) - make_consistent_header!(header, las_pc, vlrs, evlrs, user_defined_bytes) + make_consistent_header!(header, pointcloud, vlrs, evlrs, user_defined_bytes) - user_pc = isempty(other_cols) ? missing : FlexTable(NamedTuple{ (other_cols...,) }( (map(col -> getproperty(pointcloud, col), other_cols)...,) )) for col ∈ other_cols # account for potentially having an undocumented entry - in this case, don't add an ExtraBytes VLR if col != :undocumented_bytes @@ -82,15 +76,12 @@ mutable struct LASDataset # grab information about the existing ExtraBytes VLRs - need to see if we need to update them or not extra_bytes_vlr = extract_vlr_type(vlrs, LAS_SPEC_USER_ID, ID_EXTRABYTES) - @assert length(extra_bytes_vlr) ≤ 1 "Found multiple Extra Bytes VLRs in LAS file!" - if isempty(extra_bytes_vlr) + if isnothing(extra_bytes_vlr) extra_bytes_vlr = LasVariableLengthRecord(LAS_SPEC_USER_ID, ID_EXTRABYTES, "Extra Bytes", ExtraBytesCollection()) # make sure we add the VLR to our collection and update any header info push!(vlrs, extra_bytes_vlr) header.n_vlr += 1 header.data_offset += sizeof(extra_bytes_vlr) - else - extra_bytes_vlr = extra_bytes_vlr[1] end extra_bytes_data = get_extra_bytes(get_data(extra_bytes_vlr)) user_field_names = Symbol.(name.(extra_bytes_data)) @@ -121,7 +112,7 @@ mutable struct LASDataset end end end - return new(header, las_pc, user_pc, Vector{LasVariableLengthRecord}(vlrs), Vector{LasVariableLengthRecord}(evlrs), user_defined_bytes, unit_conversion) + return new(header, pointcloud, Vector{LasVariableLengthRecord}(vlrs), Vector{LasVariableLengthRecord}(evlrs), user_defined_bytes, unit_conversion) end end @@ -131,14 +122,14 @@ end Create a LASDataset from a pointcloud and optionally vlrs/evlrs/user_defined_bytes, NO header required. """ -function LASDataset(pointcloud::Table; +function LASDataset(pointcloud::AbstractVector{<:NamedTuple}; vlrs::Vector{<:LasVariableLengthRecord} = LasVariableLengthRecord[], evlrs::Vector{<:LasVariableLengthRecord} = LasVariableLengthRecord[], user_defined_bytes::Vector{UInt8} = UInt8[], - scale::Real = POINT_SCALE) + scale::Union{Real, SVector, AxisInfo} = POINT_SCALE) point_format = get_point_format(pointcloud) - spatial_info = get_spatial_info(pointcloud; scale = scale) + spatial_info = get_spatial_info(pointcloud, scale) header = LasHeader(; las_version = lasversion_for_point(point_format), @@ -158,11 +149,7 @@ end Extract point cloud data as a Table from a `LASDataset` `las` """ function get_pointcloud(las::LASDataset) - if ismissing(las._user_data) - return las.pointcloud - else - return Table(las.pointcloud, las._user_data) - end + return las.pointcloud end """ @@ -215,9 +202,11 @@ function Base.show(io::IO, las::LASDataset) println(io, "LAS Dataset") println(io, "\tNum Points: $(length(get_pointcloud(las)))") println(io, "\tPoint Format: $(point_format(get_header(las)))") - println(io, "\tPoint Channels: $(columnnames(las.pointcloud))") - if !ismissing(las._user_data) - println(io, "\tUser Fields: $(columnnames(las._user_data))") + all_cols = columnnames(las.pointcloud) + is_las = collect(all_cols .∈ Ref(RECOGNISED_LAS_COLUMNS)) + println(io, "\tPoint Channels: $(all_cols[is_las])") + if any(.!is_las) + println(io, "\tUser Fields: $(all_cols[.!is_las])") end println(io, "\tVLR's: $(length(get_vlrs(las)))") println(io, "\tEVLR's: $(length(get_evlrs(las)))") @@ -264,12 +253,16 @@ Add a `vlr` into the set of VLRs in a LAS dataset `las`. Note that this will modify the header content of `las`, including updating its LAS version to v1.4 if `vlr` is extended """ function add_vlr!(las::LASDataset, vlr::LasVariableLengthRecord) - if is_extended(vlr) && isempty(get_evlrs(las)) + if is_extended(vlr) && las_version(las) < v"1.4" # evlrs only supported in LAS 1.4 + @warn "Upgrading LAS spec version to 1.4.0 from $(las_version(las)) to support use of EVLRs" set_las_version!(get_header(las), v"1.4") end header = get_header(las) + + existing_vlr = extract_vlr_type(get_vlrs(las), get_user_id(vlr), get_record_id(vlr)) + @assert isnothing(existing_vlr) "We already have a VLR with user ID $(get_user_id(vlr)) and record ID $(get_record_id(vlr)) in the las dataset!" if is_extended(vlr) set_num_evlr!(header, number_of_evlrs(header) + 1) @@ -334,39 +327,51 @@ Add a column with name `column` and set of `values` to a `las` dataset """ function add_column!(las::LASDataset, column::Symbol, values::AbstractVector{T}) where T @assert length(values) == length(las.pointcloud) "Column size $(length(values)) inconsistent with number of points $(length(las.pointcloud))" - check_user_type(T) - if ismissing(las._user_data) - las._user_data = FlexTable(NamedTuple{ (column,) }( (values,) )) - else - # make sure if we're replacing a column we correctly update the header size - if column ∈ columnnames(las._user_data) - las.header.data_record_length -= sizeof(eltype(getproperty(las._user_data, column))) + pointcloud = get_pointcloud(las) + + is_user = column ∉ RECOGNISED_LAS_COLUMNS + if is_user + check_user_type(T) + # need to update our header information and VLRs to track this user column + if column ∈ columnnames(pointcloud) + las.header.data_record_length -= sizeof(eltype(getproperty(pointcloud, column))) end - Base.setproperty!(las._user_data, column, values) - end - las.header.data_record_length += sizeof(T) - vlrs = get_vlrs(las) - extra_bytes_vlrs = extract_vlr_type(vlrs, LAS_SPEC_USER_ID, ID_EXTRABYTES) - @assert length(extra_bytes_vlrs) ≤ 1 "Found $(length(extra_bytes_vlrs)) Extra Bytes VLRs when we can only have a max of 1" - if isempty(extra_bytes_vlrs) - extra_bytes_vlr = LasVariableLengthRecord(LAS_SPEC_USER_ID, ID_EXTRABYTES, "Extra Bytes Records", ExtraBytesCollection()) - # make sure we add it to the dataset to account for offsets in the header etc. - add_vlr!(las, extra_bytes_vlr) - else - extra_bytes_vlr = extra_bytes_vlrs[1] - end - if T <: SVector - # user field arrays have to be saved as sequential extra bytes records with names of the form "column [i]" (zero indexing encouraged) - split_col_name = split_column_name(column, length(T)) - for i ∈ 1:length(T) - add_extra_bytes!(las, split_col_name[i], eltype(T), extra_bytes_vlr) + las.header.data_record_length += sizeof(T) + vlrs = get_vlrs(las) + extra_bytes_vlr = extract_vlr_type(vlrs, LAS_SPEC_USER_ID, ID_EXTRABYTES) + if isnothing(extra_bytes_vlr) + extra_bytes_vlr = LasVariableLengthRecord(LAS_SPEC_USER_ID, ID_EXTRABYTES, "Extra Bytes Records", ExtraBytesCollection()) + # make sure we add it to the dataset to account for offsets in the header etc. + add_vlr!(las, extra_bytes_vlr) end - else - add_extra_bytes!(las, column, T, extra_bytes_vlr) + if T <: SVector + # user field arrays have to be saved as sequential extra bytes records with names of the form "column [i]" (zero indexing encouraged) + split_col_name = split_column_name(column, length(T)) + for i ∈ 1:length(T) + add_extra_bytes!(las, split_col_name[i], eltype(T), extra_bytes_vlr) + end + else + add_extra_bytes!(las, column, T, extra_bytes_vlr) + end + # make sure we keep our offset to our first EVLR consistent now we've crammed more data in + update_evlr_offset!(las) + elseif column ∉ has_columns(point_format(get_header(las))) + # we're adding a new LAS column, which will necessitate changing the point format (and possibly the version) + las_cols = filter(c -> c ∈ RECOGNISED_LAS_COLUMNS, collect(columnnames(pointcloud))) + push!(las_cols, column) + new_format = get_point_format(las_cols) + @warn "Changing point format to $(new_format) to allow the inclusion of LAS column $(column)" + + set_point_format!(las, new_format) + + # special case - if we're adding synthetic points, need to set the synthetic flag in the header + set_synthetic_return_numbers_bit!(las) end - # make sure we keep our offset to our first EVLR consistent now we've crammed more data in - update_evlr_offset!(las) - nothing + + # now actually write the values to the column + Base.setproperty!(pointcloud, column, values) + + return nothing end """ @@ -424,4 +429,76 @@ function add_extra_bytes!(las::LASDataset, col_name::Symbol, ::Type{T}, extra_by add_extra_bytes_to_collection!(get_data(extra_bytes_vlr), col_name, T) header = get_header(las) header.data_offset += sizeof(ExtraBytes) +end + +function set_point_format!(las::LASDataset, ::Type{TPoint}) where {TPoint <: LasPoint} + set_point_format!(get_header(las), TPoint) +end + +""" + $(TYPEDSIGNATURES) + +Add a collection of `points` to a `LASDataset`, `las`. Updates header information to ensure dataset consistency +""" +function add_points!(las::LASDataset, points::AbstractVector{<:NamedTuple}) + pc = get_pointcloud(las) + # make new points a FlexTable so we can add any missing columns with 0 values to ensure consistency + missing_cols = filter(c -> c ∉ columnnames(points), columnnames(pc)) + if !isempty(missing_cols) + @warn "Adding default entries for missing columns $(missing_cols)" + end + # need to make sure columns are in the same order as in pc to avoid errors from TypedTables + new_points = Table(NamedTuple{ (columnnames(pc)...,) }( (map(col -> hasproperty(points, col) ? getproperty(points, col) : zeros(eltype(getproperty(pc, col)), length(points)), columnnames(pc))...,) )) + append!(pc, new_points) + # make sure we update the header info too! + _consolidate_point_header_info!(get_header(las), pc) + return nothing +end + +""" + $(TYPEDSIGNATURES) + +Remove a set of points stored at indices `idxs` from a `las` dataset. Updates header information to ensure consistency +""" +function remove_points!(las::LASDataset, idxs::Union{AbstractUnitRange, AbstractVector{<:Integer}}) + pc = get_pointcloud(las) + deleteat!(pc, idxs) + _consolidate_point_header_info!(get_header(las), pc) + return nothing +end + +# plumb through header util functions to act on a LASDataset for convenience +for header_func ∈ ( + :las_version, + :file_source_id, + :global_encoding, + :system_id, + :software_id, + :creation_day_of_year, + :creation_year, + :header_size, + :point_data_offset, + :point_record_length, + :point_format, + :number_of_points, + :number_of_vlrs, + :number_of_evlrs, + :evlr_start, + :spatial_info, + :scale, + :num_return_channels, + :is_standard_gps, + :is_wkt, + :set_gps_standard_time_bit!, + :is_internal_waveform, + :set_waveform_internal_bit!, + :set_waveform_external_bit!, + :unset_waveform_bits!, + :set_synthetic_return_numbers_bit!, + :unset_synthetic_return_numbers_bit!, + :set_wkt_bit!, + :get_number_of_points_by_return, + :waveform_record_start + ) + @eval $header_func(las::LASDataset) = $header_func(get_header(las)) end \ No newline at end of file diff --git a/src/header.jl b/src/header.jl index cfba31f..a796f4f 100644 --- a/src/header.jl +++ b/src/header.jl @@ -476,15 +476,19 @@ spatial_info(h::LasHeader) = h.spatial_info """ $(TYPEDSIGNATURES) -Get the scale for point positions in a LAS file from a header `h`. Checks consistent scale for ALL axes. +Get the scale for point positions in a LAS file from a header `h` along an `axis` (x, y or z) """ -function scale(h::LasHeader) - - @assert (h.spatial_info.scale.x == h.spatial_info.scale.y) && (h.spatial_info.scale.y == h.spatial_info.scale.z) "We expect all axes to be scaled similarly" - - return h.spatial_info.scale.x +function scale(h::LasHeader, axis::Symbol) + return getproperty(h.spatial_info.scale, axis) end +""" + $(TYPEDSIGNATURES) + +Get the scale for point positions in a LAS file from a header `h` along all axes +""" +scale(h::LasHeader) = h.spatial_info.scale + """ $(TYPEDSIGNATURES) @@ -495,26 +499,53 @@ num_return_channels(h::LasHeader) = las_version(h) ≥ v"1.4" ? 15 : 5 """ $(TYPEDSIGNATURES) -Get the LAS version in a header `h` +Set the LAS specification version in a header `h` to version `v` """ function set_las_version!(h::LasHeader, v::VersionNumber) - if any([ - (v ≤ v"1.1") && (get_point_format_id(point_format(h)) ≥ 2), - (v ≤ v"1.2") && (get_point_format_id(point_format(h)) ≥ 4), - (v ≤ v"1.3") && (get_point_format_id(point_format(h)) ≥ 5), - (v ≤ v"1.4") && (get_point_format_id(point_format(h)) ≥ 11), - v ≥ v"1.5" - ]) - error("Incompatible LAS version $(v) with point format $(point_format(h))") - end + _point_format_version_consistent(v, point_format(h)) old_version = deepcopy(las_version(h)) h.las_version = v h.header_size = get_header_size_from_version(v) h.data_offset += (h.header_size - get_header_size_from_version(old_version)) - if (v ≥ v"1.4") && (get_point_format_id(point_format(h)) ≤ 5) - h.record_count = UInt64(h.legacy_record_count) - h.point_return_count = ntuple(i -> i ≤ 5 ? h.legacy_point_return_count[i] : 0, 15) + set_point_record_count!(h, number_of_points(h)) + set_number_of_points_by_return!(h, get_number_of_points_by_return(h)) + return nothing +end + +function _point_format_version_consistent(v::VersionNumber, ::Type{TPoint}) where {TPoint <: LasPoint} + point_format_id = get_point_format_id(TPoint) + if any([ + (v ≤ v"1.1") && (point_format_id ≥ 2), + (v ≤ v"1.2") && (point_format_id ≥ 4), + (v ≤ v"1.3") && (point_format_id ≥ 5), + (v ≤ v"1.4") && (point_format_id ≥ 11), + v ≥ v"1.5" + ]) + error("Incompatible LAS version $(v) with point format $(point_format_id)") + end +end + +""" + $(TYPEDSIGNATURES) + +Set the point format in a header `h` to a new value, `TPoint` +""" +function set_point_format!(h::LasHeader, ::Type{TPoint}) where {TPoint <: LasPoint} + v = las_version(h) + minimal_required_version = lasversion_for_point(TPoint) + old_format = point_format(h) + old_format_id = get_point_format_id(old_format) + # make sure that the LAS version in the header is consistent with the point format we want - upgrade if necessary, but let the user know + if v < minimal_required_version + @warn "Updating LAS version from $(v) to $(minimal_required_version) to accomodate changing point format from $(old_format) to $(TPoint)" + set_las_version!(h, minimal_required_version) end + _point_format_version_consistent(las_version(h), TPoint) + h.data_format_id = get_point_format_id(TPoint) + h.data_record_length += (byte_size(TPoint) - byte_size(LasPoint{Int(old_format_id)})) + set_point_record_count!(h, number_of_points(h)) + set_number_of_points_by_return!(h, get_number_of_points_by_return(h)) + return nothing end """ @@ -554,13 +585,18 @@ end Set the number of points in a LAS file with a header `header` """ function set_point_record_count!(header::LasHeader, num_points::Integer) - if las_version(header) == v"1.4" + if (las_version(header) == v"1.4") @assert num_points ≤ typemax(UInt64) "Can't have more than $(typemax(UInt64)) points for LAS v1.4" - header.record_count = UInt64(num_points) else @assert num_points ≤ typemax(UInt32) "Can't have more than $(typemax(UInt32)) points for LAS v1.0-1.3" + end + header.record_count = UInt64(num_points) + if get_point_format_id(point_format(header)) ≤ 5 header.legacy_record_count = UInt32(num_points) + else + header.legacy_record_count = zero(UInt32) end + return nothing end """ @@ -582,6 +618,15 @@ function set_num_evlr!(header::LasHeader, n::Integer) header.n_evlr = UInt64(n) end +""" + $(TYPEDSIGNATURES) + +Set the offset (in bytes) into the LAS file where the first EVLR occurs +""" +function set_evlr_start!(header, offset::Integer) + header.evlr_start = offset +end + """If true, GPS Time is standard GPS Time (satellite GPS Time) minus 1e9. If false, GPS Time is GPS Week Time. @@ -725,11 +770,14 @@ function set_number_of_points_by_return!(header::LasHeader, points_per_return::N return_limit = las_version(header) ≥ v"1.4" ? typemax(UInt64) : typemax(UInt32) @assert all(points_per_return .≤ return_limit) "Maximum allowed points per return count is $return_limit" @assert N == num_return_channels(header) "Number of returns $N doesn't match what's in header $(num_return_channels(header))" - if las_version(header) ≥ v"1.4" - header.point_return_count = points_per_return + # note - internally, we store the point return count up to 15 places, even if the version spec only needs 5 - saves having to redefine field types for header + header.point_return_count = ntuple(i -> i ≤ N ? points_per_return[i] : 0, 15) + if get_point_format_id(point_format(header)) ≤ 5 + header.legacy_point_return_count = ntuple(i -> i ≤ N ? points_per_return[i] : 0, 5) else - header.legacy_point_return_count = points_per_return + header.legacy_point_return_count = (0, 0, 0, 0, 0) end + return nothing end """ @@ -772,10 +820,10 @@ function make_consistent_header(pointcloud::AbstractVector{<:NamedTuple}, vlrs::Vector{<:LasVariableLengthRecord}, evlrs::Vector{<:LasVariableLengthRecord}, user_defined_bytes::Vector{UInt8}, - scale::Real) where {TPoint <: LasPoint} + scale::Union{Real, SVector{3, <:Real}, AxisInfo}) where {TPoint <: LasPoint} version = lasversion_for_point(point_format) - spatial_info = get_spatial_info(pointcloud; scale = scale) + spatial_info = get_spatial_info(pointcloud, scale) header = LasHeader(; las_version = version, @@ -804,12 +852,8 @@ function make_consistent_header!(header::LasHeader, point_data_offset = header_size + vlr_size + length(user_defined_bytes) set_point_data_offset!(header, point_data_offset) - set_spatial_info!(header, get_spatial_info(pointcloud; scale = scale(header))) - set_point_record_count!(header, length(pointcloud)) - returns = (:returnnumber ∈ columnnames(pointcloud)) ? pointcloud.returnnumber : ones(Int, length(pointcloud)) - points_per_return = ntuple(r -> count(returns .== r), num_return_channels(header)) - set_number_of_points_by_return!(header, points_per_return) + _consolidate_point_header_info!(header, pointcloud) if !isempty(vlrs) set_num_vlr!(header, length(vlrs)) @@ -833,4 +877,24 @@ function make_consistent_header!(header::LasHeader, # don't want waveform bits set for point formats that don't have waveform data unset_waveform_bits!(header) end + + return nothing +end + +function _consolidate_point_header_info!(header::LasHeader, pointcloud::AbstractVector{<:NamedTuple}) + set_spatial_info!(header, get_spatial_info(pointcloud, scale(header))) + + old_point_count = number_of_points(header) + set_point_record_count!(header, length(pointcloud)) + # make sure we update the EVLR start value if we change the point count + set_evlr_start!(header, evlr_start(header) + (length(pointcloud) - old_point_count)) + returns = (:returnnumber ∈ columnnames(pointcloud)) ? pointcloud.returnnumber : ones(Int, length(pointcloud)) + points_per_return = ntuple(r -> count(returns .== r), num_return_channels(header)) + set_number_of_points_by_return!(header, points_per_return) + if :synthetic ∈ columnnames(pointcloud) + set_synthetic_return_numbers_bit!(header) + else + unset_synthetic_return_numbers_bit!(header) + end + return nothing end \ No newline at end of file diff --git a/src/parse_points.jl b/src/parse_points.jl index 6965a20..cd7d658 100644 --- a/src/parse_points.jl +++ b/src/parse_points.jl @@ -6,8 +6,8 @@ $(METHODLIST) function laspoint end laspoint(::Type{TPoint}, p::TPoint, xyz::SpatialInfo) where TPoint = p -byte_size(point::TPoint) where {N, TPoint <: LasPoint{N}} = byte_size(TPoint) -byte_size(::Type{TPoint}) where {N, TPoint <: LasPoint{N}} = sum(sizeof.(eltype.(fieldtypes(TPoint)))) + +byte_size(::Type{TPoint}) where {TPoint <: LasPoint} = sum(sizeof.(eltype.(fieldtypes(get_point_format(TPoint))))) byte_size(vector::Type{SVector{N,T}}) where {N,T} = sizeof(T) * N @@ -16,9 +16,9 @@ byte_size(vector::Type{SVector{N,T}}) where {N,T} = sizeof(T) * N Get the minimum point format that is compatible with the contents of a point cloud in a `table` """ -function get_point_format(table::AbstractVector{<:NamedTuple}) - - columns = columnnames(table) +get_point_format(table::AbstractVector{<:NamedTuple}) = get_point_format(columnnames(table)) + +function get_point_format(columns::Union{AbstractVector{Symbol}, NTuple{N, Symbol}}) where N # these are the known formats - the ones without wave packets columns_per_format = map(n -> has_columns(get_point_format(LasPoint{n})), collect(0:9)) diff --git a/src/points.jl b/src/points.jl index abcbeec..7dbb022 100644 --- a/src/points.jl +++ b/src/points.jl @@ -227,6 +227,7 @@ Get the concrete point format struct from an abstract `LasPoint` type $(METHODLIST) """ +get_point_format(::Type{TPoint}) where {TPoint <: LasPoint} = TPoint get_point_format(::Type{LasPoint{0}}) = LasPoint0 get_point_format(::Type{LasPoint{1}}) = LasPoint1 get_point_format(::Type{LasPoint{2}}) = LasPoint2 diff --git a/src/read.jl b/src/read.jl index 6a89d25..fa1d5f4 100644 --- a/src/read.jl +++ b/src/read.jl @@ -107,8 +107,7 @@ function read_las_data(io::TIO, required_columns::TTuple=DEFAULT_LAS_COLUMNS; user_defined_bytes = read(io, header.data_offset - pos) extra_bytes_vlr = extract_vlr_type(vlrs, LAS_SPEC_USER_ID, ID_EXTRABYTES) - @assert length(extra_bytes_vlr) ≤ 1 "Found multiple extra bytes columns!" - extra_bytes = isempty(extra_bytes_vlr) ? ExtraBytes[] : get_extra_bytes(get_data(extra_bytes_vlr[1])) + extra_bytes = isnothing(extra_bytes_vlr) ? ExtraBytes[] : get_extra_bytes(get_data(extra_bytes_vlr)) this_format = record_format(header, extra_bytes) xyz = spatial_info(header) diff --git a/src/spatial_info.jl b/src/spatial_info.jl index 9c0b165..ff48cbf 100644 --- a/src/spatial_info.jl +++ b/src/spatial_info.jl @@ -68,7 +68,7 @@ Base.write(io::IO, info::SpatialInfo) = write_struct(io, info) const DEFAULT_SPATIAL_INFO = SpatialInfo(AxisInfo(POINT_SCALE, POINT_SCALE, POINT_SCALE), AxisInfo(0.0, 0.0, 0.0), AxisInfo(Range(Inf, -Inf), Range(Inf, -Inf), Range(Inf, -Inf))) -function bounding_box(points::Vector{SVector{3, T}}) where {T <: Real} +function bounding_box(points::AbstractVector{SVector{3, T}}) where {T <: Real} x_min = typemax(T) x_max = typemin(T) y_min = typemax(T) @@ -103,28 +103,29 @@ function bounding_box(points::Vector{SVector{3, T}}) where {T <: Real} return (; xmin = x_min, ymin = y_min, zmin = z_min, xmax = x_max, ymax = y_max, zmax = z_max) end -function get_spatial_info(points::Vector{SVector{3, T}}; scale::T = T(POINT_SCALE)) where {T <: Real} +get_spatial_info(points::AbstractVector{SVector{3, T}}, scale::T = T(POINT_SCALE)) where {T <: Real} = get_spatial_info(points, AxisInfo{T}(scale, scale, scale)) +get_spatial_info(points::AbstractVector{SVector{3, T}}, scale::SVector{3, T}) where {T <: Real} = get_spatial_info(points, AxisInfo{T}(scale.x, scale.y, scale.z)) + +function get_spatial_info(points::AbstractVector{SVector{3, T}}, scale::AxisInfo{T}) where {T <: Real} bb = bounding_box(points) - x_offset = determine_offset(bb.xmin, bb.xmax, scale) - y_offset = determine_offset(bb.ymin, bb.ymax, scale) - z_offset = determine_offset(bb.zmin, bb.zmax, scale) + x_offset = determine_offset(bb.xmin, bb.xmax, scale.x) + y_offset = determine_offset(bb.ymin, bb.ymax, scale.y) + z_offset = determine_offset(bb.zmin, bb.zmax, scale.z) offset = AxisInfo(x_offset, y_offset, z_offset) - x_min = scale * round((bb.xmin / scale) - 0.5) - y_min = scale * round((bb.ymin / scale) - 0.5) - z_min = scale * round((bb.zmin / scale) - 0.5) - x_max = scale * round((bb.xmax / scale) + 0.5) - y_max = scale * round((bb.ymax / scale) + 0.5) - z_max = scale * round((bb.zmax / scale) + 0.5) - - scale_info = AxisInfo(scale, scale, scale) + x_min = scale.x * round((bb.xmin / scale.x) - 0.5) + y_min = scale.y * round((bb.ymin / scale.y) - 0.5) + z_min = scale.z * round((bb.zmin / scale.z) - 0.5) + x_max = scale.x * round((bb.xmax / scale.x) + 0.5) + y_max = scale.y * round((bb.ymax / scale.y) + 0.5) + z_max = scale.z * round((bb.zmax / scale.z) + 0.5) - return SpatialInfo(scale_info, offset, AxisInfo(Range(x_max, x_min), Range(y_max, y_min), Range(z_max, z_min))) + return SpatialInfo(scale, offset, AxisInfo(Range(x_max, x_min), Range(y_max, y_min), Range(z_max, z_min))) end -get_spatial_info(pc::AbstractVector{<:NamedTuple}; kwargs...) = get_spatial_info(pc.position; kwargs...) +get_spatial_info(pc::Union{AbstractVector{<:NamedTuple}, FlexTable}, scale::Union{Real, SVector, AxisInfo} = POINT_SCALE) = get_spatial_info(pc.position, scale) function determine_offset(min_value, max_value, scale; threshold=10^7) s = round(Int64, ((min_value + max_value) / 2) / scale / threshold) diff --git a/src/vlrs.jl b/src/vlrs.jl index 770de04..2f9a4ad 100644 --- a/src/vlrs.jl +++ b/src/vlrs.jl @@ -259,9 +259,14 @@ is_srs(vlr::LasVariableLengthRecord) = vlr.record_id in ( """ $(TYPEDSIGNATURES) -Extract all VLRs with a `user_id` and `record_id` from a collection of VLRs, `vlrs` +Extract the VLR with a `user_id` and `record_id` from a collection of VLRs, `vlrs` """ function extract_vlr_type(vlrs::Vector{<:LasVariableLengthRecord}, user_id::String, record_id::Integer) - matches_ids = findall(vlr -> (get_user_id(vlr) == user_id) && (get_record_id(vlr) == record_id), vlrs) - return vlrs[matches_ids] + # there should only ever be 1 VLR that has a particular user/record ID combination + matches_ids = findfirst(vlr -> (get_user_id(vlr) == user_id) && (get_record_id(vlr) == record_id), vlrs) + if isnothing(matches_ids) + return nothing + else + return vlrs[matches_ids] + end end \ No newline at end of file diff --git a/src/write.jl b/src/write.jl index e555712..a8f47c7 100644 --- a/src/write.jl +++ b/src/write.jl @@ -96,7 +96,10 @@ function write_las(io::IO, las::LASDataset) this_point_format = point_format(header) xyz = spatial_info(header) - user_fields = ismissing(las._user_data) ? () : filter(c -> c != :undocumented_bytes, columnnames(las._user_data)) + cols = collect(columnnames(pc)) + these_are_las_cols = cols .∈ Ref(RECOGNISED_LAS_COLUMNS) + other_cols = cols[.!these_are_las_cols] + user_fields = filter(c -> c != :undocumented_bytes, other_cols) write(io, header) @@ -151,10 +154,10 @@ function get_record_bytes(records::StructVector{TRecord}, vlrs::Vector{LasVariab if user_field_bytes > 0 # need to write the extra bytes fields in the same order as they appear in the VLR - extra_bytes_vlrs = extract_vlr_type(vlrs, LAS_SPEC_USER_ID, ID_EXTRABYTES) - @assert length(extra_bytes_vlrs) == 1 "Expected to find 1 Extra Bytes VLR, instead found $(length(extra_bytes_vlrs))" + extra_bytes_vlr = extract_vlr_type(vlrs, LAS_SPEC_USER_ID, ID_EXTRABYTES) + @assert !isnothing(extra_bytes_vlr) "Expected to find an Extra Bytes VLR but found none!" # get the order they appear in the VLR - user_field_names = unique(get_base_field_name.(Symbol.(name.(get_extra_bytes(get_data(extra_bytes_vlrs[1])))))) + user_field_names = unique(get_base_field_name.(Symbol.(name.(get_extra_bytes(get_data(extra_bytes_vlr)))))) # create a mapping between the order in the VLR and the order in the record per_record_user_field_names = get_user_field_names(TRecord) user_field_idxs = indexin(user_field_names, collect(per_record_user_field_names)) diff --git a/test/dataset.jl b/test/dataset.jl index 4cdc98c..a5d444a 100644 --- a/test/dataset.jl +++ b/test/dataset.jl @@ -58,9 +58,7 @@ this_header = get_header(las) @test point_record_length(this_header) == LASDatasets.byte_size(LasPoint1) + 10 # our user fields should be populated in the dataset - @test las._user_data isa FlexTable - @test length(las._user_data) == num_points - @test sort(collect(columnnames(las._user_data))) == [:other_thing, :thing] + @test sort(collect(columnnames(las.pointcloud))) == [:classification, :gps_time, :id, :other_thing, :position, :thing] this_pc = get_pointcloud(las) @test this_pc.thing == spicy_pc.thing @test this_pc.other_thing == spicy_pc.other_thing @@ -119,9 +117,11 @@ "Text Area Description", TextAreaDescription("This is the new dataset description") ) - add_vlr!(las, new_desc) # mark the old one as superseded set_superseded!(las, desc) + # and add the new one + add_vlr!(las, new_desc) + vlrs = get_vlrs(las) @test length(vlrs) == 3 @test vlrs[3] == new_desc @@ -183,4 +183,40 @@ remove_vlr!(las, superseded_comment) @test number_of_evlrs(get_header(las)) == 1 @test get_evlrs(las) == [new_commment] + + # test modifying point formats and versions + las = LASDataset(pc) + # if we add a LAS column that isn't covered by the current format, the point format (and possibly LAS version) should be updated in the header + add_column!(las, :overlap, falses(length(pc))) + @test point_format(get_header(las)) == LasPoint6 + @test las_version(get_header(las)) == v"1.4" + # we should be able to add new points in too + new_points = Table( + id = (length(pc) + 1):(length(pc) + 10), + position = rand(SVector{3, Float64}, 10), + classification = rand(UInt8, 10), + gps_time = rand(10), + overlap = falses(10) + ) + add_points!(las, new_points) + pointcloud = get_pointcloud(las) + # check the point contents to make sure our new points are there + @test length(pointcloud) == 110 + # equality checks are annoying when the column order gets switched internally, so check each column individually + for col ∈ columnnames(new_points) + @test getproperty(pointcloud, col)[101:end] == getproperty(new_points, col) + end + orig_pc = deepcopy(pointcloud) + # also make sure our header information is correctly set + @test number_of_points(las) == 110 + @test spatial_info(las) == LASDatasets.get_spatial_info(pointcloud) + + # we can also delete these points if we want + remove_points!(las, 101:110) + @test number_of_points(las) == 100 + pc = get_pointcloud(las) + for col ∈ columnnames(pc) + @test getproperty(pc, col) == getproperty(orig_pc, col)[1:100] + end + end \ No newline at end of file diff --git a/test/header.jl b/test/header.jl index 5c970aa..5bc3c19 100644 --- a/test/header.jl +++ b/test/header.jl @@ -6,7 +6,9 @@ creation_year = UInt16(2024), data_format_id = 0x01, data_record_length = UInt16(28), - record_count = UInt64(100) + record_count = UInt64(100), + point_return_count = UInt64.((100, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)), + legacy_point_return_count = UInt32.((100, 0, 0, 0, 0)) ) # make sure we can get and set properties ok @@ -27,6 +29,35 @@ @test number_of_points(header) == 100 set_point_record_count!(header, 150) @test number_of_points(header) == 150 + @test get_number_of_points_by_return(header) == (100, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) + + # modifying record type and las version + set_las_version!(header, v"1.3") + @test las_version(header) == v"1.3" + @test header_size(header) == 235 + @test point_data_offset(header) == 236 + @test point_format(header) == LasPoint1 + set_point_format!(header, LasPoint0) + @test point_format(header) == LasPoint0 + # should now be 20 bytes (for las point 0) plus 2 undocumented extra bytes + @test point_record_length(header) == 22 + # the point return counts should match the legacy counts here + @test get_number_of_points_by_return(header) == (100, 0, 0, 0, 0) + + # if we request setting to a higher point format than the LAS version in the header supports, the version should be updated to the minimal compatible one + set_point_format!(header, LasPoint7) + @test las_version(header) == v"1.4" + @test header_size(header) == 375 + @test point_data_offset(header) == 376 + + # now, we should get the same return counts, but set the legacy counts internally to 0 + @test get_number_of_points_by_return(header) == (100, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) + @test header.legacy_point_return_count == (0, 0, 0, 0, 0) + + set_point_format!(header, LasPoint1) + @test point_format(header) == LasPoint1 + @test point_record_length(header) == 30 + @test number_of_vlrs(header) == 0 set_num_vlr!(header, 2) @@ -36,7 +67,6 @@ @test number_of_evlrs(header) == 1 @test num_return_channels(header) == 15 - @test get_number_of_points_by_return(header) == Tuple(zeros(Int, 15)) points_per_return = ntuple(i -> 10, 15) set_number_of_points_by_return!(header, points_per_return) @test get_number_of_points_by_return(header) == points_per_return