r-lidar / lidR

Airborne LiDAR data manipulation and visualisation for forestry application
https://CRAN.R-project.org/package=lidR
GNU General Public License v3.0
587 stars 132 forks source link

normalize_intensity Internal error: access to coordinates below 0 in the array #608

Closed candelas762 closed 2 years ago

candelas762 commented 2 years ago

Not very sure what is the problem here. I can not find other issues or questions regarding this error.

You can download the LAS file that produces the error here.

library(lidR)
las = readLAS("C:/Users/plot899error.las", filter = "-drop_class 7")
sensor = track_sensor(las, Roussel2020(pmin = 500), multi_pulse = T, thin_pulse_with_time=0)
las_dI = normalize_intensity(las, range_correction(sensor, Rs = 1300)) # Rs = average flight height
#> 3 coordinates detected in the sensor positions, parameter 'elevation' is not considered.
#> Error: Internal error: access to coordinates below 0 in the array. Please report this bug
Jean-Romain commented 2 years ago

What I can say right now is that there is a conditional statement somewhere to protect against serious crash for a case that is not supposed to happen. Yet it happened and your R session did not crash so it is a success :smile: Now I need to figure out why it happened.

Jean-Romain commented 2 years ago

Without investigating yet, here what I think. You cannot apply normalized_intensity() using the output of track_sensor() on the same point cloud. You need a buffer to get extra sensor position and context. Exactly like you need a buffer to compute a DTM otherwise you have edge artifacts. For a DTM, with our without buffer it works the same but in the later case your DTM is incorrect and inaccurate. Here this is the same idea but the fact that you missed some spatial context is more damageable. The code is expected to handle this case and apply a poor interpolation instead of nothing actually but I guess you reached a limit cases not supported. I'll try to fix it but in does not invalidate my point. You need a buffer to compute a valid normalization.

candelas762 commented 2 years ago

Since I am trying to make the process as computationally-fast as possible I added a very small buffer (it needs to be bigger than 0 so I used 0.5... probably too small 😅 ) to just process the catalog tiles covering target plots. If I wanted to keep this strategy, should I apply track_sensor() to the selected catalog tiles and apply normalize_intensity() to a smaller clipped cloud around the plot? I tried this and it solves the problem but I guess that is not a final solution for this error.

My problem is that it seems like I have the worst data for intensity normalization. So far I have came across al the documented errors and now not documented ones. I found that is easier to select the tile/s covering the target plots and fix them one at a time rather than trying to fix the whole catalog problems.

Jean-Romain commented 2 years ago

There are two parts. The computation of the positions (track_sensor()) and the normalization from the positions (normalize_intensity(). In both you have potential problem with edges.

The doc of track_sensor()

For LAScatalog processing it is recommended to use large chunks and large buffers (e.g. a swath width). The point cloud must not be normalized.

But obviously at the edges of your file collection you can't have a buffer and the positions are necessarily weaker. You can't do anything against that because you don't have extra data. Then you apply normalize_intensity() and you have problems at the edges. The only correct options is to reduce slightly your area and normalize only your collection minus e.g. 500 meters. But in practice it should not fail.

I guess that is not a final solution for this error.

It is the correct way if you want to be sure all your values are decently correct. At the very edges it will always be poor. Yet, as I said, it should not fail computationally speaking. So I'll try to investigate on the failure to return at least a more informative message.

My problem is that it seems like I have the worst data for intensity normalization. So far I have came across al the documented errors and now not documented ones. I found that is easier to select the tile/s covering the target plots and fix them one at a time rather than trying to fix the whole catalog problems.

The number of way to fail in these two functions is virtually infinite. This is problematic on my side too. Having correct XYZ coordinates in a point cloud is trivial and the most important part. Yet many point clouds are corrupted and invalid because data providers do not all pay attention to the specifications. So imagine, here, you need to have XYZ, return number, number of returns, gpstime and sometime UserData valid and consistent. This should not be a problem but in practice it is. This is why there are so many tests and so many possible error (there is more error handling code than actual processing code).

candelas762 commented 2 years ago

By selecting the tile (~0.5 km^2) coverig the target plot and the adjacent tiles to use track_sensor() on and then applying normalize_intensity to just a 250 m buffer around the target plot the amount of errors have diminish considerably. Ocasionally I have to adjust the gpstime" range and theReturnNumber` on the clouds, but that is a problem of my badly populated data. Thanks for the tips.

I still get very rarely the error regarding "Internal error: access to coordinates below 0 in the array [...]".

Jean-Romain commented 2 years ago

Let look at the following picture

psid844 = filter_poi(las, PointSourceID == 844)
plot(sensor["PointSourceID"], pch = 19)
plot(sf::st_as_sf(psid844), add = TRUE, col = "black", cex = 0.25, pch = 19)

The black points are lidar points for lightline 844

Rplot02

You have a single point for flightline 844. There is no surprise it crashes because it is impossible to find the position of the sensor at an instant t even with interpolation. Indeed interpolation is not possible. I'll investigate more to figure out where I can put warning and safegards and where I can catch this case. And I don't understand why you have a single position. It seems there are enough point to get more positions

Jean-Romain commented 2 years ago

This one gave me trouble but I fixed it by rebuilding the range function from scratch. The problem was the one I already described. You have only one sensor position for some flightline. This was problematic and should not arise if you pay attention at using large enough area that encompasses more than your working area.

That being said the bug is solved. Now if there is only one sensor position it uses this sensor position without interpolation no matter how far it is. It does not fail but be careful the points can be very far from the sensor position and return incorrect range. The best is to have a set of sensor position that encompass more than your working area. If you do not have that you need at least 2 positions to get a valid interpolation. See the examples below

library(lidR)
las = readLAS("plot899error.las", filter = "-drop_class 7")
sensor = track_sensor(las, Roussel2020(pmin = 500), multi_pulse = T, thin_pulse_with_time=0)
r = get_range(las, sensor)
hist(r)
las  = add_attribute(las, r, "range")
plot(las, color = "range")
library(lidR)
las = readLAS("plot899error.las", filter = "-drop_class 7")
sensor = track_sensor(las, Roussel2020(pmin = 100), multi_pulse = T, thin_pulse_with_time=0)
r = get_range(las, sensor)
hist(r)
las  = add_attribute(las, r, "range")
plot(las, color = "range")
candelas762 commented 2 years ago

That explains another error that I was having regarding "An high range R has been computed" error. It was interpolating using very far away sensor position points that hardly had to do anything with my target plot. I just dropped those points.

Thanks a lot!