example syntax for filtering using a created raster?

[Crghilardi]

I'm working through a situation where I need to filter a point cloud to where the points are +/- a tolerance value of some surface I have created.

For instance, if I have some height layer, I want to subset the points to ones within 5m of the surface value. What is the correct way to pass that to filter? I have tried a few ways but couldn't get it working.

using PointCloudRasterizers
using LazIO
using GeoArrays
using Statistics

lazfn = joinpath(dirname(pathof(LazIO)), "..", "test/libLAS_1.2.laz")
pointcloud = LazIO.open(lazfn)

cellsizes = (1.,1.)
raster_index = index(pointcloud, cellsizes)

raster = reduce(raster_index, field=:Z, reducer=median)

# try to filter to points within 5 units of the surface
tolerance = 5.0

#not working?
within_tol(p) =  p.Z < raster + tolerance && p.Z > raster - tolerance

#?
filter!(raster_index, within_tol)

[visr]

Hmm yeah this usage currently falls between the cracks. reduce reduces it to a grid, and filter filters the points based only on the point itself. What you want to do is like filter, but based not only on the point but also based on the value of a raster on that location. I tried to prototype something quickly, but it's not working yet:

function filter_raster(index::PointCloudRasterizers.PointCloudIndex, raster::GeoArray, condition)
    @showprogress 1 "Reducing points..." for (i, p) in enumerate(index.ds)
        @inbounds ind = index.index[i]
        ind == 0 && continue  # filtered point
        raster_value = raster[ind]

        if ~condition(p, raster_value)
            @inbounds index.counts.A[ind] -= 1
            @inbounds index.index[i] = 0
        end
    end
end

within_tol(p, raster_value) = isapprox(p.Z, raster_value, atol=5.0)
filter_raster(raster_index, raster, within_tol)

Problem is that I get on the first point already a missing raster_value, not sure how that is possible. But haven't had the time to look into it more.

Such a function may be a nice addition.

[Crghilardi]

Thanks for taking a look. I remember the example data above (libLAS_1.2.laz) has many missing cells for some reason.I tried the code above with a different point cloud and got the same result. I don't know if I have any dense enough point clouds to test without missing cells.

I'll keep trying things on my end. If you have any ideas, feel free to list them here.

[visr]

Yeah I'm not so surprised about the presence of missing cells, since any cell in the bounding box of the point cloud that doesn't have points will have missing. But in my code I use the cell index associated to a point, so it can by definition not be empty. I'm probably making a mistake though.

[Crghilardi]

It looks like the issue was missings have some comparison logic rules and just return missing rather than a boolean back to you.

isapprox(1800,missing,atol=5.0) #returns missing
isapprox(1800,1801,atol=5.0) #returns true

I was able to get around it by adding a method to Base.isapprox and having it return false. I kept the type signature to just the three arguments, not sure if that helps to quarantine this without running into issues.

function Base.isapprox(x::Number, y::Missing; atol::Real=0)
    false
end

After that, your code above runs.

[visr]

Ah but my issue is that raster_value becomes missing, and I cannot see how that can be, since the raster should not be missing if it has one or more points in the cell. And we know it has at least one point in the cell because the current point p is in that cell.

Regarding the isapprox, it's probably better not to redefine Base methods for Base types, a lot of things can break as a result. I would just put it in a new function with a new name.

[Crghilardi]

Ahhh OK, I understand your question now. It didn't click right away...

I tried to find some specific examples where raster_index.counts does not match up with raster using the example data above.

#does not match, raster value with no points counted
julia> raster[500,500]
1-element Array{Union{Missing, Float64},1}:
 842.3000000000001

julia> raster_index.counts[500,500]
1-element Array{Int64,1}:
 0

#does not match, missing raster but points counted
julia> raster[3237,359]
1-element Array{Union{Missing, Float64},1}:
 missing

julia> raster_index.counts[3237,359]
1-element Array{Int64,1}:
 2

#matches missing and no points counted
julia> raster_index.counts[700,700]
1-element Array{Int64,1}:
 0

julia> raster[700,700]
1-element Array{Union{Missing, Float64},1}:
 missing

[Crghilardi]

Got it! The issue was line 133 in PointCloudRasterizers.jl, there was a image flip on the reduced output but not index.counts so they ended up out of sync. I commented the flipud! out and then the counts layer matches the reduced layer and thefilter_raster snippet above works without needing to do any of the isapprox hacking I did before.

ga = GeoArray(output, index.counts.f, index.counts.crs)
#GeoArrays.flipud!(ga)  # move to GeoArrays

I am not sure if it is better to get rid of the flip or add another flip on counts as well?

[visr]

Yeah that sounds logical, nice find!

I think ideally we git rid of flips and just set up our PointCloudIndex such that it always matches the y direction that you want your grid in, such that there is no need for any flips.

[evetion]

I agree that flipud should be removed, I put it in a comment on the same line (move to GeoArrays.jl). This was mainly done because GDAL can't always work with positive y directions (can't include such tiff files in a .vrt for example)

This packages tries to keep a PointCloudIndex consistent (hence the struct), but assumes that products of a reduce are not linked to that index anymore as it cannot guarantee that.

I like @visr's approach of filter_raster (or even filter with a GeoArray). Although the same checks as in evetion/GeoArrays.jl#19 would have to apply, only the same size/bounds as the .count array.

[Crghilardi]

OK, I can a PR together for all this sometime soon.

This also explains the issue I ran into a while ago where raster = reduce(raster_index, field=:Z, reducer = length) did not match the raster_index.counts

Does it make sense to add a test to check that all indices with missing have a count of 0?