Merge pull request #31 from fugro-oss/ReverseUnitConv

Reverse coordinate system conversion on write

Project.toml

@@ -1,7 +1,7 @@
 name = "LAS"
 uuid = "cc498e2a-d443-4943-8f26-2a8a0f3c7cdb"
 authors = ["BenCurran98 <b.curran@fugro.com>"]
-version = "0.2.0"
+version = "0.2.1"
 
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"

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

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)

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)))")

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,17 +263,21 @@ 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)])
             conversion = conversion_from_vlrs(ogc_wkt, convert_x_y_units = convert_x_y_units, convert_z_units = convert_z_units)
             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
 
 """

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,18 +43,17 @@ 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
 
 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)
 

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