[pdal] EPT points in polygon with the C++ API

[Balázs Dukai]

Hi,

I've been testing the C++ API for reading an Entwine Point Tile and I 
have a working implementation for reading, filtering the points in one 
go. I've been relying on the docs, the PDAL unit tests and source code 
for guidance. Each of them are great on its own, but especially thanks 
for the clean and well structured code and tests which makes it 
relatively easy to pick up PDAL.

My goal is to use PDAL as a C++ library to get the points from an EPT 
that are inside a given polygon. There are millions of polygons that the 
application needs to process, thus millions of queries against the EPT. 
The EPT covers a country (AOI), uses LAS and about 2TB in size.

*Background:*

There is a working implementation with LASlib, in which we have several 
LAZ files (tiles if you want) covering the AOI, and run an application 
process for each LAZ file. Due to polygons that span across LAZ files, 
we need to re-tile the LAZ files to include some buffer for each file. 
My assumption is that with the EPT+PDAL combo we could omit the costly 
re-tiling process and just query the EPT directly, and have "good" 
performance in retrieving the points for each polygon.

*What I did:*

Below is some pesudocode of the relevant parts of my current, not 
working, implementation.

polygons = a vector of 2D polygons

dirpath = path to the EPT directory

function ( polygons, dirpath ) {

     pdal::EptReader *reader*;
     pdal::PointTable *eptTable*;
     std::string *eptPath* = "ept://" + dirpath;
     std::vector<pdal::BOX2D> *poly_bboxes*;

     // Compute the bounding box for each polygon and create a 
pdal::BOX2D from it, which will be passed as "bounds" to the EptReader

     for ( *polygon* in *polygons* ) {

*poly_bboxes*.push_back( pdal::BOX2D( compute_bbox(*polygon *) );

     }

     for ( *box2d* in *poly_bboxes* ) {

         // I'm re-creating the Options for each polygon, because I 
think I need to update the "bounds" for each of them in order to get the 
points only withing their bounding box.

         pdal::Options *options*;
*options*.add( "filename", *eptPath* );
*options*.add( "bounds", *box2d* );

         // I'm using setOptions to clear all previously set Options and 
sets the new ones. Essentially, to replace the "bounds" of the reader. 
This is the essence of my question, how to reset the "bounds" on the 
reader in a loop?

*reader*.setOptions(*options* );
*reader*.prepare(*eptTable* );         // I'm suspecting that I should 
do this only once in the beginning (and not in the loop), but the docs 
say that .prepare() must be called at the terminal stage, after all 
options has been set.

         const auto *set*(*reader*.execute(*eptTable *) );

         for (const pdal::PointViewPtr& *view* : *set*) {

             for (pdal::point_count_t p(0); p < *view*->size(); ++p) {

                 // I'm processing the point here. This works fine.

             }

         }

     }

}

*Questons:*

1. My primary, technical question is how to reset the "bounds" Option on 
a Reader in a loop?

2. A more "meta" question is if the process above is sensible and 
efficient in terms of PDAL API use for a retrieving the points within a 
given polygon, for many many polygons :-)?

3. For the time being we have our own point-in-polygon implementation, 
hence I'm only querying the EPT for the bounding boxes of the polygons. 
Once the points in the bbox are retrieved, we do the actual 
point-in-polygon test. But I've been looking into the Crop Filter to 
eventually completely replace our point-in-polygon implementation. As 
far as I understand I would need to feed the Crop Filter with the WKT of 
the polygon, which would have some overhead for casting back-and-forth 
to WKT. Is the Crop Filter a reasonable choice for replacing the 
point-in-polygon query in the process above?

I hope someone could give me a few pointers even though the niche 
question :-)

Thanks,

Balázs