How do I calculate a surface roughness index of multiple lidar point clouds with lidR rumple_index?
library(lidR)
files <- list.files(path= "/pathToFiles", pattern= "*.laz", full.names = TRUE, recursive = FALSE)
O = lapply(files, function(x) {
las = readLAS(x, select = "xyzc", filter = "-drop_z_below 0")
chm = grid_canopy(las, 0.2, p2r())
roughness <- rumple_index(chm)
return(data.frame(file=x,roughness=roughness))
})
O = do.call(rbind,O)
However this code returns error:
Error in C_rasterize(las, layout, subcircle, 1L) :
C++ unexpected internal error in 'rasterize': point out of raster.
In addition: Warning message:
In showSRID(SRS_string, format = "PROJ", multiline = "NO") :
If I run the following code for a single file, it gives reasonable result.
las = readLAS("/input.laz", select = "xyzicnrRGB", filter = "keep_first -drop_z_below 0")
chm = grid_canopy(las, 0.5, p2r())
rumple_index(chm)
I found the file which was causing this error, below is a lasinfo report of that file.
lasinfo (190812) report for 'infile.las'
reporting all LAS header entries:
file signature: 'LASF'
file source ID: 0
global_encoding: 0
project ID GUID data 1-4: 00000000-0000-0000-0000-000000000000
version major.minor: 1.4
system identifier: 'LAStools (c) by rapidlasso GmbH'
generating software: 'las2las (version 190623)'
file creation day/year: 0/0
header size: 375
offset to point data: 469
number var. length records: 1
point data format: 6
point data record length: 30
number of point records: 0
number of points by return: 0 0 0 0 0
scale factor x y z: 0.00025 0.00025 0.00025
offset x y z: 698292.56999999995 6104003.0700000003 0
min x y z: 698607.00800 6102989.34350 -0.03375
max x y z: 698698.55000 6103195.27825 19.08325
start of waveform data packet record: 0
start of first extended variable length record: 0
number of extended_variable length records: 0
extended number of point records: 363189
extended number of points by return: 363189 0 0 0 0 0 0 0 0 0 0 0 0 0 0
variable length header record 1 of 1:
reserved 0
user ID 'LASF_Projection'
record ID 34735
length after header 40
description 'by LAStools of rapidlasso GmbH'
GeoKeyDirectoryTag version 1.1.0 number of keys 4
key 1024 tiff_tag_location 0 count 1 value_offset 1 - GTModelTypeGeoKey: ModelTypeProjected
key 3072 tiff_tag_location 0 count 1 value_offset 28355 - ProjectedCSTypeGeoKey: GDA94 / MGA 55S
key 3076 tiff_tag_location 0 count 1 value_offset 9001 - ProjLinearUnitsGeoKey: Linear_Meter
key 4099 tiff_tag_location 0 count 1 value_offset 9001 - VerticalUnitsGeoKey: Linear_Meter
reporting minimum and maximum for all LAS point record entries ...
X 1257752 1623920
Y -4054906 -3231167
Z -135 76333
intensity 0 42532
return_number 1 1
number_of_returns 1 1
edge_of_flight_line 0 1
scan_direction_flag 0 0
classification 1 1
scan_angle_rank 0 0
user_data 0 0
point_source_ID 0 0
gps_time 164539.903446 168324.056522
extended_return_number 1 1
extended_number_of_returns 1 1
extended_classification 1 1
extended_scan_angle 0 0
extended_scanner_channel 0 0
number of first returns: 363189
number of intermediate returns: 0
number of last returns: 363189
number of single returns: 363189
overview over extended number of returns of given pulse: 363189 0 0 0 0 0 0 0 0 0 0 0 0 0 0
histogram of classification of points:
363189 unclassified (1)
I am using lidR_3.0.1 raster_3.1-5 sp_1.4-2
unexpected internal error
. This message speaks for itself. One of your point-clouds fall into an unexpected behavior partially handled internally. Please find which one, make a reproducible example and report this issue. Has an hint to make a minimal reproducible example: it comes fromp2r()
and is not related torumple_index()
. Also please mention the version and ensure you run a recent version of the package. This bug has already been fixed. It might be a new one but maybe not.print(x)
in your function and you will find the file that failsprint(x)
didn't find the file which was failing. Thanks for your help.