Flexibly Modify Point Contents

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!, merge_column!, add_vlr!, remove_vlr!, set_superseded!
 
 # I/O methods
 export load_las, load_pointcloud, save_las, load_header, load_vlrs

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,
                         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
@@ -121,7 +115,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
 
@@ -158,11 +152,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 +205,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)))")
@@ -334,39 +326,52 @@ 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_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
-    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)
     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,58 @@ 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)
+    @assert all(columnnames(points) .∈ Ref(columnnames(pc))) "Point fields being appended $(columnnames(points)) must match those present in the dataset, $(columnnames(pc))"
+    append!(pc, points)
+    # make sure we update the header info too!
+    _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

src/header.jl

@@ -495,26 +495,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 +581,16 @@ 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)
     end
+    return nothing
 end
 
 """
@@ -725,10 +755,12 @@ 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
 end
 
@@ -804,12 +836,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 +861,16 @@ 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 = 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)
+    return nothing
 end