<div dir="ltr"><div dir="ltr">On Fri, Sep 11, 2020 at 7:10 AM Peder Axensten <<a href="mailto:Peder.Axensten@slu.se">Peder.Axensten@slu.se</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">> On 10 Sep 2020, at 16:27, Andrew Bell <<a href="mailto:andrew.bell.ia@gmail.com" target="_blank">andrew.bell.ia@gmail.com</a>> wrote:<br>
><br>
> On Thu, Sep 10, 2020 at 9:43 AM Peder Axensten <<a href="mailto:Peder.Axensten@slu.se" target="_blank">Peder.Axensten@slu.se</a>> wrote:<br>
> I want to implement a writer that takes a point cloud file and rasterises it (percentiles and other statistics):<br>
> class Raster_metrics : public pdal::Writer;<br>
><br>
> Metadata is generally stored with each stage. If you want the metadata that was read, you have to walk stages back until you find the one you're looking for, call getMetadata(), and then extract any information you're looking for. We don't usually use file metadata during processing, with the exception of spatial reference, which is stored separately. We don't use metadata-provided bounding boxes because 1) they have been known to be wrong 2) PDAL may aggregate and drop points in between the time that data is read and the time that a filter is run.<br>
<br>
That is a serious bummer… Our processing chain is based on the fact that the input files (there are tens of thousands of them that we process one by one) are aligned to each other in a regular grid. This ensures that we produce rasters that will align and have neither gaps nor overlaps.</blockquote><div><br></div><div>If you were to use writers.gdal, you can pass whatever bounds you like to the stage as an option and they will be respected, regardless of the data. You can run pdal info to extract the metadata and pass that as bounds to the writer if you wish. I believe some users do this.<br></div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"> Points that are outside their file’s bbox are outliers (!) and ignored. To us, the file’s bbox as given by the header info is at least as important as the spatial reference. Such an access function could simply return an empty bbox for file formats with no such header info:<br>
table.headerBounds( bbox );<br>
if( bbox.empty() || argDontTrustHeaderBbox )table.calculateBounds( bbox );<br></blockquote><div><br></div><div>PDAL just doesn't work like this. The input may be many files with many processing steps before writing a raster. Your processing model of a single file without data modification before calculating a raster is a specialization.</div><div><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">I would rather not use time to scan through all points for min/max values when the file header already supplies that info.</blockquote><div><br></div><div>It's generally trivial compared to other processing.</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"> In a previous processing stage we make sure that the header info is correct. Recalculating the bbox at every step (scanning through all points an extra time) translates to a rather longer processing time as there are many Tbytes of point files… </blockquote><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<br>
It seems that pdal info somehow gets the header info – but I don’t see/understand how this is implemented in code?<br></blockquote><div><br></div><div>pdal info operates on a single file. It emits data from the file header as stored in metadata.</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">I’m trying really hard to present strong arguments for the usefulness of a pdal::PointViewPtr::headerBounds() member function – am I making progress? :-)<br></blockquote><div><br></div><div>You're looking at a PointView as a representation of a LAS file, and it may be that, but it may well NOT be that. Furthermore, many (most?) file formats don't provide such information. PDAL supports 20+ formats.</div><div><br></div><div>If you're wanting to do special processing using the PDAL code, go for it. It's open-source and you're free to modify it for your own purposes. You could modify the LAS reader to put the bounds from the header in the PointView and then use that in your raster creation code, but that's not something that's going to be generally applicable.</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">> Further I should process the actual points in<br>
> void Raster_metrics::write( pdal::PointViewPtr view ) — for chunks of points.<br>
> Correct?<br>
> I want to be able to process points also in a serialised manner – in what member function should I do that?<br>
><br>
> I don't know what you mean. You can simply loop through all the points in the view and process them one at a time. If you don't need all the points prior to beginning processing, you can implement a stream mode filter by inheriting from Streamable and implement processOne().  That said, if you need bounds, you need all the points, which means that your filter won't work in stream mode.<br>
<br>
Another argument for a pdal::PointViewPtr::headerBounds() member function?<br>
I imagine there are more tools than mine that could be made Streamable using this info?<br></blockquote><div><br></div><div>See above. Also, this information wouldn't be sufficient to permit most of the non-streamable filters to become streaming.</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">> And finally I should calculate the statistics and output the these as rasters [using gdal] in<br>
> void Raster_metrics::done( pdal::PointTableRef table )<br>
> Correct?<br>
><br>
> This is up to you. Using the new functionality in 2.2, you can call view.createRaster(). It will return a raster to which you can write data. I would do this in filter() or run().<br>
<br>
That sounds very useful! I will check it out when I get to implement this part.<br>
<br>
To summarise, your suggestion is that I should do everything (setup, consuming the points, and saving the rasters to file) in one member function, either in filter() or run()? This would certainly simplify things. What are the principle differences between the two? What would be the more natural fit for my Stage?<br></blockquote><div><br></div><div>If you're creating a raster as final output, you should create a Writer and implement the write() function.</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Would it be more natural to implement my tool as a Filter, rather than a Writer?<br>
Being new to pdal I don’t really understand how the principle differences between the two applies to a tool such as mine.<br></blockquote><div><br></div><div>For all practical purposes, a filter and writer are the same, but writers can be created/inferred from the output filenames in pipelines.</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Are filter()/run() called once for each input file?<br></blockquote><div><br></div><div>They are called once for each PointView. Each input file is initially placed in its own PointView. Perhaps this is helpful: <a href="https://pdal.io/development/overview.html">https://pdal.io/development/overview.html</a> </div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">In my case I would probably throw an exception at the second call, as merging files does not really make sense here...<br></blockquote><div><br></div><div>The PDAL "engine" makes these calls as appropriate. You don't need to call them (and in fact can't without modifying the code, as they're private). Your code should call prepare() and execute().</div><div><br></div><div>--</div></div><div dir="ltr" class="gmail_signature">Andrew Bell<br><a href="mailto:andrew.bell.ia@gmail.com" target="_blank">andrew.bell.ia@gmail.com</a></div></div>