在以下代码中，可以导入具有这些属性的部分 LAS 文件版本 1.1：
列表。项目格式 尺寸要求
- X长4字节*
- Y长4字节*
- Z长4字节*
- 强度无符号短 2 字节
- 返回编号 3 位（位 0、1、2） 3 位 *
- 返回数（给定脉冲） 3 位（位 3、4、5） 3 位 *
- 扫描方向标志 1 位 (位 6) 1 位 *
- 飞行线边缘 1 位（位 7） 1 位 *
- 分类 unsigned char 1 字节 *
- 扫描角度等级（-90 至 +90） – 左侧无符号字符 1 字节 *
- 用户数据 unsigned char 1 字节
- 点源 ID 无符号短 2 字节 *
- GPS 时间双 8 字节 *
并返回 x、y、z 和强度（上面列表的第 1 行到第 4 行）：
allbytes <- matrix(readBin(con, "raw", n = pointDataRecordLength * numberPointRecords, size = 1, endian = "little"),
ncol= pointDataRecordLength, nrow = numberPointRecords, byrow = TRUE)
close(con)
mm <- matrix(readBin(t(allbytes[,1:(3*4)]), "integer", size = 4, n = 3 * numberPointRecords, endian = "little"), ncol = 3, byrow = TRUE)
mm[,1] <- mm[ ,1] * xyzScaleOffset[1,1] + xyzScaleOffset[1, 2]
mm[,2] <- mm[ ,2] * xyzScaleOffset[2,1] + xyzScaleOffset[2, 2]
mm[,3] <- mm[ ,3] * xyzScaleOffset[3,1] + xyzScaleOffset[3, 2]
colnames(mm) <- c("x", "y", "z")
intensity <- readBin(t(allbytes[, 13:14]), "double", size = 2, n = numberPointRecords, signed = FALSE, endian = "little")
我想扩展代码以读取“返回号码”和“返回号码”（表的第 5 行和第 6 行）。我做错了什么？
returnNumber <- readBin(t(allbytes[, 15]), "integer", size = 1, n = numberPointRecords, signed = FALSE, endian = "little")
预期回报应为 1 到 5 之间的整数。提前致谢！
示例文件：
link http://www.filefactory.com/file/68sdgepk12pd/n/lidar_las
完整代码：
publicHeaderDescription <- function() {
hd <- structure(list(Item = c("File Signature (\"LASF\")",
"(1.1) File Source ID", "(1.1) Global Encoding",
"(1.1) Project ID - GUID data 1", "(1.1) Project ID - GUID data 2",
"(1.1) Project ID - GUID data 3", "(1.1) Project ID - GUID data 4",
"Version Major", "Version Minor", "(1.1) System Identifier",
"Generating Software", "(1.1) File Creation Day of Year",
"(1.1) File Creation Year", "Header Size", "Offset to point data",
"Number of variable length records",
"Point Data Format ID (0-99 for spec)", "Point Data Record Length",
"Number of point records", "Number of points by return",
"X scale factor", "Y scale factor", "Z scale factor", "X offset",
"Y offset", "Z offset", "Max X", "Min X", "Max Y", "Min Y", "Max Z",
"Min Z"), Format = c("char[4]", "unsigned short", "unsigned short",
"unsigned long", "unsigned short", "unsigned short",
"unsigned char[8]", "unsigned char", "unsigned char", "char[32]",
"char[32]", "unsigned short", "unsigned short", "unsigned short",
"unsigned long", "unsigned long", "unsigned char", "unsigned short",
"unsigned long", "unsigned long[5]", "double", "double", "double",
"double", "double", "double", "double", "double", "double", "double",
"double", "double"), Size = c("4 bytes", "2 bytes", "2 bytes",
"4 bytes", "2 byte", "2 byte", "8 bytes", "1 byte", "1 byte",
"32 bytes", "32 bytes", "2 bytes", "2 bytes", "2 bytes", "4 bytes",
"4 bytes", "1 byte", "2 bytes", "4 bytes", "20 bytes", "8 bytes",
"8 bytes", "8 bytes", "8 bytes", "8 bytes", "8 bytes", "8 bytes",
"8 bytes", "8 bytes", "8 bytes", "8 bytes", "8 bytes"), Required =
c("*", "*", "*", "*", "*", "*", "*", "*", "*", "*", "*", "*", "*",
"*", "*", "*", "*", "*", "*", "*", "*", "*", "*", "*", "*", "*", "*",
"*", "*", "*", "*", "*")), .Names = c("Item", "Format", "Size",
"Required"), row.names = 2:33, class = "data.frame")
hd$what <- ""
hd$what[grep("unsigned", hd$Format)] <- "integer"
hd$what[grep("char", hd$Format)] <- "raw"
hd$what[grep("short", hd$Format)] <- "integer"
hd$what[grep("long", hd$Format)] <- "integer"
hd$what[grep("double", hd$Format)] <- "numeric"
hd$signed <- TRUE
hd$signed[grep("unsigned", hd$Format)] <- FALSE
## number of values in record
hd$n <- as.numeric(gsub("[[:alpha:][:punct:]]", "", hd$Format))
hd$n[hd$what == "character"] <- 1
hd$n[is.na(hd$n)] <- 1
## size of record
hd$Hsize <- as.numeric(gsub("[[:alpha:]]", "", hd$Size))
## size of each value in record
hd$Rsize <- hd$Hsize / hd$n
hd$Rsize[hd$what == "raw"] <- 1
hd$n[hd$what == "raw"] <- hd$Hsize[hd$what == "raw"]
hd
}
readLAS <-
function(lasfile, skip = 0, nrows = NULL, returnSP = FALSE, returnHeaderOnly = FALSE) {
hd <- publicHeaderDescription()
pheader <- vector("list", nrow(hd))
names(pheader) <- hd$Item
con <- file(lasfile, open = "rb")
isLASFbytes <- readBin(con, "raw", size = 1, n = 4, endian = "little")
pheader[[hd$Item[1]]] <- readBin(isLASFbytes, "character", size = 4, endian = "little")
if (! pheader[[hd$Item[1]]] == "LASF") {
stop("Not a valid LAS file")
}
for (i in 2:nrow(hd)) {
pheader[[hd$Item[i]]] <- readBin(con, what = hd$what[i], signed = hd$signed[i], size = hd$Rsize[i], endian = "little", n = hd$n[i])
#print(names(pheader)[i])
#print(pheader[[hd$Item[i]]])
}
close(con)
## read the data
numberPointRecords <- pheader[["Number of point records"]]
offsetToPointData <- pheader[["Offset to point data"]]
pointDataRecordLength <-pheader[["Point Data Record Length"]]
xyzScaleOffset <- cbind(unlist(pheader[c("X scale factor", "Y scale factor", "Z scale factor")]),
unlist(pheader[c("X offset", "Y offset", "Z offset")]))
if (returnHeaderOnly) return(pheader)
con <- file(lasfile, open = "rb")
junk <- readBin(con, "raw", size = 1, n = offsetToPointData)
## deal with rows to skip and max rows to be read
if (skip > 0) {
## seek is unreliable on windows, or I'm using it incorrectly
## so we junk the bytes to skip
junk <- readBin(con, "raw", size = 1, n = pointDataRecordLength * skip)
numberPointRecords <- numberPointRecords - skip
#pos <- seek(con, where = pointDataRecordLength * skip)
# print(c(pos = seek(con), skip = skip, where = pointDataRecordLength * skip))
}
if (!is.null(nrows)) {
if (numberPointRecords > nrows) numberPointRecords <- nrows
}
if (numberPointRecords < 1) stop("no records left to read")
# include a loop to read just points inside the x and y coordinates
allbytes <- matrix(readBin(con, "raw", n = pointDataRecordLength * numberPointRecords, size = 1, endian = "little"),
ncol= pointDataRecordLength, nrow = numberPointRecords, byrow = TRUE)
close(con)
mm <- matrix(readBin(t(allbytes[,1:(3*4)]), "integer", size = 4, n = 3 * numberPointRecords, endian = "little"), ncol = 3, byrow = TRUE)
gpstime <- NULL
if (ncol(allbytes) == 28) gpstime <- readBin(t(allbytes[ , 21:28]), "numeric", size = 8, n = numberPointRecords, endian = "little")
intensity <- readBin(t(allbytes[, 13:14]), "double", size = 2, n = numberPointRecords, signed = FALSE, endian = "little")
mm[,1] <- mm[ ,1] * xyzScaleOffset[1,1] + xyzScaleOffset[1, 2]
mm[,2] <- mm[ ,2] * xyzScaleOffset[2,1] + xyzScaleOffset[2, 2]
mm[,3] <- mm[ ,3] * xyzScaleOffset[3,1] + xyzScaleOffset[3, 2]
colnames(mm) <- c("x", "y", "z")
returnNumber <- readBin(t(allbytes[,15]), "integer", size = 1, n = numberPointRecords, signed = FALSE, endian = "little")
require(bitops)
if (returnSP) {
require(sp)
SpatialPoints(cbind(mm, gpstime, intensity))
} else {
cbind(mm, gpstime, intensity)
}
}
