[R-sig-Geo] read function for LAS data

[Michael Sumner]

Hello,

I'm looking for interest in this functionality and eventually working
it into a package.

I don't actually use LAS data much, and only have one example file so
I'm hoping others who do or know others who would be interested can
help. I have previously posted this to r-spatial-devel.

R's support for the "raw" type provides a very neat way of reading
data from binary files, such as LAS files. These files have a header
containing metadata that is used to specify reading from the bulk data
in the file. Essentially, you can read large sets of records at once
into a raw matrix with readBin and then use the indexing tools to
extract the relevant bytes for each element of the record. This means
that the bulk read is fast, and processing is also vectorized as the
raw matrix can then be read in index pieces by readBin.

I've done this in a simplistic way in the attached functions

  o "publicHeaderDescription"  - function returns a list generated to hold
      information about the header

  o "readLAS" - read file function


It works for a LAS 1.0 file that I have (see lasinfo summary below)
and returns a matrix of record values.  Certainly it works fast and
efficiently enough for smallish files. I really don't have experience
with the likely range of sizes of these files.

The remaining work I see involves:

 - write as a "chunked" read so large files can be processed in parts,
    not as one big slurp
 - return only desired record values, and / or subset of records
 - extract all components of the records, and generalize for LAS 1.1, 1.2
 - re-organize the header read and arrangement (what I've done was
    easy enough for me to understand, but it's a bit of a mess)
 (- no doubt many more general improvements I haven't thought of)

To use:

source("http://staff.acecrc.org.au/~mdsumner/las/readLAS.R")

## I cannot provide LAS data, this is just an illustration

f <- "lfile.las"
## [1] 53.83393 Mb

file.info(f)$size/1e6

system.time({
  d <- readLAS(f)
})
# user  system elapsed
#   1.01    0.21    1.21

dim(d)
# [1] 1922632       5

colnames(d)
#[1] "x"         "y"         "z"         "gpstime"   "intensity"

Optionally the list of header components (some extracted as actual
data, some left in their "raw" form) can be returned instead of data
using "returnHeaderOnly=TRUE".

###
## example: simple plot with rgl
library(rgl)
scl <- function(x) (x - min(x))/diff(range(x))
plot3d(d[,1:3], col = topo.colors(256)[scl(d[,5]) * 256])

I hope this is of interest.

Cheers, Mike.




lasinfo output on "lfile.las"
---------------------------------------------------------
 Header Summary
---------------------------------------------------------
 File Name: lfile.las
 Version:                    1.0
 Source ID:                  0
 Reserved:                   0
 Project ID/GUID:           '00000000-0000-0000-0000-000000000000'
 System Identifier:         ''
 Generating Software:       'TerraScan'
 File Creation Day/Year:    0/0
 Header Size                227
 Offset to Point Data       229
 Number Var. Length Records 0
 Point Data Format          1
 Point Data Record Length   28
 Number of Point Records    1922632
 Number of Points by Return 1572923 324305 24659 740 5
 Scale Factor X Y Z         0.01 0.01 0.01
##  etc...