3

I have two .las datasets; one has been pre-processed but has limited spatial coverage and the other is the raw data with full coverage.

My analysis only requires points between 10 and 40 m, hence I applied the following filter:

lasfilt = lasfilter(las, Z > 10 & Z < 40)

Returns the error message when plotted:

Error in plot.LAS(x, y, color, colorPalette, bg, trim, backend, clear_artifacts, : Cannot display an empty point cloud

The function works on the dataset that has been pre-processed but not on the raw data. I have read the accompanying documentation but cannot work out why it does not work as I use exactly the same code. In fact the filter does not work with any condition. The main difference between the datasets is the vegetation is classified in the pre-processed data and it is retiled.

Unfortunately, I do not have access to the software used to pre-process the data. How can I apply this filter to the raw data?

Update:

npoints(las)
Returns a value of: 17586563

This is the summary:

class        : LAS (LASF v1.2)
point format : 1
memory       : 402.5 Mb 
extent       :414000, 415000, 558000, 559000 (xmin, xmax, ymin, ymax)
coord. ref.  : NA 
area         : 1 kunits²
points       : 17.59 million points
density      : 17.67 points/units²
names        : X Y Z 
File signature:           LASF 
File source ID:           0 
Global encoding:
 - GPS Time Type: GPS Week Time 
 - Synthetic Return Numbers: no 
 - Well Know Text: CRS is GeoTIFF 
 - Aggregate Model: false 
Project ID - GUID:        00000000-0000-0000-0000-000000000000 
Version:                  1.2
System identifier:         
Generating software:      TerraScan 
File creation d/y:        0/0
header size:              227 
Offset to point data:     229 
Num. var. length record:  0 
Point data format:        1 
Point data record length: 28 
Num. of point records:    17586563 
Num. of points by return: 14588291 2474276 485330 38666 0 
Scale factor X Y Z:       0.01 0.01 0.01 
Offset X Y Z:             0 0 0 
min X Y Z:                414000 558000 -610.28 
max X Y Z:                415000 559000 478.67 
Variable length records:  void

The npoints from the filtered dataset is 0.

This is just one of the .las files I am focusing on to get it working on from a catalog of 32 files.

Applying summary statistics to the Z values:

f = function(x) {list(mean = mean(x), max = max(x))}
grid_metrics(las1, ~f(Z), res = 20)

Returns:

class      : RasterBrick 
dimensions : 50, 50, 2500, 2  (nrow, ncol, ncell, nlayers)
resolution : 20, 20  (x, y)
extent     : 414000, 415000, 558000, 559000  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      :     mean,      max 
min values : 132.1497, 140.6000 
max values : 225.5666, 478.6700

I suspect given the minimum is greater than what I am filtering, this is why the point cloud is returned empty?

For completeness I ran the same as above on some of the pre-processed data:

npoints(las_prepro)
summary(las_prepro)
grid_metrics(las_prepro, ~f(Z), res = 20)


[1] 26583034


class        : LAS (LASF v1.2)
point format : 1
memory       : 608.4 Mb 
extent       :413000, 414000, 558071.2, 559000 (xmin, xmax, ymin, ymax)
coord. ref.  : NA 
area         : 340146.1 units²
points       : 26.58 million points
density      : 78.15 points/units²
names        : X Y Z 
File signature:           LASF 
File source ID:           0 
Global encoding:
 - GPS Time Type: GPS Week Time 
 - Synthetic Return Numbers: no 
 - Well Know Text: CRS is GeoTIFF 
 - Aggregate Model: false 
Project ID - GUID:        00000000-0000-0000-0000-000000000000 
Version:                  1.2
System identifier:         
Generating software:      TerraScan 
File creation d/y:        298/2012
header size:              227 
Offset to point data:     229 
Num. var. length record:  0 
Point data format:        1 
Point data record length: 28 
Num. of point records:    26583034 
Num. of points by return: 17664187 7698339 1179992 40516 0 
Scale factor X Y Z:       0.01 0.01 0.01 
Offset X Y Z:             0 0 0 
min X Y Z:                413000 558071.2 -110.81 
max X Y Z:                414000 559000 423.38 
Variable length records:  void


class      : RasterBrick 
dimensions : 47, 50, 2350, 2  (nrow, ncol, ncell, nlayers)
resolution : 20, 20  (x, y)
extent     : 413000, 414000, 558060, 559000  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      :       mean,        max 
min values : -0.1692308, -0.1300000 
max values :   14.31257,  423.38000 

Script ran from answer:

> # Use a RasterLater as DTM and normalize points
> dtm <- grid_terrain(las, 1, kriging(k = 10L))

Error in grid_terrain.LAS(las, 1, kriging(k = 10L)) : 
  LAS object does not contain 'Classification' data

# Attempting to classify the ground points 

> las_groundcsf <- lasground(las, csf())

Warning message:
'ReturnNumber' and/or 'NumberOfReturns' not found. Cannot use the option 'last_returns', all the points will be used.

# Method number two
ws  <- seq(3,12, 3)
th  <- seq(0.1, 1.5, length.out = length(ws))
las <- lasground(las, pmf(ws, th))

'ReturnNumber' and/or 'NumberOfReturns' not found. Cannot use the option 'last_returns', all the points will be used.
  • 1
    When you say it doesn't work, what happens? Error message? No points returned? Wrong points returned? If you cant release the data so we can replicate this, can you at least show us the summary output you get when you print them? – Spacedman Jan 3 at 22:13
  • The line of code above produces no error messages but in the environment you can see the file size is exactly the same when it should be much smaller. Then when you try plot(lasfilt) the following error message is returned: Error in plot.LAS(x, y, color, colorPalette, bg, trim, backend, clear_artifacts, : Cannot display an empty point cloud – Emma Jan 3 at 22:50
  • 1
    What does npoints say on your lidar objects? What do the basic summary statistics say? What does summary say? Please edit your question to show all this as well as the information in your previous comment. – Spacedman Jan 3 at 22:54
  • 1
    Thank you very much for your help so far, I have added the extra information to my question – Emma Jan 3 at 23:42
  • You can inspect the Z values using lasfilt$Z, and so count how many should be in that range with something like sum(las$Z > 10 & las$Z<40). If that's zero then there's no points in that range. Do a few more tests to see what the range of Z is, try a histogram and so on. But that would seem to be the problem. – Spacedman Jan 3 at 23:49
1

You need to normalize the raw point cloud data so that points are measured above ground level rather than above sea level. Your current filter lasfilter(las, Z > 10 & Z < 40) is not being applied to your point cloud data which has the following elevation values:

names      :     mean,      max 
min values : 132.1497, 140.6000 
max values : 225.5666, 478.6700

The lidR documentation (p.64) shows an example of normalizing point cloud data using lasnormalize():

library(lidR)

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

# Use a RasterLayer as DTM and normalize points
dtm <- grid_terrain(las, 1, kriging(k = 10L))
las_n <- lasnormalize(las, dtm)

# Check to see if las is normalized now
lascheck(las_n)
plot(las_n) 
| improve this answer | |
  • Thank you for this; that makes sense. I have run the script, however as the .las is unclassified there is an error message in the dtm line. I have copied the script ran and message returned to my quesiton if you wouldn't mind having a look. – Emma Jan 4 at 13:02
  • @Emma Have you tried using lasground() to classify the points into ground and non-ground points? – Aaron Jan 8 at 15:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.