Las1.4 support with all point formats

.github/workflows/CompatHelper.yml

@@ -21,4 +21,4 @@ jobs:
       - name: CompatHelper.main()
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
-        run: julia -e 'using CompatHelper; CompatHelper.main()'
+        run: julia -e 'using CompatHelper; CompatHelper.main()'

.github/workflows/TagBot.yml

@@ -8,4 +8,4 @@ jobs:
     steps:
       - uses: JuliaRegistries/TagBot@v1
         with:
-          token: ${{ secrets.GITHUB_TOKEN }}
+          token: ${{ secrets.GITHUB_TOKEN }}

.gitignore

@@ -1,3 +1,8 @@
+.vscode
+
 *.jl.cov
 *.jl.*.cov
 *.jl.mem
+
+Manifest.toml
+deps/build.log

CHANGELOG.md

@@ -1,6 +1,9 @@
 # Changelog
 
 ## [Unreleased]
+## [0.4.0] - 2018-10-15
+### Changed
+- Support all LAS/LAZ versions and point formats (without waveform data)
 
 ## [0.3.0] - 2018-10-15
 ### Changed

Project.toml

@@ -3,7 +3,7 @@ uuid = "570499db-eae3-5eb6-bdd5-a5326f375e68"
 keywords = ["LAS", "lidar", "IO"]
 license = "MIT"
 desc = "Read and write LAS lidar point cloud files"
-version = "0.3.6"
+version = "0.4.0"
 
 [deps]
 ColorTypes = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
@@ -12,6 +12,7 @@ FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 FixedPointNumbers = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
 GeometryTypes = "4d00f742-c7ba-57c2-abde-4428a4b178cb"
 Mmap = "a63ad114-7e13-5084-954f-fe012c677804"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
 [compat]

README.md

@@ -5,7 +5,9 @@
 
 Julia package for reading and writing the LAS lidar format.
 
-This is a pure Julia package for reading and writing ASPRS `.las` files. Currently only LAS versions 1.1 - 1.3 and point formats 0 - 3 are supported. For LAZ support see below.
+This is a pure Julia package for reading and writing ASPRS `.las` files. Currently all LAS versions 1.1 - 1.4 and point formats 0 - 10 are semi-supported. By semi-supported, we mean that we do not read or write the waveform data.
+
+TODO - Support for Waveform data is future work.
 
 If the file fits into memory, it can be loaded using
 
@@ -28,9 +30,8 @@ where `points` is now a memory mapped `PointVector{LasPoint3}` which behaves in
 See `test/runtests.jl` for other usages.
 
 ## LAZ support
-We advise to use [LazIO](https://github.com/evetion/LazIO.jl), which works out of the box and is compatible with LasIO.
-
-The compressed LAZ format is supported by LasIO itself, but requires the user to make sure the `laszip` executable can be found in the PATH. LAZ files are piped through `laszip` to provide reading and writing capability. `laszip` is not distributed with this package. One way to get it is to download `LAStools` from https://rapidlasso.com/. The LAStools ZIP file already contains `laszip.exe` for Windows, for Linux or Mac it needs to be compiled first. When this is done this should work just like with LAS:
+LasIO comes with laszip which will be used to read/write laz files just like LAS file. There is no need for LazIO anymore.
+TODO - build the `laszip` instead of having the executable sitting there.
 
 ```julia
 using FileIO, LasIO

deps/build.jl

@@ -0,0 +1,25 @@
+@info "Downloading LasTools"
+import Pkg
+if Sys.isapple()
+    Pkg.add("LibGit2")
+    using LibGit2
+    resource_path = joinpath(dirname(@__DIR__), "resources")
+    lastools_path = joinpath(dirname(@__DIR__), "LAStools")
+    lastools_build_path = joinpath(lastools_path, "build")
+    lastools_install_path = joinpath(lastools_build_path, "install")
+    lastools_executable_path = joinpath(lastools_install_path, "bin")
+    if !isdir(lastools_path)
+        LibGit2.clone("https://github.com/LAStools/LAStools", lastools_path)
+    end
+
+    mkpath(lastools_build_path)
+    cd(lastools_build_path)
+    run(`cmake -DCMAKE_INSTALL_PREFIX=$(lastools_install_path) ../`)
+    run(`cmake --build . --target install --config Release`)
+
+    # find the las executable
+    laszip_executables = filter(x -> startswith(x, "laszip"), readdir(lastools_executable_path))
+    length(laszip_executables) == 0 && error("Unable to build a laszip executable for $(Sys.MACHINE)") 
+    cp(joinpath(lastools_executable_path, laszip_executables[1]), joinpath(resource_path, "laszip"), force=true)
+    rm(lastools_path, recursive=true)
+end

src/LasIO.jl

@@ -13,10 +13,19 @@ export
     LasHeader,
     LasVariableLengthRecord,
     LasPoint,
+    LasPoint_0_5,
+    LasPoint_6_10,
     LasPoint0,
     LasPoint1,
     LasPoint2,
     LasPoint3,
+    LasPoint4,
+    LasPoint5,
+    LasPoint6,
+    LasPoint7,
+    LasPoint8,
+    LasPoint9,
+    LasPoint10,
     PointVector,
 
     # Functions on LasHeader
@@ -27,12 +36,14 @@ export
     # Functions on LasPoint
     return_number,
     number_of_returns,
+    scanner_channel,
     scan_direction,
     edge_of_flight_line,
     classification,
     synthetic,
     key_point,
     withheld,
+    overlap,
     xcoord,
     ycoord,
     zcoord,
@@ -43,6 +54,8 @@ export
     gps_time,
     raw_classification,
     flag_byte,
+    flag_byte_1,
+    flag_byte_2,
 
     # extended from ColorTypes
     red,
@@ -63,6 +76,8 @@ function __init__()
     # https://github.com/JuliaIO/FileIO.jl/blob/master/src/registry.jl
     add_format(format"LAS", "LASF", ".las", [:LasIO])
     add_format(format"LAZ", (), ".laz", [:LasIO])
+    add_loader(format"LAZ", :LasIO)
+    add_loader(format"LAZ", :LasIO)
 end
 
 end # module

src/fileio.jl

@@ -1,5 +1,17 @@
 using Mmap
 
+function get_laszip_executable_path()
+    if Sys.iswindows()
+        return joinpath(dirname(@__DIR__), "resources", "laszip.exe")
+    else # point to the linux build of the executable
+        return joinpath(dirname(@__DIR__), "resources", "laszip")
+    end
+end
+
+function get_record_count(header::LasHeader)
+    return header.extended_number_of_point_records > 0 ? Int(header.extended_number_of_point_records) : Int(header.records_count)
+end
+
 function pointformat(header::LasHeader)
     id = header.data_format_id
     if id == 0x00
@@ -10,6 +22,20 @@ function pointformat(header::LasHeader)
         return LasPoint2
     elseif id == 0x03
         return LasPoint3
+    elseif id == 0x04
+        return LasPoint4
+    elseif id == 0x05
+        return LasPoint5
+    elseif id == 0x06
+        return LasPoint6
+    elseif id == 0x07
+        return LasPoint7
+    elseif id == 0x08
+        return LasPoint8
+    elseif id == 0x09
+        return LasPoint9
+    elseif id == 0x0a
+        return LasPoint10
     else
         error("unsupported point format $(Int(id))")
     end
@@ -28,11 +54,13 @@ end
 function load(s::Base.AbstractPipe)
     skiplasf(s)
     header = read(s, LasHeader)
-
-    n = header.records_count
+    n = get_record_count(header)
     pointtype = pointformat(header)
+
+    @info "Reading $(n) '$(pointtype)' points"
     pointdata = Vector{pointtype}(undef, n)
     for i=1:n
+        i%1000000 == 0 && @info("Read $(i)/$(n) points")
         pointdata[i] = read(s, pointtype)
     end
     header, pointdata
@@ -41,17 +69,18 @@ end
 function load(s::Stream{format"LAS"}; mmap=false)
     skiplasf(s)
     header = read(s, LasHeader)
-
-    n = header.records_count
+    n = get_record_count(header)
     pointtype = pointformat(header)
 
+    @info "Reading $(n) '$(pointtype)' points"
     if mmap
         pointsize = Int(header.data_record_length)
         pointbytes = Mmap.mmap(s.io, Vector{UInt8}, n*pointsize, position(s))
         pointdata = PointVector{pointtype}(pointbytes, pointsize)
     else
         pointdata = Vector{pointtype}(undef, n)
         for i=1:n
+            i%1000000 == 0 && @info("Read $(i)/$(n) points")
             pointdata[i] = read(s, pointtype)
         end
     end
@@ -61,8 +90,10 @@ end
 
 function load(f::File{format"LAZ"})
     # read las from laszip, which decompresses to stdout
-    open(`laszip -olas -stdout -i $(filename(f))`) do s
-        load(s)
+    open(`$(get_laszip_executable_path()) -olas -stdout -i $(filename(f))`) do s
+        h,p = load(s)
+        read(s)
+        return h, p
     end
 end
 
@@ -86,9 +117,11 @@ end
 
 function save(s::Stream{format"LAS"}, header::LasHeader, pointdata::AbstractVector{<:LasPoint})
     # checks
-    header_n = header.records_count
+    header_n = get_record_count(header)
     n = length(pointdata)
-    msg = "number of records in header ($header_n) does not match data length ($n)"
+    msg = "Number of records in header ($header_n) does not match data length ($n)"
+    @info "Writing $(n) '$(typeof(pointdata[1]))' points"
+
     @assert header_n == n msg
 
     # write header
@@ -103,7 +136,7 @@ end
 
 function save(f::File{format"LAZ"}, header::LasHeader, pointdata::AbstractVector{<:LasPoint})
     # pipes las to laszip to write laz
-    open(`laszip -olaz -stdin -o $(filename(f))`, "w") do s
+    open(`$(get_laszip_executable_path()) -olaz -stdin -o $(filename(f))`, "w") do s
         savebuf(s, header, pointdata)
     end
 end
@@ -113,10 +146,11 @@ end
 # but it speeds up a lot when the result is piped to laszip.
 function savebuf(s::IO, header::LasHeader, pointdata::AbstractVector{<:LasPoint})
     # checks
-    header_n = header.records_count
+    header_n = get_record_count(header)
     n = length(pointdata)
     msg = "number of records in header ($header_n) does not match data length ($n)"
     @assert header_n == n msg
+    @info "Writing $(n) '$(typeof(pointdata[1]))' points to LAZ file"
 
     # write header
     write(s, magic(format"LAS"))
@@ -125,7 +159,7 @@ function savebuf(s::IO, header::LasHeader, pointdata::AbstractVector{<:LasPoint}
     # 2048 points seemed to be an optimum for the libLAS_1.2.las testfile
     npoints_buffered = 2048
     bufsize = header.data_record_length * npoints_buffered
-    buf = IOBuffer(bufsize)
+    buf = IOBuffer(sizehint=bufsize)
     # write points
     for (i, p) in enumerate(pointdata)
         write(buf, p)

src/header.jl

@@ -47,6 +47,11 @@ mutable struct LasHeader
     y_min::Float64
     z_max::Float64
     z_min::Float64
+    start_of_waveform_data_packet_record::UInt64
+    start_of_first_extended_variable_length_record::UInt64
+    number_of_extended_variable_length_records::UInt32
+    extended_number_of_point_records::UInt64
+    extended_number_of_points_by_return::Vector{UInt64}
     variable_length_records::Vector{LasVariableLengthRecord}
     user_defined_bytes::Vector{UInt8}
 end
@@ -61,8 +66,8 @@ function Base.getproperty(h::LasHeader, s::Symbol)
 end
 
 function Base.show(io::IO, header::LasHeader)
-    n = Int(header.records_count)
-    println(io, "LasHeader with $n points.")
+    n = header.extended_number_of_point_records > 0 ? Int(header.extended_number_of_point_records) : Int(header.records_count)
+    println(io, "LasHeader with $(n) points.")
     println(io, string("\tfile_source_id = ", header.file_source_id))
     println(io, string("\tglobal_encoding = ", header.global_encoding))
     println(io, string("\tguid_1 = ", header.guid_1))
@@ -81,7 +86,7 @@ function Base.show(io::IO, header::LasHeader)
     println(io, string("\tdata_format_id = ", header.data_format_id))
     println(io, string("\tdata_record_length = ", header.data_record_length))
     println(io, string("\trecords_count = ", header.records_count))
-    println(io, string("\tpoint_return_count = ", header.point_return_count))
+    println(io, string("\tpoint_return_count = ", [Int(i) for i in header.point_return_count]))
     println(io, string("\tx_scale = ", header.x_scale))
     println(io, string("\ty_scale = ", header.y_scale))
     println(io, string("\tz_scale = ", header.z_scale))
@@ -95,10 +100,15 @@ function Base.show(io::IO, header::LasHeader)
     println(io, @sprintf "\tz_max = %.7f" header.z_max)
     println(io, @sprintf "\tz_min = %.7f" header.z_min)
 
+    if header.version_minor > 3
+        println(io, string("\tnumber_of_extended_variable_length_records = ", header.number_of_extended_variable_length_records))
+        println(io, string("\textended_number_of_point_records = ", header.extended_number_of_point_records))
+        println(io, string("\textended_number_of_points_by_return = ", [Int(i) for i in header.extended_number_of_points_by_return]))
+    end
     if !isempty(header.variable_length_records)
         nrecords = min(10, size(header.variable_length_records, 1))
 
-        println(io, string("\tvariable_length_records (max 10) = "))
+        println(io, string("\tvariable_length_records (showing first $(nrecords)/$(size(header.variable_length_records, 1))) = "))
         for vlr in header.variable_length_records[1:nrecords]
             println(io, "\t\t($(vlr.user_id), $(vlr.record_id)) => ($(vlr.description), $(sizeof(vlr.data)) bytes...)")
         end
@@ -163,10 +173,21 @@ function Base.read(io::IO, ::Type{LasHeader})
     z_max = read(io, Float64)
     z_min = read(io, Float64)
     lasversion = VersionNumber(version_major, version_minor)
+    start_of_waveform_data_packet_record = UInt64(0)
+    start_of_first_extended_variable_length_record = UInt64(0)
+    number_of_extended_variable_length_records = UInt32(0)
+    extended_number_of_point_records = UInt64(0)
+    extended_number_of_points_by_return = fill(UInt32(0), 15)
+
     if lasversion >= v"1.3"
-         # start of waveform data record (unsupported)
-        _ = read(io, UInt64)
+        # start of waveform data record (unsupported)
+        start_of_waveform_data_packet_record = read(io, UInt64)
+        start_of_first_extended_variable_length_record = read(io, UInt64)
+        number_of_extended_variable_length_records = read(io, UInt32)
+        extended_number_of_point_records = read(io, UInt64) # this is a UInt64 which is different from records_count which is UInt32
+        extended_number_of_points_by_return = read!(io, Vector{UInt64}(undef, 15))
     end
+
     vlrs = [read(io, LasVariableLengthRecord, false) for i=1:n_vlr]
 
     # From here until the data_offset everything is read in
@@ -209,8 +230,13 @@ function Base.read(io::IO, ::Type{LasHeader})
         y_min,
         z_max,
         z_min,
+        start_of_waveform_data_packet_record,
+        start_of_first_extended_variable_length_record,
+        number_of_extended_variable_length_records,
+        extended_number_of_point_records,
+        extended_number_of_points_by_return,
         vlrs,
-        user_defined_bytes
+        user_defined_bytes,
     )
 end
 
@@ -251,7 +277,11 @@ function Base.write(io::IO, h::LasHeader)
     lasversion = VersionNumber(h.version_major, h.version_minor)
     if lasversion >= v"1.3"
         # start of waveform data record (unsupported)
-        write(io, UInt64(0))
+        write(io, h.start_of_waveform_data_packet_record) # start_of_waveform_data_packet_record
+        write(io, h.start_of_first_extended_variable_length_record) # start_of_first_extended_variable_length_record
+        write(io, h.number_of_extended_variable_length_records) # number_of_extended_variable_length_records
+        write(io, h.extended_number_of_point_records)  # extended_number_of_point_records
+        write(io, h.extended_number_of_points_by_return) # extended_number_of_points_by_return
     end
     for i in 1:h.n_vlr
         write(io, h.variable_length_records[i])