diff --git a/Project.toml b/Project.toml index d4f35b1..5d57a32 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "LAS" uuid = "cc498e2a-d443-4943-8f26-2a8a0f3c7cdb" authors = ["BenCurran98 "] -version = "0.2.0" +version = "0.2.1" [deps] ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3" diff --git a/src/LAS.jl b/src/LAS.jl index 6b4e8a7..ff984b4 100644 --- a/src/LAS.jl +++ b/src/LAS.jl @@ -53,7 +53,7 @@ export get_horizontal_unit, get_vertical_unit, get_wkt_string 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 +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! # I/O methods diff --git a/src/constants.jl b/src/constants.jl index b3cd013..c608c03 100644 --- a/src/constants.jl +++ b/src/constants.jl @@ -18,6 +18,8 @@ const UNIT_CONVERSION = Dict{String, Float64}( "us-in" => 0.0254000508, "cm" => 100.0, "mm" => 1000.0) +const NO_CONVERSION = SVector{3, Float64}(1.0, 1.0, 1.0) + const LAS_SPEC_USER_ID = "LASF_Spec" const LAS_PROJ_USER_ID = "LASF_Projection" const ID_GEOKEYDIRECTORYTAG = UInt16(34735) diff --git a/src/dataset.jl b/src/dataset.jl index d6ef443..208eed7 100644 --- a/src/dataset.jl +++ b/src/dataset.jl @@ -24,11 +24,15 @@ mutable struct LasDataset """Extra user bytes packed between the Header block and the first VLR of the source LAS file""" const 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} + function LasDataset(header::LasHeader, pointcloud::Table, vlrs::Vector{<:LasVariableLengthRecord}, evlrs::Vector{<:LasVariableLengthRecord}, - user_defined_bytes::Vector{UInt8}) + user_defined_bytes::Vector{UInt8}, + unit_conversion::SVector{3, Float64} = NO_CONVERSION) # do a few checks to make sure everything is consistent between the header and other data point_format_from_table = get_point_format(pointcloud) @@ -50,6 +54,9 @@ mutable struct LasDataset num_evlrs_in_header = number_of_evlrs(header) @assert num_evlrs == num_evlrs_in_header "Number of EVLR's in header $(num_evlrs_in_header) doesn't match number of EVLR's supplied $(num_evlrs)" + # make sure our unit conversions are non-zero to avoid NaN's + @assert all(unit_conversion .> 0) "Unit conversion factors must be non-zero! Got $(unit_conversion)" + # 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))) @@ -114,7 +121,7 @@ mutable struct LasDataset end end end - return new(header, las_pc, user_pc, Vector{LasVariableLengthRecord}(vlrs), Vector{LasVariableLengthRecord}(evlrs), user_defined_bytes) + return new(header, las_pc, user_pc, Vector{LasVariableLengthRecord}(vlrs), Vector{LasVariableLengthRecord}(evlrs), user_defined_bytes, unit_conversion) end end @@ -160,6 +167,13 @@ Extract the set of user-defined bytes from a `LasDataset` `las` """ get_user_defined_bytes(las::LasDataset) = las.user_defined_bytes +""" + $(TYPEDSIGNATURES) + +Get the unit factor conversion that was applied to this dataset when ingested +""" +get_unit_conversion(las::LasDataset) = las.unit_conversion + function Base.show(io::IO, las::LasDataset) println(io, "LAS Dataset") println(io, "\tNum Points: $(length(get_pointcloud(las)))") diff --git a/src/read.jl b/src/read.jl index 8ef4056..0f17057 100644 --- a/src/read.jl +++ b/src/read.jl @@ -119,13 +119,15 @@ function read_las_data(io::TIO, required_columns::TTuple=DEFAULT_LAS_COLUMNS; as_table = make_table(records, required_columns, xyz) - if convert_to_metres + conversion = if convert_to_metres convert_units!(as_table, vlrs, convert_x_y_units, convert_z_units) + else + NO_CONVERSION end evlrs = Vector{LasVariableLengthRecord}(map(_ -> read(io, LasVariableLengthRecord, true), 1:number_of_evlrs(header))) - LasDataset(header, as_table, vlrs, evlrs, user_defined_bytes) + LasDataset(header, as_table, vlrs, evlrs, user_defined_bytes, conversion) end """ @@ -261,7 +263,7 @@ function convert_units!(pointcloud::AbstractVector{<:NamedTuple}, vlrs::Vector{L these_are_wkts = is_ogc_wkt_record.(vlrs) # we are not requesting unit conversion and there is no OGC WKT VLR if ismissing(convert_x_y_units) && ismissing(convert_z_units) && count(these_are_wkts) == 0 - return nothing + return NO_CONVERSION else @assert count(these_are_wkts) == 1 "Expected to find 1 OGC WKT VLR, instead found $(count(these_are_wkts))" ogc_wkt = get_data(vlrs[findfirst(these_are_wkts)]) @@ -269,9 +271,13 @@ function convert_units!(pointcloud::AbstractVector{<:NamedTuple}, vlrs::Vector{L if !ismissing(conversion) && any(conversion .!= 1.0) @info "Positions converted to meters using conversion $(conversion)" pointcloud = pointcloud.position .= map(p -> p .* conversion, pointcloud.position) + return conversion + else + return NO_CONVERSION end end end + return NO_CONVERSION end """ diff --git a/src/write.jl b/src/write.jl index 89dc59e..91991be 100644 --- a/src/write.jl +++ b/src/write.jl @@ -12,8 +12,6 @@ Saves a pointcloud to LAS or LAZ. The appropriate LAS version and point format i * `evlrs` : Collection of Extended Variable Length Records to write to the LAS file, default `LasVariableLengthRecord[]` * `user_defined_bytes` : Any user-defined bytes to write in between the VLRs and point records, default `UInt8[]` * `scale` : Scaling factor applied to points on writing, default `LAS.POINT_SCALE` -* `compressed` : Whether or not data is to be written to a compressed .laz file, default false - --- $(METHODLIST) """ @@ -25,14 +23,14 @@ function save_las(file_name::AbstractString, pointcloud::AbstractVector{<:NamedT kwargs...) open_func = get_open_func(file_name) open_func(file_name, "w") do io - write_las(io, pointcloud, vlrs, evlrs, user_defined_bytes, scale, is_laz(file_name)) + write_las(io, pointcloud, vlrs, evlrs, user_defined_bytes, scale) end end function save_las(file_name::AbstractString, las::LasDataset) open_func = get_open_func(file_name) open_func(file_name, "w") do io - write_las(io, las, is_laz(file_name)) + write_las(io, las) end end @@ -45,7 +43,7 @@ function save_las(file_name::AbstractString, scale::Real = POINT_SCALE) where {TRecord <: LasRecord} open_func = get_open_func(file_name) open_func(file_name, "w") do io - write_las(io, header, point_records, vlrs, evlrs, user_defined_bytes, scale, is_laz(file_name)) + write_las(io, header, point_records, vlrs, evlrs, user_defined_bytes, scale) end end @@ -53,10 +51,9 @@ function write_las(io::IO, pointcloud::AbstractVector{<:NamedTuple}, vlrs::Vector{<:LasVariableLengthRecord} = LasVariableLengthRecord[], evlrs::Vector{<:LasVariableLengthRecord} = LasVariableLengthRecord[], user_defined_bytes::Vector{UInt8} = UInt8[], - scale::Real = POINT_SCALE, - compressed::Bool = false) + scale::Real = POINT_SCALE) point_format = get_point_format(pointcloud) - write_las(io, pointcloud, point_format, vlrs, evlrs, user_defined_bytes, scale, compressed) + write_las(io, pointcloud, point_format, vlrs, evlrs, user_defined_bytes, scale) end @@ -72,26 +69,30 @@ Write a pointcloud and additional VLR's and user-defined bytes to an IO stream i * `evlrs` : Collection of Extended Variable Length Records to write to `io` * `user_defined_bytes` : Any user-defined bytes to write in between the VLRs and point records * `scale` : Scaling factor applied to points on writing -* `compressed` : Whether or not data is to be written to a compressed .laz file, default false """ function write_las(io::IO, pointcloud::AbstractVector{<:NamedTuple}, point_format::Type{TPoint}, vlrs::Vector{<:LasVariableLengthRecord}, evlrs::Vector{<:LasVariableLengthRecord}, user_defined_bytes::Vector{UInt8}, - scale::Real, - compressed::Bool = false) where {TPoint} + scale::Real) where {TPoint} # automatically construct a header that's consistent with the data and point format we've supplied header = make_consistent_header(pointcloud, point_format, vlrs, evlrs, user_defined_bytes, scale) - write_las(io, LasDataset(header, pointcloud, vlrs, evlrs, user_defined_bytes), compressed) + write_las(io, LasDataset(header, pointcloud, vlrs, evlrs, user_defined_bytes)) end -function write_las(io::IO, las::LasDataset, compressed::Bool = false) +function write_las(io::IO, las::LasDataset) header = get_header(las) vlrs = get_vlrs(las) pc = get_pointcloud(las) + # reverse our unit conversion done on ingest so the positions match the header / VLR info + if get_unit_conversion(las) != NO_CONVERSION + pc.position .= map(p -> p ./ get_unit_conversion(las), pc.position) + end + + this_point_format = point_format(header) xyz = spatial_info(header) diff --git a/test/file_io.jl b/test/file_io.jl index 2e3c56d..d54d98c 100644 --- a/test/file_io.jl +++ b/test/file_io.jl @@ -382,7 +382,6 @@ end ) tmpdir = mktempdir() - # tmpdir = "/home/msb/git/FugroLAS.jl/test_folder" output_file_path = joinpath(tmpdir, "test_las_write_laspoint6.las") wkt_str = """GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]""" @@ -514,4 +513,17 @@ end end end end +end + +@testset "Unit Conversion" begin + ogc_ft_file = joinpath(@__DIR__, "test_files/ogc_ft.las") + las = load_las(ogc_ft_file; convert_to_metres = true) + # make sure we get the write unit conversion - should be all feet + @test get_unit_conversion(las) == SVector{3, Float64}(0.304800609601219, 0.304800609601219, 0.304800609601219) + mktempdir() do tmp + # make sure that when we read a LAS file, do a unit conversion then write it again, we get the same dataset back + out_file = joinpath(tmp, "pc.las") + save_las(out_file, las) + @test load_las(ogc_ft_file, [:position]; convert_to_metres = false) == load_las(out_file, [:position]; convert_to_metres = false) + end end \ No newline at end of file diff --git a/test/test_files/ogc_ft.las b/test/test_files/ogc_ft.las new file mode 100644 index 0000000..374c9dd Binary files /dev/null and b/test/test_files/ogc_ft.las differ