[WIP] First DataStreams version

README.md

@@ -7,7 +7,53 @@ Julia package for reading and writing the LAS lidar format.
 
 This is a pure Julia alternative to [LibLAS.jl](https://github.com/visr/LibLAS.jl) or [Laszip.jl](https://github.com/joa-quim/Laszip.jl). Currently only LAS versions 1.1 - 1.3 and point formats 0 - 3 are supported. For LAZ support see below.
 
-If the file fits into memory, it can be loaded using
+## Usage
+
+```julia
+julia> using LasIO
+julia> s = LasIO.Source("test/srs.las")
+LasIO.Source(Data.Schema:
+rows: 10  cols: 10
+Columns:
+ "x"                   Int32  
+ "y"                   Int32  
+ "z"                   Int32  
+ "intensity"           UInt16 
+ "flag_byte"           UInt8  
+ "raw_classification"  UInt8  
+ "scan_angle"          Int8   
+ "user_data"           UInt8  
+ "pt_src_id"           UInt16 
+ "gps_time"            Float64, LasHeader with 10 points.
+, IOStream(<file test/srs.las>), "test/srs.las", 759)
+
+julia> using DataStreams, DataFrames
+julia> d = Data.stream!(s, DataFrame);
+julia> Data.close!(d)
+10×10 DataFrames.DataFrame
+│ Row │ x         │ y         │ z      │ intensity │ flag_byte │ raw_classification │ scan_angle │ user_data │ pt_src_id │ gps_time  │
+├─────┼───────────┼───────────┼────────┼───────────┼───────────┼────────────────────┼────────────┼───────────┼───────────┼───────────┤
+│ 1   │ 2.89814e5 │ 4.32098e6 │ 170.76 │ 0x0104    │ 0x30      │ 0x02               │ 0          │ 0x00      │ 0x0000    │ 4.99451e5 │
+│ 2   │ 2.89815e5 │ 4.32098e6 │ 170.76 │ 0x0118    │ 0x30      │ 0x02               │ 0          │ 0x00      │ 0x0000    │ 4.99451e5 │
+│ 3   │ 2.89815e5 │ 4.32098e6 │ 170.75 │ 0x0118    │ 0x30      │ 0x02               │ 0          │ 0x00      │ 0x0000    │ 4.99451e5 │
+│ 4   │ 2.89816e5 │ 4.32098e6 │ 170.74 │ 0x0118    │ 0x30      │ 0x02               │ 0          │ 0x00      │ 0x0000    │ 4.99451e5 │
+│ 5   │ 2.89816e5 │ 4.32098e6 │ 170.68 │ 0x0104    │ 0x30      │ 0x02               │ 0          │ 0x00      │ 0x0000    │ 4.99451e5 │
+│ 6   │ 2.89817e5 │ 4.32098e6 │ 170.66 │ 0x00f0    │ 0x30      │ 0x02               │ 0          │ 0x00      │ 0x0000    │ 4.99451e5 │
+│ 7   │ 289817.0  │ 4.32098e6 │ 170.63 │ 0x00f0    │ 0x30      │ 0x02               │ 0          │ 0x00      │ 0x0000    │ 4.99451e5 │
+│ 8   │ 2.89818e5 │ 4.32098e6 │ 170.62 │ 0x0118    │ 0x30      │ 0x02               │ 0          │ 0x00      │ 0x0000    │ 4.99451e5 │
+│ 9   │ 289818.0  │ 4.32098e6 │ 170.61 │ 0x0118    │ 0x30      │ 0x02               │ 0          │ 0x00      │ 0x0000    │ 4.99451e5 │
+│ 10  │ 2.89819e5 │ 4.32098e6 │ 170.58 │ 0x0104    │ 0x30      │ 0x02               │ 0          │ 0x00      │ 0x0000    │ 4.99451e5 │
+
+julia> Data.reset!(s)
+julia> d = Data.stream!(s, LasIO.Sink, "test_final.las");
+julia> Data.close!(d)
+LasIO.Sink{LasIO.LasPoint1}(IOStream(<file test_final.las>), LasHeader with 10 points.
+, LasIO.LasPoint1)
+```
+
+### Legacy API
+
+For backwards compatibility, the legacy API is still provided.
 
 ```julia
 using FileIO, LasIO

REQUIRE

@@ -3,3 +3,5 @@ ColorTypes 0.3.3
 FixedPointNumbers 0.3.2
 FileIO 0.3.0
 GeometryTypes 0.4.0
+DataStreams
+Parameters

src/LasIO.jl

@@ -2,10 +2,12 @@ __precompile__()
 
 module LasIO
 
+using DataStreams
 using FileIO
 using FixedPointNumbers
 using ColorTypes
 using GeometryTypes # for conversion
+using Parameters
 
 export
     # Types
@@ -19,6 +21,7 @@ export
 
     # Functions on LasHeader
     update!,
+    pointtype,
 
     # Functions on LasPoint
     return_number,
@@ -52,6 +55,7 @@ include("point.jl")
 include("util.jl")
 include("fileio.jl")
 include("srs.jl")
+include("datastreams.jl")
 
 function __init__()
     # these should eventually go in

src/datastreams.jl

@@ -0,0 +1,189 @@
+########
+# SOURCE
+########
+
+const minimum_coordinate = typemin(Int32)
+const maximum_coordinate = typemax(Int32)
+
+mutable struct Source <: Data.Source
+    schema::Data.Schema
+    header::LasHeader
+    io::IO
+    fullpath::String
+    datapos::Int
+end
+
+dict_of_struct(T) = Dict((String(fieldname(typeof(T), i)), getfield(T, i)) for i = 1:nfields(T))
+
+function Source(f::AbstractString)
+    !isfile(f) && error("Please provide valid path.")
+
+    io = is_windows() ? open(f) : IOBuffer(Mmap.mmap(f))
+
+    skiplasf(io)
+    header = read(io, LasHeader)
+
+    n = header.records_count
+    pointtype = pointformat(header)
+
+    # Schema
+    columns = Vector{String}(fieldnames(pointtype))
+    types = Vector{Type}([fieldtype(pointtype, i) for i in 1:nfields(pointtype)])
+    types[1:3] = [Float64, Float64, Float64]  # Convert XYZ coordinates
+    sch = Data.Schema(types, columns, n, dict_of_struct(header))
+
+    return Source(sch, header, io, f, position(io))
+end
+
+Data.reset!(s::LasIO.Source) = (seek(s.io, s.datapos); return nothing)
+Data.schema(source::LasIO.Source) = source.schema
+Data.accesspattern(::Type{LasIO.Source}) = Data.Sequential
+Data.isdone(source::LasIO.Source, row::Int, col::Int) = eof(source.io) || (row, col) > Data.size(Data.schema(source))
+Data.isdone(source::LasIO.Source, row, col, rows, cols) = eof(source.io) || row > rows || col > cols
+Data.streamtype(::Type{<:LasIO.Source}, ::Type{Data.Field}) = true
+
+# Data.streamfrom(source::LasIO.Source, st::Type{Data.Row}, t::Type{T}, row::Int) where {T} = read(source.io, pointformat(source.header))
+function Data.streamfrom(source::LasIO.Source, ::Type{Data.Field}, ::Type{T}, row::Int, col::Int) where {T}
+
+    # XYZ => scale stored Integers to Floats using header information
+    if col == 1
+        return muladd(read(source.io, Int32), source.header.x_scale, source.header.x_offset)
+    elseif col == 2
+        return muladd(read(source.io, Int32), source.header.y_scale, source.header.y_offset)
+    elseif col == 3
+        return muladd(read(source.io, Int32), source.header.z_scale, source.header.z_offset)
+
+    # Otherwise return read value
+    else
+        v = read(source.io, T)
+    end
+end
+
+######
+# SINK
+######
+
+mutable struct Sink{T} <: Data.Sink where T <: LasPoint
+    stream::IO
+    header::LasHeader
+    pointformat::Type{T}
+    bbox::Vector{Float64}
+    returncount::Vector{UInt32}
+end
+
+# setup header and empty pointvector
+function Sink(sch::Data.Schema, S::Type{Data.Field}, append::Bool, fullpath::AbstractString, bbox::Vector{Float64}; scale::Real=0.01, epsg=nothing)
+
+    # validate input
+    length(bbox) != 6 && error("Provide bounding box as (xmin, ymin, zmin, xmax, ymax, zmax)")
+
+    s = open(fullpath, "w")
+
+    # determine point version and derivatives
+    pointtype = gettype(Data.header(sch))
+    data_format_id = pointformat(pointtype)
+    data_record_length = packed_size(pointtype)
+    n = Data.size(sch, 1)  # rows
+
+    # setup and validate scaling based on bbox and scale
+    x_offset = determine_offset(bbox[1], bbox[4], scale)
+    y_offset = determine_offset(bbox[2], bbox[5], scale)
+    z_offset = determine_offset(bbox[3], bbox[6], scale)
+
+    # create header
+    header = LasHeader(data_format_id=data_format_id, data_record_length=data_record_length, records_count=n,
+        x_max=bbox[4], x_min=bbox[1], y_max=bbox[5], y_min=bbox[2], z_max=bbox[6], z_min=bbox[3],
+        x_scale=scale, y_scale=scale, z_scale=scale,
+        x_offset=x_offset, y_offset=y_offset, z_offset=z_offset)
+    epsg != nothing && epsg_code!(header, epsg)
+
+    # write header
+    write(s, magic(format"LAS"))
+    write(s, header)
+
+    # empty return count
+    return_count = Array{UInt32}([0, 0, 0, 0, 0])
+
+    # return stream, position is correct for writing points
+    return Sink(s, header, pointtype, bbox, return_count)
+end
+
+# Update existing Sink
+# function Sink(sink::LasIO.Sink, sch::Data.Schema, S::Type{StreamType}, append::Bool; reference::Vector{UInt8}=UInt8[])
+    # return Sink
+# end
+
+"""Determine LAS versions based on specific columns."""
+function gettype(columns::Vector{String})
+    # LAS versions only differ in gps_time and rgb color information
+    has_gps, has_color = false, false
+    ("gps_time" in columns) && (has_gps = true)
+    ("red" in columns || "green" in columns || "blue" in columns) && (has_color = true)
+
+    has_gps && has_color && return LasPoint3
+    has_color && return LasPoint2
+    has_gps && return LasPoint1
+    return LasPoint0
+end
+
+Data.streamtypes(::Type{LasIO.Sink}) = [Data.Field, Data.Row]
+Data.cleanup!(sink::LasIO.Sink) = nothing
+Data.weakrefstrings(::Type{LasIO.Sink}) = false
+
+function update_sink(sink::LasIO.Sink, col::Integer, val::T) where {T}
+    if col == 1
+        (val > sink.bbox[1] || sink.bbox[1] == 0.0) && (setindex!(sink.bbox, val, 1))  # xmax
+        (val < sink.bbox[2] || sink.bbox[2] == 0.0) && (setindex!(sink.bbox, val, 2))  # xmin
+    end
+    if col == 2
+        (val > sink.bbox[3] || sink.bbox[3] == 0.0) && (setindex!(sink.bbox, val, 3))  # ymax
+        (val < sink.bbox[4] || sink.bbox[4] == 0.0) && (setindex!(sink.bbox, val, 4))  # ymin
+    end
+    if col == 3
+        (val > sink.bbox[5] || sink.bbox[5] == 0.0) && (setindex!(sink.bbox, val, 5))  # zmax
+        (val < sink.bbox[6] || sink.bbox[6] == 0.0) && (setindex!(sink.bbox, val, 6))  # zmin  
+    end
+    if col == 5
+        return_number = val & 0b00000111
+        return_number < 5 && (sink.returncount[return_number+1] += 1)
+    end
+end
+
+# actually write points to our pointvector
+function Data.streamto!(sink::LasIO.Sink, S::Type{Data.Field}, val, row, col)
+    # TODO(evetion) check if stream position is at row*col
+    update_sink(sink, col, val)
+
+    # XYZ => scale given Floats to Integers using header information
+    if col == 1
+        write(sink.stream, round(Int32, ((val - sink.header.x_offset) / sink.header.x_scale)))
+    elseif col == 2
+        write(sink.stream, round(Int32, ((val - sink.header.y_offset) / sink.header.y_scale)))
+    elseif col == 3
+        write(sink.stream, round(Int32, ((val - sink.header.z_offset) / sink.header.z_scale)))
+    else
+        write(sink.stream, val)
+    end
+end
+
+# save file
+function Data.close!(sink::LasIO.Sink)
+
+    # update header
+    header = sink.header
+    header.x_max = sink.bbox[1]
+    header.x_min = sink.bbox[2]
+    header.y_max = sink.bbox[3]
+    header.y_min = sink.bbox[4]
+    header.z_max = sink.bbox[5]
+    header.z_min = sink.bbox[6]
+    header.point_return_count = sink.returncount
+
+    # seek back to beginning and write header
+    seekstart(sink.stream)
+    skiplasf(sink.stream)
+    write(sink.stream, header)
+
+    close(sink.stream)
+    return sink
+end

src/fileio.jl

@@ -14,8 +14,22 @@ function pointformat(header::LasHeader)
     end
 end
 
+function pointformat(t::Type{T}) where T <: LasPoint
+    if t == LasPoint0
+        return 0x00
+    elseif t == LasPoint1
+        return 0x01
+    elseif t == LasPoint2
+        return 0x02
+    elseif t == LasPoint3
+        return 0x03
+    else
+        error("unsupported point format $t")
+    end
+end
+
 # skip the LAS file's magic four bytes, "LASF"
-skiplasf(s::Union{Stream{format"LAS"}, Stream{format"LAZ"}, IO}) = read(s, UInt32)
+skiplasf(s::Union{Stream{format"LAS"}, Stream{format"LAZ"}, IO}) = skip(s, sizeof(UInt32))
 
 function load(f::File{format"LAS"})
     open(f) do s
@@ -45,8 +59,7 @@ end
 
 function read_header(f::AbstractString)
     open(f) do s
-        skiplasf(s)
-        read(s, LasHeader)
+        read_header(s::IO)
     end
 end
 

src/header.jl

@@ -14,40 +14,40 @@ Backward compatibility with LAS 1.1 – LAS 1.3 when payloads consist of only le
 content
 =#
 
-mutable struct LasHeader
-    file_source_id::UInt16
-    global_encoding::UInt16
-    guid_1::UInt32
-    guid_2::UInt16
-    guid_3::UInt16
-    guid_4::AbstractString
-    version_major::UInt8
-    version_minor::UInt8
-    system_id::AbstractString
-    software_id::AbstractString
-    creation_doy::UInt16
-    creation_year::UInt16
-    header_size::UInt16
-    data_offset::UInt32
-    n_vlr::UInt32
-    data_format_id::UInt8
-    data_record_length::UInt16
-    records_count::UInt32
-    point_return_count::Vector{UInt32}
-    x_scale::Float64
-    y_scale::Float64
-    z_scale::Float64
-    x_offset::Float64
-    y_offset::Float64
-    z_offset::Float64
-    x_max::Float64
-    x_min::Float64
-    y_max::Float64
-    y_min::Float64
-    z_max::Float64
-    z_min::Float64
-    variable_length_records::Vector{LasVariableLengthRecord}
-    user_defined_bytes::Vector{UInt8}
+@with_kw mutable struct LasHeader
+    file_source_id::UInt16 = UInt16(0)  # no known flightline
+    global_encoding::UInt16 = UInt16(0)
+    guid_1::UInt32 = UInt32(0)  # project IDs
+    guid_2::UInt16 = UInt16(0)
+    guid_3::UInt16 = UInt16(0)
+    guid_4::AbstractString = rpad("", 8)
+    version_major::UInt8 = UInt8(1)
+    version_minor::UInt8 = UInt8(0)
+    system_id::AbstractString = rpad("LasIO.jl datastream", 32)  # small char array would be nicer
+    software_id::AbstractString = rpad("LasIO.jl", 32)
+    creation_doy::UInt16 = UInt16(Dates.dayofyear(now()))
+    creation_year::UInt16 = UInt16(Dates.year(now()))
+    header_size::UInt16 = UInt16(227)
+    data_offset::UInt32 = UInt32(227)
+    n_vlr::UInt32 = UInt32(0)
+    data_format_id::UInt8 = UInt8(pointformat(LasPoint0))
+    data_record_length::UInt16 = UInt16(packed_size(LasPoint0))
+    records_count::UInt32 = UInt32(0)
+    point_return_count::Vector{UInt32} = Vector{UInt32}([0,0,0,0,0])
+    x_scale::Float64 = 0.01
+    y_scale::Float64 = 0.01
+    z_scale::Float64 = 0.01
+    x_offset::Float64 = 0.0
+    y_offset::Float64 = 0.0
+    z_offset::Float64 = 0.0
+    x_max::Float64 = 0.0
+    x_min::Float64 = 0.0
+    y_max::Float64 = 0.0
+    y_min::Float64 = 0.0
+    z_max::Float64 = 0.0
+    z_min::Float64 = 0.0
+    variable_length_records::Vector{LasVariableLengthRecord} = Vector{LasVariableLengthRecord}()
+    user_defined_bytes::Vector{UInt8} = Vector{UInt8}()
 end
 
 function Base.show(io::IO, header::LasHeader)