r-lidar / rlas

R package to read and write las and laz files used to store LiDAR data
https://cran.r-project.org/package=rlas
GNU General Public License v3.0
34 stars 14 forks source link

write.las writes 0 point for point data format 8 #60

Open candelas762 opened 1 year ago

candelas762 commented 1 year ago

After adding R-G-B-NIR values to the cloud I can not read the cloud back. It appears like the cloud is empty although the file itself has size.

A reproducible example is this:

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las = readLAS(LASfile)

dir.out= "path to folder"

# Generate dataframe with random RGBI values
rgbi.info=data.frame(Red=  as.integer(sample(seq(from=0,to=100,by=5),size=nrow(las@data),replace=TRUE)),
                Green=  as.integer(sample(seq(from=0,to=100,by=5),size=nrow(las@data),replace=TRUE)),
                Blue=  as.integer(sample(seq(from=0,to=100,by=5),size=nrow(las@data),replace=TRUE)),
                IR=as.integer(sample(seq(from=0,to=100,by=5),size=nrow(las@data),replace=TRUE)))

# Add rgbi values with  
las = lidR::add_lasrgb(las, rgbi.info$Red, rgbi.info$Green, rgbi.info$Blue)
las = lidR::add_lasnir(las, rgbi.info$IR)

# Export
writeLAS(las, paste0(dir.out,"/Megaplot_rgbi.las"))

# Read back into R
las.rgbi = readLAS(paste0(dir.out,"/Megaplot_rgbi.las"))
head(las.rgbi@data)

EDIT:

The problem relies in the add_lasnir() function. If only RGB values are stored, there is no problem.

Jean-Romain commented 1 year ago

It seems to be an issue in the package rlas that wrote 0 point in the file. The reading is correct. The file is written with 0 point.

Jean-Romain commented 1 year ago

Minimal reproducible example without lidR (moving to rlas)

library(rlas)
LASfile <- system.file("extdata", "example.laz", package="rlas")
las = read.las(LASfile)
head = read.lasheader(LASfile)

pdf3 = tempfile(fileext = ".las")
pdf8 = tempfile(fileext = ".las")

# Generate dataframe with random RGBI values
las$Red = as.integer(runif(nrow(las), 0, 2^16-1))
las$Green = as.integer(runif(nrow(las), 0, 2^16-1))
las$Blue = as.integer(runif(nrow(las), 0, 2^16-1))
las$NIR = as.integer(runif(nrow(las), 0, 2^16-1))

head3 = head
head3$`Point Data Format ID` = 3L

head8 = head
head8$`Point Data Format ID` = 8L

write.las(pdf3, head3, las)
#> Warning: Invalid file: the data contains a 'NIR' attribute but point data format
#> is not set to 8 or 10
write.las(pdf8, head8, las)

las3 = read.las(pdf3)
las8 = read.las(pdf8)

head(las3)
#>           X       Y       Z  gpstime Intensity ReturnNumber NumberOfReturns
#> 1: 339002.9 5248001 975.589 269347.3        82            1               1
#> 2: 339003.0 5248000 974.778 269347.3        54            1               1
#> 3: 339002.9 5248000 974.471 269347.3        27            2               2
#> 4: 339002.9 5248000 974.025 269347.3        55            2               2
#> 5: 339003.6 5248000 974.298 269347.3       117            1               1
#> 6: 339003.5 5248000 974.985 269347.3        81            1               1
#>    ScanDirectionFlag EdgeOfFlightline Classification Synthetic_flag
#> 1:                 1                1              1          FALSE
#> 2:                 1                0              1          FALSE
#> 3:                 1                0              1          FALSE
#> 4:                 1                0              1          FALSE
#> 5:                 0                0              1          FALSE
#> 6:                 0                0              1          FALSE
#>    Keypoint_flag Withheld_flag ScanAngleRank UserData PointSourceID R G B
#> 1:         FALSE         FALSE           -21       32            17 0 0 0
#> 2:         FALSE         FALSE           -21       32            17 0 0 0
#> 3:         FALSE         FALSE           -21       32            17 0 0 0
#> 4:         FALSE         FALSE           -21       32            17 0 0 0
#> 5:         FALSE         FALSE           -21       32            17 0 0 0
#> 6:         FALSE         FALSE           -21       32            17 0 0 0
head(las8)
#> Empty data.table (0 rows and 20 cols): X,Y,Z,gpstime,Intensity,ReturnNumber...

Created on 2023-02-02 by the reprex package (v2.0.1)