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...
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
