brycefrank / pyfor

Tools for analyzing aerial point clouds of forest data.
MIT License
93 stars 19 forks source link

normalization produces blank cloud (missing data?) #29

Closed bw4sz closed 5 years ago

bw4sz commented 5 years ago

I've noticed that on occasion the normalization function can strip the entire cloud of points.

For example here it works.

import pyfor
from matplotlib import pyplot as plt
pc=pyfor.cloud.Cloud("/Users/ben/Documents/DeepLidar/data/SJER/SJER_006.laz")

pc.plot()
pc.normalize(3)
pc.plot()

figure_5 figure_4

But here it doesn't

import pyfor
from matplotlib import pyplot as plt
pc=pyfor.cloud.Cloud("/Users/ben/Documents/DeepLidar/data/SJER/SJER_004.laz")

pc.plot()
pc.normalize(3)
pc.plot()

figure_1 figure_2

These were crops from a larger tile (I can reprocess if I must). It looks like some missing data in the corner might be to blame? Just a guess. But then again

pc.data.z.isna().sum()
0
pc.data.points.z
0      NaN
1      NaN
2      NaN
3      NaN
4      NaN
5      NaN
6      NaN

data can be found here: https://github.com/weecology/NeonTreeEvaluation/tree/master/SJER

brycefrank commented 5 years ago

Interesting. I have a guess at what is wrong, and I think your theory about the missing points may be valid. I will take a look at the data and get back to you tonight.

brycefrank commented 5 years ago

Ben,

The default ground filter is KrausPfeifer1998 which uses a larger pixel size for the intermediate surfaces than other filters. As an artifact of this, some interesting things are happening when computing the bare earth model for such a small area. It's not to say it isn't broken, but I think there are better options.

I think what may be best for your immediate use case is using the Zhang2003 ground filter, which is about 80% implemented (it is in fact on my to do list). It is more robust to smaller plot areas because it can effectively filter using finer resolution intermediate surfaces. However, its implementation has been a bit tricky in some fringe-cases, which is why it isn't complete. I will see if I can't crunch it through here in the next few hours. If not today it will be tomorrow.

Edit

I have put a little time into Zhang2003, there were some obvious mistakes that a bit of time away from the algorithm has clarified. Unfortunately, I am a bit stacked with other things tomorrow. If I get the chance in the morning I will see if I can cobble together a sensible fix for you.

brycefrank commented 5 years ago

Here is what I came up with (note that I converted to las just to get it working on a machine without lastools). Please pull from master to get the update via https://github.com/brycefrank/pyfor/commit/dd346bde3facf195bac0cb8761dce7923a68a11c

Note that I have not tested the new fix thoroughly, keep an eye out for any weirdness.

pc = pyfor.cloud.Cloud('/home/bryce/Desktop/pyfor_testing/SJER_004.las')
pc.plot()

image

zhang_filter=pyfor.ground_filter.Zhang2003(pc, cell_size=1)
zhang_filter.normalize()

pc.plot()

image

Note that the Zhang2003 filter must be instantiated manually (that is the previous code block). zhang_filter.normalize() filters the cloud in place. This is (hopefully) consistent with the default function.

For more details on the Zhang filter here is the original paper and my somewhat scant documentation on it. It may take a little tuning.

http://users.cis.fiu.edu/~chens/PDF/TGRS.pdf http://brycefrank.com/pyfor/html/source/pyfor.ground_filter.html#pyfor.ground_filter.Zhang2003

As for the KrausPfeifer1998 issues I will leave this issue open as I look into it. I should have time this week.

bw4sz commented 5 years ago

Looking good so far. Thanks.

brycefrank commented 5 years ago

KrausPfeifer1998 is giving me some headaches. I am not sure if it is an artifact of the tile size anymore. It does have two pesky hyperparameters, (g and w) that may be leading to some weird behavior I am getting. I may switch back the default to Zhang2003 now that it is fixed in 0.3.2. I will keep this open while I continue to debug and ponder.

brycefrank commented 5 years ago

Ben,

A funny story. The other night, after working on this for a good chunk of the day, I woke up. It was dark. It must have been 2 am. I had the solution! When debugging KrausPfeifer1998 I switched from filling a 3D matrix of Z values with nans (for missing) to filling a 3D matrix of Z values with 0s (for missing). This causes the residuals computed in the filter to be incorrect. I quickly fell back into a deep slumber.

Today, I suddenly remembered my midnight revelation, and the fix has been implemented in 0.3.2.

There is still some small edge effect issues for small tiles, but this is really more of a shortcoming of the filter than anything else. It could be alleviated by using a buffered cloud to compute the DEM. I left a small TODO if I have time to look at it later.