Add COPC reading support

[tmontaigu]

This adds the CopcReader class to laspy
that understands the COPC LAZ format, enabling to
do spatial and / or LOD based queries.

The class require lazrs to work,
and optionaly, if requests is installed
COPC over HTTP is supported.

[hobu]

In addition to filtering by a bounds, it would also be really useful to filter by a resolution or tree level. PDAL has some examples of this resolution-based filtering at https://github.com/PDAL/PDAL/tree/master/io/private/copc and https://github.com/PDAL/PDAL/blob/master/io/CopcReader.cpp#L241

[tmontaigu]

In addition to filtering by a bounds, it would also be really useful to filter by a resolution or tree level. PDAL has some examples of this resolution-based filtering at https://github.com/PDAL/PDAL/tree/master/io/private/copc and https://github.com/PDAL/PDAL/blob/master/io/CopcReader.cpp#L241

The query function is defined as

    def query(
        self,
        bounds: Optional[Bounds] = None,
        level: Optional[Union[int, range]] = None,
    ) -> ScaleAwarePointRecord:

So you can filter by bounds or level or both at the same time , some wrapper functions also exists

    def spatial_query(self, bounds: Bounds) -> ScaleAwarePointRecord:
        return self.query(bounds=bounds, level=None)

    def level_query(self, level: Union[int, range]) -> ScaleAwarePointRecord:
        return self.query(bounds=None, level=level)

[hobu]

A helper function to allow users to convert level->resolution would be very useful.

[tmontaigu]

Ok I got confused, resolution != level

[hobu]

Ok I got confused, resolution != level

Yes. They're directly related, but we have found with EPT and COPC it is much easier to say "select everything at or below 3m" instead of having to calculate based on the dimensions of the octree what level would correspond to everything 3m or below.

[hobu]

Some more files to play with

https://viewer.copc.io/?copc=https://s3.amazonaws.com/hobu-lidar/montreal-2015.copc.laz
https://viewer.copc.io/?copc=https://s3.amazonaws.com/hobu-lidar/melbourne-2018.copc.laz

[tmontaigu]

resolution was added as a query parameter

[hobu]

I attempted to test this but got the following:

AttributeError: module 'lazrs' has no attribute 'decompress_points_with_chunk_table'

I presume this is due to lazrs not being released with chunk table support quite yet?

[tmontaigu]

Yes, I plan on releasing lazrs later this week, after that this PR will be testable easily

[tmontaigu]

The new version of lazrs was pushed to pypi so it should be easier to test now

[hobu]

OSX arm wheels for lazrs didn't post https://pypi.org/project/lazrs/#files

[hobu]

I installed my own lazrs.

It seems a bit slower than PDAL while fetching over the network. Are you reading the entire octree down to the client?

time pdal info --all https://s3.amazonaws.com/hobu-lidar/melbourne-2018.copc.laz   
1.59s user 1.48s system 110% cpu 2.786 total

time python copc.py https://s3.amazonaws.com/hobu-lidar/melbourne-2018.copc.laz
copc_reader_ready for https://s3.amazonaws.com/hobu-lidar/melbourne-2018.copc.laz
Spatial Query gave: <ScaleAwarePointRecord(fmt: <PointFormat(7, 0 bytes of extra dims)>, len: 77486, point size: 36)>
0.021911293026868
Spatial Query gave: <ScaleAwarePointRecord(fmt: <PointFormat(7, 0 bytes of extra dims)>, len: 101208, point size: 36)>
0.028619339553767862
It took: 10.786364793777466  seconds
python copc.py https://s3.amazonaws.com/hobu-lidar/melbourne-2018.copc.laz  2.86s user 2.64s system 50% cpu 10.950 total

[tmontaigu]

OSX arm wheels for lazrs didn't post https://pypi.org/project/lazrs/#files

I have not posted these wheels yet

It seems a bit slower than PDAL while fetching over the network. Are you reading the entire octree down to the client?

Yes the entire octree hierarchy info is fetched at the begnning.
The idea was that when doing multiple queries on the same ressource, we would not have to refetch the octree info.
I'll change it to only fetch progressively the node informations and see how much it improves.

The pdal info only reads the header and does not fetch any octree info nor points ?

The python example reads the whole hierarchy and then queries two quarter of the bounds with a Some resolution

[tmontaigu]

I just pushed changes that makes it so that the octree hierarchies are fetched when actually needed (and not fetch every thing on init like before)

[hobu]

I just pushed changes that makes it so that the octree hierarchies are fetched when actually needed (and not fetch every thing on init like before)

Confirmed and behaving as I would expect now.

I have not posted these wheels yet

I didn't realize these weren't automated. laz-rs/laz-rs-python#9 implied they were, but maybe that isn't quite done yet.

[hobu]

Sorry, my pdal info call was this:

time pdal info --all https://s3.amazonaws.com/hobu-lidar/melbourne-2018.copc.laz --readers.copc.resolution=20.0

That reads across the entire domain and fetches points to a 20m resolution (376303 points)

[tmontaigu]

Confirmed and behaving as I would expect now.

Ok nice

I didn't realize these weren't automated. laz-rs/laz-rs-python#9 implied they were, but maybe that isn't quite done yet.

Yeah its a thing that I want to do, it should not be hard, but CI is not the most interesting thing to do in projects 😅