Skip to content

Commit ee471bf

Browse files
alexbruywonder-sk
alexbruy
authored and
wonder-sk
committedApr 13, 2023
update pdal_wrench
·
ltr-3_40final-3_32_0
1 parent c4774d6 commit ee471bf

File tree

16 files changed

+826
-170
lines changed

16 files changed

+826
-170
lines changed
 

‎external/pdal_wrench/alg.cpp

Lines changed: 16 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@
1818
#include <thread>
1919

2020
#include <pdal/QuickInfo.hpp>
21+
#include <pdal/util/Bounds.hpp>
2122

2223
using namespace pdal;
2324

@@ -79,14 +80,15 @@ bool runAlg(std::vector<std::string> args, Alg &alg)
7980

8081

8182
bool Alg::parseArgs(std::vector<std::string> args)
82-
{
83+
{
8384
pdal::Arg* argInput = nullptr;
8485
if (hasSingleInput)
8586
{
8687
argInput = &programArgs.add("input,i", "Input point cloud file", inputFile);
8788
}
8889

8990
(void)programArgs.add("filter,f", "Filter expression for input data", filterExpression);
91+
(void)programArgs.add("bounds", "Filter by rectangle", filterBounds);
9092

9193
addArgs(); // impl in derived
9294

@@ -112,6 +114,19 @@ bool Alg::parseArgs(std::vector<std::string> args)
112114
return false;
113115
}
114116

117+
if (!filterBounds.empty())
118+
{
119+
try
120+
{
121+
parseBounds(filterBounds);
122+
}
123+
catch (pdal::Bounds::error err)
124+
{
125+
std::cerr << "invalid bounds: " << err.what() << std::endl;
126+
return false;
127+
}
128+
}
129+
115130
if (!checkArgs()) // impl in derived class
116131
return false;
117132

‎external/pdal_wrench/alg.hpp

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -42,6 +42,8 @@ struct Alg
4242

4343
std::string filterExpression; // optional argument to limit input points
4444

45+
std::string filterBounds; // optional clipping rectangle for input (pdal::Bounds)
46+
4547
bool needsSingleCrs = true; // most algs assume that all input files in VPC are in the same CRS,
4648
// and only few exceptions (e.g. info) tolerate mixture of multiple CRS
4749

‎external/pdal_wrench/boundary.cpp

Lines changed: 29 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -68,6 +68,32 @@ static std::unique_ptr<PipelineManager> pipeline(ParallelJobInfo *tile, double r
6868

6969
Stage& r = manager->makeReader(tile->inputFilenames[0], "");
7070

71+
Stage *last = &r;
72+
73+
// filtering
74+
if (!tile->filterBounds.empty())
75+
{
76+
Options filter_opts;
77+
filter_opts.add(pdal::Option("bounds", tile->filterBounds));
78+
79+
if (readerSupportsBounds(r))
80+
{
81+
// Reader of the format can do the filtering - use that whenever possible!
82+
r.addOptions(filter_opts);
83+
}
84+
else
85+
{
86+
// Reader can't do the filtering - do it with a filter
87+
last = &manager->makeFilter( "filters.crop", *last, filter_opts);
88+
}
89+
}
90+
if (!tile->filterExpression.empty())
91+
{
92+
Options filter_opts;
93+
filter_opts.add(pdal::Option("expression", tile->filterExpression));
94+
last = &manager->makeFilter( "filters.expression", *last, filter_opts);
95+
}
96+
7197
// TODO: what edge size? (by default samples 5000 points if not specified
7298
// TODO: set threshold ? (default at least 16 points to keep the cell)
7399
// btw. if threshold=0, there are still missing points because of simplification (smooth=True)
@@ -78,13 +104,8 @@ static std::unique_ptr<PipelineManager> pipeline(ParallelJobInfo *tile, double r
78104
hexbin_opts.add(pdal::Option("edge_size", resolution));
79105
}
80106
hexbin_opts.add(pdal::Option("threshold", pointsThreshold));
107+
(void)manager->makeFilter( "filters.hexbin", *last, hexbin_opts );
81108

82-
if (!tile->filterExpression.empty())
83-
{
84-
hexbin_opts.add(pdal::Option("where", tile->filterExpression));
85-
}
86-
87-
(void)manager->makeFilter( "filters.hexbin", r, hexbin_opts );
88109
return manager;
89110
}
90111

@@ -99,14 +120,14 @@ void Boundary::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& p
99120

100121
for (const VirtualPointCloud::File& f : vpc.files)
101122
{
102-
ParallelJobInfo tile(ParallelJobInfo::FileBased, BOX2D(), filterExpression);
123+
ParallelJobInfo tile(ParallelJobInfo::FileBased, BOX2D(), filterExpression, filterBounds);
103124
tile.inputFilenames.push_back(f.filename);
104125
pipelines.push_back(pipeline(&tile, resolution, pointsThreshold));
105126
}
106127
}
107128
else
108129
{
109-
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression);
130+
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression, filterBounds);
110131
tile.inputFilenames.push_back(inputFile);
111132
pipelines.push_back(pipeline(&tile, resolution, pointsThreshold));
112133
}

‎external/pdal_wrench/clip.cpp

Lines changed: 27 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -118,21 +118,38 @@ static std::unique_ptr<PipelineManager> pipeline(ParallelJobInfo *tile, const pd
118118

119119
Stage& r = manager->makeReader( tile->inputFilenames[0], "");
120120

121-
Stage& f = manager->makeFilter( "filters.crop", r, crop_opts );
121+
Stage *last = &r;
122122

123-
pdal::Options writer_opts;
124-
writer_opts.add(pdal::Option("forward", "all"));
125-
126-
Stage& w = manager->makeWriter( tile->outputFilename, "", f, writer_opts);
123+
// filtering
124+
if (!tile->filterBounds.empty())
125+
{
126+
Options filter_opts;
127+
filter_opts.add(pdal::Option("bounds", tile->filterBounds));
127128

129+
if (readerSupportsBounds(r))
130+
{
131+
// Reader of the format can do the filtering - use that whenever possible!
132+
r.addOptions(filter_opts);
133+
}
134+
else
135+
{
136+
// Reader can't do the filtering - do it with a filter
137+
last = &manager->makeFilter( "filters.crop", *last, filter_opts);
138+
}
139+
}
128140
if (!tile->filterExpression.empty())
129141
{
130142
Options filter_opts;
131-
filter_opts.add(pdal::Option("where", tile->filterExpression));
132-
f.addOptions(filter_opts);
133-
w.addOptions(filter_opts);
143+
filter_opts.add(pdal::Option("expression", tile->filterExpression));
144+
last = &manager->makeFilter( "filters.expression", *last, filter_opts);
134145
}
135146

147+
last = &manager->makeFilter( "filters.crop", *last, crop_opts );
148+
149+
pdal::Options writer_opts;
150+
writer_opts.add(pdal::Option("forward", "all"));
151+
manager->makeWriter( tile->outputFilename, "", *last, writer_opts);
152+
136153
return manager;
137154
}
138155

@@ -175,7 +192,7 @@ void Clip::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& pipel
175192
std::cout << "using " << f.filename << std::endl;
176193
}
177194

178-
ParallelJobInfo tile(ParallelJobInfo::FileBased, BOX2D(), filterExpression);
195+
ParallelJobInfo tile(ParallelJobInfo::FileBased, BOX2D(), filterExpression, filterBounds);
179196
tile.inputFilenames.push_back(f.filename);
180197

181198
// for input file /x/y/z.las that goes to /tmp/hello.vpc,
@@ -190,7 +207,7 @@ void Clip::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& pipel
190207
}
191208
else
192209
{
193-
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression);
210+
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression, filterBounds);
194211
tile.inputFilenames.push_back(inputFile);
195212
tile.outputFilename = outputFile;
196213
pipelines.push_back(pipeline(&tile, crop_opts));

‎external/pdal_wrench/density.cpp

Lines changed: 65 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -73,20 +73,53 @@ std::unique_ptr<PipelineManager> Density::pipeline(ParallelJobInfo *tile) const
7373
readers.push_back(&manager->makeReader(f, ""));
7474
}
7575

76-
if (tile->mode == ParallelJobInfo::Spatial)
76+
std::vector<Stage*> last = readers;
77+
78+
// find out what will be the bounding box for this job
79+
// (there could be also no bbox if there's no "bounds" filter and no tiling)
80+
BOX2D filterBox = !tile->filterBounds.empty() ? parseBounds(tile->filterBounds).to2d() : BOX2D();
81+
BOX2D box = intersectTileBoxWithFilterBox(tile->box, filterBox);
82+
83+
if (box.valid())
7784
{
85+
// We are going to do filtering of points based on 2D box. Ideally we want to do
86+
// the filtering in the reader (if the reader can do it efficiently like copc/ept),
87+
// otherwise we have to add filters.crop stage to filter points after they were read
88+
7889
for (Stage* reader : readers)
7990
{
80-
// with COPC files, we can also specify bounds at the reader
81-
// that will only read the required parts of the file
82-
if (reader->getName() == "readers.copc")
91+
if (readerSupportsBounds(*reader))
8392
{
93+
// add "bounds" option to reader
8494
pdal::Options copc_opts;
8595
copc_opts.add(pdal::Option("threads", 1));
86-
copc_opts.add(pdal::Option("bounds", box_to_pdal_bounds(tile->box)));
96+
copc_opts.add(pdal::Option("bounds", box_to_pdal_bounds(box)));
8797
reader->addOptions(copc_opts);
8898
}
8999
}
100+
101+
if (!allReadersSupportBounds(readers) && !tile->filterBounds.empty())
102+
{
103+
// At least some readers can't do the filtering - do it with a filter
104+
Options filter_opts;
105+
filter_opts.add(pdal::Option("bounds", tile->filterBounds));
106+
Stage *filterCrop = &manager->makeFilter( "filters.crop", filter_opts);
107+
for (Stage *s : last)
108+
filterCrop->setInput(*s);
109+
last.clear();
110+
last.push_back(filterCrop);
111+
}
112+
}
113+
114+
if (!tile->filterExpression.empty())
115+
{
116+
Options filter_opts;
117+
filter_opts.add(pdal::Option("expression", tile->filterExpression));
118+
Stage *filterExpr = &manager->makeFilter( "filters.expression", filter_opts);
119+
for (Stage *s : last)
120+
filterExpr->setInput(*s);
121+
last.clear();
122+
last.push_back(filterExpr);
90123
}
91124

92125
pdal::Options writer_opts;
@@ -98,9 +131,9 @@ std::unique_ptr<PipelineManager> Density::pipeline(ParallelJobInfo *tile) const
98131
writer_opts.add(pdal::Option("gdalopts", "TILED=YES"));
99132
writer_opts.add(pdal::Option("gdalopts", "COMPRESS=DEFLATE"));
100133

101-
if (tile->box.valid())
134+
if (box.valid())
102135
{
103-
BOX2D box2 = tile->box;
136+
BOX2D box2 = box;
104137
// fix tile size - PDAL's writers.gdal adds one pixel (see GDALWriter::createGrid()),
105138
// because it probably expects that that the bounds and resolution do not perfectly match
106139
box2.maxx -= resolution;
@@ -109,20 +142,13 @@ std::unique_ptr<PipelineManager> Density::pipeline(ParallelJobInfo *tile) const
109142
writer_opts.add(pdal::Option("bounds", box_to_pdal_bounds(box2)));
110143
}
111144

112-
if (!tile->filterExpression.empty())
113-
{
114-
writer_opts.add(pdal::Option("where", tile->filterExpression));
115-
}
116-
117145
// TODO: "writers.gdal: Requested driver 'COG' does not support file creation.""
118146
// writer_opts.add(pdal::Option("gdaldriver", "COG"));
119147

120148
pdal::StageCreationOptions opts{ tile->outputFilename, "", nullptr, writer_opts, "" };
121149
Stage& w = manager->makeWriter( opts );
122-
for (Stage *stage : readers)
123-
{
124-
w.setInput(*stage); // connect all readers to the writer
125-
}
150+
for (Stage *stage : last)
151+
w.setInput(*stage);
126152

127153
return manager;
128154
}
@@ -172,7 +198,14 @@ void Density::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& pi
172198
// to avoid empty areas in resulting rasters
173199
tileBox.clip(gridBounds);
174200

175-
ParallelJobInfo tile(ParallelJobInfo::Spatial, tileBox, filterExpression);
201+
if (!filterBounds.empty() && !intersectionBox2D(tileBox, parseBounds(filterBounds).to2d()).valid())
202+
{
203+
if (verbose)
204+
std::cout << "skipping tile " << iy << " " << ix << " -- " << tileBox.toBox() << std::endl;
205+
continue;
206+
}
207+
208+
ParallelJobInfo tile(ParallelJobInfo::Spatial, tileBox, filterExpression, filterBounds);
176209
for (const VirtualPointCloud::File & f: vpc.overlappingBox2D(tileBox))
177210
{
178211
tile.inputFilenames.push_back(f.filename);
@@ -226,7 +259,14 @@ void Density::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& pi
226259
{
227260
BOX2D tileBox = t.boxAt(ix, iy);
228261

229-
ParallelJobInfo tile(ParallelJobInfo::Spatial, tileBox, filterExpression);
262+
if (!filterBounds.empty() && !intersectionBox2D(tileBox, parseBounds(filterBounds).to2d()).valid())
263+
{
264+
if (verbose)
265+
std::cout << "skipping tile " << iy << " " << ix << " -- " << tileBox.toBox() << std::endl;
266+
continue;
267+
}
268+
269+
ParallelJobInfo tile(ParallelJobInfo::Spatial, tileBox, filterExpression, filterBounds);
230270
tile.inputFilenames.push_back(inputFile);
231271

232272
// create temp output file names
@@ -245,7 +285,7 @@ void Density::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& pi
245285
{
246286
// single input LAS/LAZ - no parallelism
247287

248-
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression);
288+
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression, filterBounds);
249289
tile.inputFilenames.push_back(inputFile);
250290
tile.outputFilename = outputFile;
251291
pipelines.push_back(pipeline(&tile));
@@ -256,8 +296,13 @@ void Density::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& pi
256296

257297
void Density::finalize(std::vector<std::unique_ptr<PipelineManager>>& pipelines)
258298
{
259-
if (pipelines.size() > 1)
299+
if (!tileOutputFiles.empty())
260300
{
261301
rasterTilesToCog(tileOutputFiles, outputFile);
302+
303+
// clean up the temporary directory
304+
fs::path outputParentDir = fs::path(outputFile).parent_path();
305+
fs::path outputSubdir = outputParentDir / fs::path(outputFile).stem();
306+
fs::remove_all(outputSubdir);
262307
}
263308
}

‎external/pdal_wrench/merge.cpp

Lines changed: 41 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -55,27 +55,59 @@ static std::unique_ptr<PipelineManager> pipeline(ParallelJobInfo *tile)
5555

5656
std::unique_ptr<PipelineManager> manager( new PipelineManager );
5757

58-
pdal::Options options;
59-
options.add(pdal::Option("forward", "all"));
60-
61-
if (!tile->filterExpression.empty())
58+
std::vector<Stage*> readers;
59+
for (const std::string& f : tile->inputFilenames)
6260
{
63-
options.add(pdal::Option("where", tile->filterExpression));
61+
readers.push_back(&manager->makeReader(f, ""));
6462
}
6563

66-
Stage& w = manager->makeWriter(tile->outputFilename, "", options);
64+
std::vector<Stage*> last = readers;
6765

68-
for (const std::string& f : tile->inputFilenames)
66+
// filtering
67+
if (!tile->filterBounds.empty())
6968
{
70-
w.setInput(manager->makeReader(f, ""));
69+
Options filter_opts;
70+
filter_opts.add(pdal::Option("bounds", tile->filterBounds));
71+
72+
if (allReadersSupportBounds(readers))
73+
{
74+
// Reader of the format can do the filtering - use that whenever possible!
75+
for (Stage *r : readers)
76+
r->addOptions(filter_opts);
77+
}
78+
else
79+
{
80+
// Reader can't do the filtering - do it with a filter
81+
Stage *filterCrop = &manager->makeFilter( "filters.crop", filter_opts);
82+
for (Stage *s : last)
83+
filterCrop->setInput(*s);
84+
last.clear();
85+
last.push_back(filterCrop);
86+
}
7187
}
88+
if (!tile->filterExpression.empty())
89+
{
90+
Options filter_opts;
91+
filter_opts.add(pdal::Option("expression", tile->filterExpression));
92+
Stage *filterExpr = &manager->makeFilter( "filters.expression", filter_opts);
93+
for (Stage *s : last)
94+
filterExpr->setInput(*s);
95+
last.clear();
96+
last.push_back(filterExpr);
97+
}
98+
99+
pdal::Options options;
100+
options.add(pdal::Option("forward", "all"));
101+
Stage* writer = &manager->makeWriter(tile->outputFilename, "", options);
102+
for (Stage *s : last)
103+
writer->setInput(*s);
72104

73105
return manager;
74106
}
75107

76108
void Merge::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& pipelines, const BOX3D &, point_count_t &totalPoints)
77109
{
78-
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression);
110+
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression, filterBounds);
79111
tile.inputFilenames = inputFiles;
80112
tile.outputFilename = outputFile;
81113

‎external/pdal_wrench/thin.cpp

Lines changed: 30 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -99,33 +99,48 @@ static std::unique_ptr<PipelineManager> pipeline(ParallelJobInfo *tile, std::str
9999

100100
Stage& r = manager->makeReader( tile->inputFilenames[0], "");
101101

102-
Stage *filterPtr = nullptr;
102+
Stage *last = &r;
103+
104+
// filtering
105+
if (!tile->filterBounds.empty())
106+
{
107+
Options filter_opts;
108+
filter_opts.add(pdal::Option("bounds", tile->filterBounds));
109+
110+
if (readerSupportsBounds(r))
111+
{
112+
// Reader of the format can do the filtering - use that whenever possible!
113+
r.addOptions(filter_opts);
114+
}
115+
else
116+
{
117+
// Reader can't do the filtering - do it with a filter
118+
last = &manager->makeFilter( "filters.crop", *last, filter_opts);
119+
}
120+
}
121+
if (!tile->filterExpression.empty())
122+
{
123+
Options filter_opts;
124+
filter_opts.add(pdal::Option("expression", tile->filterExpression));
125+
last = &manager->makeFilter( "filters.expression", *last, filter_opts);
126+
}
103127

104128
if (mode == "every-nth")
105129
{
106130
pdal::Options decim_opts;
107131
decim_opts.add(pdal::Option("step", stepEveryN));
108-
filterPtr = &manager->makeFilter( "filters.decimation", r, decim_opts );
132+
last = &manager->makeFilter( "filters.decimation", *last, decim_opts );
109133
}
110134
else if (mode == "sample")
111135
{
112136
pdal::Options sample_opts;
113137
sample_opts.add(pdal::Option("cell", stepSample));
114-
filterPtr = &manager->makeFilter( "filters.sample", r, sample_opts );
138+
last = &manager->makeFilter( "filters.sample", *last, sample_opts );
115139
}
116140

117141
pdal::Options writer_opts;
118142
writer_opts.add(pdal::Option("forward", "all")); // TODO: maybe we could use lower scale than the original
119-
120-
Stage& w = manager->makeWriter( tile->outputFilename, "", *filterPtr, writer_opts);
121-
122-
if (!tile->filterExpression.empty())
123-
{
124-
Options filter_opts;
125-
filter_opts.add(pdal::Option("where", tile->filterExpression));
126-
filterPtr->addOptions(filter_opts);
127-
w.addOptions(filter_opts);
128-
}
143+
manager->makeWriter( tile->outputFilename, "", *last, writer_opts);
129144

130145
return manager;
131146
}
@@ -153,7 +168,7 @@ void Thin::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& pipel
153168

154169
for (const VirtualPointCloud::File& f : vpc.files)
155170
{
156-
ParallelJobInfo tile(ParallelJobInfo::FileBased, BOX2D(), filterExpression);
171+
ParallelJobInfo tile(ParallelJobInfo::FileBased, BOX2D(), filterExpression, filterBounds);
157172
tile.inputFilenames.push_back(f.filename);
158173

159174
// for input file /x/y/z.las that goes to /tmp/hello.vpc,
@@ -168,7 +183,7 @@ void Thin::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& pipel
168183
}
169184
else
170185
{
171-
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression);
186+
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression, filterBounds);
172187
tile.inputFilenames.push_back(inputFile);
173188
tile.outputFilename = outputFile;
174189
pipelines.push_back(pipeline(&tile, mode, stepEveryN, stepSample));

‎external/pdal_wrench/tile/tile.cpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -684,7 +684,7 @@ static void tilingPass2(BaseInfo &m_b, TileGrid &m_grid, FileInfo &m_srsFileInfo
684684
{
685685
size_t nPointsToRead = std::min(readChunkSize, ptCnt-i);
686686
fileReader.read(contentPtr, nPointsToRead*m_b.pointSize);
687-
for (size_t a = 0; a < nPointsToRead; ++a)
687+
for (std::size_t a = 0; a < nPointsToRead; ++a)
688688
{
689689
char *base = contentPtr + a * m_b.pointSize;
690690
for (const untwine::FileDimInfo& fdi : dims)

‎external/pdal_wrench/to_raster.cpp

Lines changed: 87 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -74,7 +74,9 @@ bool ToRaster::checkArgs()
7474
}
7575

7676

77-
static std::unique_ptr<PipelineManager> pipeline(ParallelJobInfo *tile, double resolution, std::string attribute)
77+
78+
79+
static std::unique_ptr<PipelineManager> pipeline(ParallelJobInfo *tile, double resolution, std::string attribute, double collarSize)
7880
{
7981
std::unique_ptr<PipelineManager> manager( new PipelineManager );
8082

@@ -84,20 +86,64 @@ static std::unique_ptr<PipelineManager> pipeline(ParallelJobInfo *tile, double r
8486
readers.push_back(&manager->makeReader(f, ""));
8587
}
8688

87-
if (tile->mode == ParallelJobInfo::Spatial)
89+
std::vector<Stage*> last = readers;
90+
91+
// find out what will be the bounding box for this job
92+
// (there could be also no bbox if there's no "bounds" filter and no tiling)
93+
BOX2D filterBox = !tile->filterBounds.empty() ? parseBounds(tile->filterBounds).to2d() : BOX2D();
94+
BOX2D box = intersectTileBoxWithFilterBox(tile->box, filterBox);
95+
BOX2D boxWithCollar;
96+
97+
// box with collar is used for filtering of data on the input
98+
// bow without collar is used for output bounds
99+
100+
if (box.valid())
88101
{
102+
BOX2D filterBoxWithCollar = filterBox;
103+
if (filterBoxWithCollar.valid())
104+
filterBoxWithCollar.grow(collarSize);
105+
boxWithCollar = tile->box;
106+
boxWithCollar.grow(collarSize);
107+
boxWithCollar = intersectTileBoxWithFilterBox(boxWithCollar, filterBoxWithCollar);
108+
109+
// We are going to do filtering of points based on 2D box. Ideally we want to do
110+
// the filtering in the reader (if the reader can do it efficiently like copc/ept),
111+
// otherwise we have to add filters.crop stage to filter points after they were read
112+
89113
for (Stage* reader : readers)
90114
{
91-
// with COPC files, we can also specify bounds at the reader
92-
// that will only read the required parts of the file
93-
if (reader->getName() == "readers.copc")
115+
if (readerSupportsBounds(*reader))
94116
{
117+
// add "bounds" option to reader
95118
pdal::Options copc_opts;
96119
copc_opts.add(pdal::Option("threads", 1));
97-
copc_opts.add(pdal::Option("bounds", box_to_pdal_bounds(tile->boxWithCollar)));
120+
copc_opts.add(pdal::Option("bounds", box_to_pdal_bounds(boxWithCollar)));
98121
reader->addOptions(copc_opts);
99122
}
100123
}
124+
125+
if (!allReadersSupportBounds(readers) && !tile->filterBounds.empty())
126+
{
127+
// At least some readers can't do the filtering - do it with a filter
128+
Options filter_opts;
129+
filter_opts.add(pdal::Option("bounds", box_to_pdal_bounds(filterBoxWithCollar)));
130+
Stage *filterCrop = &manager->makeFilter( "filters.crop", filter_opts);
131+
for (Stage *s : last)
132+
filterCrop->setInput(*s);
133+
last.clear();
134+
last.push_back(filterCrop);
135+
}
136+
}
137+
138+
if (!tile->filterExpression.empty())
139+
{
140+
Options filter_opts;
141+
filter_opts.add(pdal::Option("expression", tile->filterExpression));
142+
Stage *filterExpr = &manager->makeFilter( "filters.expression", filter_opts);
143+
for (Stage *s : last)
144+
filterExpr->setInput(*s);
145+
last.clear();
146+
last.push_back(filterExpr);
101147
}
102148

103149
pdal::Options writer_opts;
@@ -109,9 +155,9 @@ static std::unique_ptr<PipelineManager> pipeline(ParallelJobInfo *tile, double r
109155
writer_opts.add(pdal::Option("gdalopts", "TILED=YES"));
110156
writer_opts.add(pdal::Option("gdalopts", "COMPRESS=DEFLATE"));
111157

112-
if (tile->box.valid())
158+
if (box.valid())
113159
{
114-
BOX2D box2 = tile->box;
160+
BOX2D box2 = box;
115161
// fix tile size - PDAL's writers.gdal adds one pixel (see GDALWriter::createGrid()),
116162
// because it probably expects that that the bounds and resolution do not perfectly match
117163
box2.maxx -= resolution;
@@ -120,20 +166,13 @@ static std::unique_ptr<PipelineManager> pipeline(ParallelJobInfo *tile, double r
120166
writer_opts.add(pdal::Option("bounds", box_to_pdal_bounds(box2)));
121167
}
122168

123-
if (!tile->filterExpression.empty())
124-
{
125-
writer_opts.add(pdal::Option("where", tile->filterExpression));
126-
}
127-
128169
// TODO: "writers.gdal: Requested driver 'COG' does not support file creation.""
129170
// writer_opts.add(pdal::Option("gdaldriver", "COG"));
130171

131172
pdal::StageCreationOptions opts{ tile->outputFilename, "", nullptr, writer_opts, "" };
132173
Stage& w = manager->makeWriter( opts );
133-
for (Stage *stage : readers)
134-
{
135-
w.setInput(*stage); // connect all readers to the writer
136-
}
174+
for (Stage *stage : last)
175+
w.setInput(*stage);
137176

138177
return manager;
139178
}
@@ -184,13 +223,20 @@ void ToRaster::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& p
184223
// to avoid empty areas in resulting rasters
185224
tileBox.clip(gridBounds);
186225

187-
ParallelJobInfo tile(ParallelJobInfo::Spatial, tileBox, filterExpression);
226+
if (!filterBounds.empty() && !intersectionBox2D(tileBox, parseBounds(filterBounds).to2d()).valid())
227+
{
228+
if (verbose)
229+
std::cout << "skipping tile " << iy << " " << ix << " -- " << tileBox.toBox() << std::endl;
230+
continue;
231+
}
232+
233+
ParallelJobInfo tile(ParallelJobInfo::Spatial, tileBox, filterExpression, filterBounds);
188234

189235
// add collar to avoid edge effects
190-
tile.boxWithCollar = tileBox;
191-
tile.boxWithCollar.grow(collarSize);
236+
BOX2D boxWithCollar = tileBox;
237+
boxWithCollar.grow(collarSize);
192238

193-
for (const VirtualPointCloud::File & f: vpc.overlappingBox2D(tile.boxWithCollar))
239+
for (const VirtualPointCloud::File & f: vpc.overlappingBox2D(boxWithCollar))
194240
{
195241
tile.inputFilenames.push_back(f.filename);
196242
totalPoints += f.count;
@@ -206,7 +252,7 @@ void ToRaster::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& p
206252

207253
tileOutputFiles.push_back(tile.outputFilename);
208254

209-
pipelines.push_back(pipeline(&tile, resolution, attribute));
255+
pipelines.push_back(pipeline(&tile, resolution, attribute, collarSize));
210256
}
211257
}
212258
}
@@ -232,12 +278,19 @@ void ToRaster::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& p
232278
{
233279
BOX2D tileBox = t.boxAt(ix, iy);
234280

235-
ParallelJobInfo tile(ParallelJobInfo::Spatial, tileBox, filterExpression);
281+
if (!filterBounds.empty() && !intersectionBox2D(tileBox, parseBounds(filterBounds).to2d()).valid())
282+
{
283+
if (verbose)
284+
std::cout << "skipping tile " << iy << " " << ix << " -- " << tileBox.toBox() << std::endl;
285+
continue;
286+
}
287+
288+
ParallelJobInfo tile(ParallelJobInfo::Spatial, tileBox, filterExpression, filterBounds);
236289
tile.inputFilenames.push_back(inputFile);
237290

238291
// add collar to avoid edge effects
239-
tile.boxWithCollar = tileBox;
240-
tile.boxWithCollar.grow(collarSize);
292+
BOX2D boxWithCollar = tileBox;
293+
boxWithCollar.grow(collarSize);
241294

242295
// create temp output file names
243296
// for tile (x=2,y=3) that goes to /tmp/hello.tif,
@@ -247,27 +300,32 @@ void ToRaster::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& p
247300

248301
tileOutputFiles.push_back(tile.outputFilename);
249302

250-
pipelines.push_back(pipeline(&tile, resolution, attribute));
303+
pipelines.push_back(pipeline(&tile, resolution, attribute, collarSize));
251304
}
252305
}
253306
}
254307
else
255308
{
256309
// single input LAS/LAZ - no parallelism
257310

258-
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression);
311+
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression, filterBounds);
259312
tile.inputFilenames.push_back(inputFile);
260313
tile.outputFilename = outputFile;
261-
pipelines.push_back(pipeline(&tile, resolution, attribute));
314+
pipelines.push_back(pipeline(&tile, resolution, attribute, 0));
262315
}
263316

264317
}
265318

266319

267320
void ToRaster::finalize(std::vector<std::unique_ptr<PipelineManager>>& pipelines)
268321
{
269-
if (pipelines.size() > 1)
322+
if (!tileOutputFiles.empty())
270323
{
271324
rasterTilesToCog(tileOutputFiles, outputFile);
325+
326+
// clean up the temporary directory
327+
fs::path outputParentDir = fs::path(outputFile).parent_path();
328+
fs::path outputSubdir = outputParentDir / fs::path(outputFile).stem();
329+
fs::remove_all(outputSubdir);
272330
}
273331
}

‎external/pdal_wrench/to_raster_tin.cpp

Lines changed: 95 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -70,7 +70,7 @@ bool ToRasterTin::checkArgs()
7070
}
7171

7272

73-
std::unique_ptr<PipelineManager> pipeline(ParallelJobInfo *tile, double resolution)
73+
std::unique_ptr<PipelineManager> pipeline(ParallelJobInfo *tile, double resolution, double collarSize)
7474
{
7575
std::unique_ptr<PipelineManager> manager( new PipelineManager );
7676

@@ -80,44 +80,90 @@ std::unique_ptr<PipelineManager> pipeline(ParallelJobInfo *tile, double resoluti
8080
readers.push_back(&manager->makeReader(f, ""));
8181
}
8282

83-
if (tile->mode == ParallelJobInfo::Spatial)
83+
std::vector<Stage*> last = readers;
84+
85+
// find out what will be the bounding box for this job
86+
// (there could be also no bbox if there's no "bounds" filter and no tiling)
87+
BOX2D filterBox = !tile->filterBounds.empty() ? parseBounds(tile->filterBounds).to2d() : BOX2D();
88+
BOX2D box = intersectTileBoxWithFilterBox(tile->box, filterBox);
89+
BOX2D boxWithCollar;
90+
91+
// box with collar is used for filtering of data on the input
92+
// bow without collar is used for output bounds
93+
94+
if (box.valid())
8495
{
96+
BOX2D filterBoxWithCollar = filterBox;
97+
if (filterBoxWithCollar.valid())
98+
filterBoxWithCollar.grow(collarSize);
99+
boxWithCollar = tile->box;
100+
boxWithCollar.grow(collarSize);
101+
boxWithCollar = intersectTileBoxWithFilterBox(boxWithCollar, filterBoxWithCollar);
102+
103+
// We are going to do filtering of points based on 2D box. Ideally we want to do
104+
// the filtering in the reader (if the reader can do it efficiently like copc/ept),
105+
// otherwise we have to add filters.crop stage to filter points after they were read
106+
85107
for (Stage* reader : readers)
86108
{
87-
// with COPC files, we can also specify bounds at the reader
88-
// that will only read the required parts of the file
89-
if (reader->getName() == "readers.copc")
109+
if (readerSupportsBounds(*reader))
90110
{
111+
// add "bounds" option to reader
91112
pdal::Options copc_opts;
92113
copc_opts.add(pdal::Option("threads", 1));
93-
copc_opts.add(pdal::Option("bounds", box_to_pdal_bounds(tile->boxWithCollar)));
114+
copc_opts.add(pdal::Option("bounds", box_to_pdal_bounds(boxWithCollar)));
94115
reader->addOptions(copc_opts);
95116
}
96117
}
97-
}
98118

99-
Stage &delaunay = manager->makeFilter("filters.delaunay");
100-
for (Stage *stage : readers)
101-
{
102-
delaunay.setInput(*stage); // connect all readers to the writer
119+
if (!allReadersSupportBounds(readers) && !tile->filterBounds.empty())
120+
{
121+
// At least some readers can't do the filtering - do it with a filter
122+
Options filter_opts;
123+
filter_opts.add(pdal::Option("bounds", box_to_pdal_bounds(filterBoxWithCollar)));
124+
Stage *filterCrop = &manager->makeFilter( "filters.crop", filter_opts);
125+
for (Stage *s : last)
126+
filterCrop->setInput(*s);
127+
last.clear();
128+
last.push_back(filterCrop);
129+
}
103130
}
104131

105132
if (!tile->filterExpression.empty())
106133
{
107134
Options filter_opts;
108-
filter_opts.add(pdal::Option("where", tile->filterExpression));
109-
delaunay.addOptions(filter_opts);
135+
filter_opts.add(pdal::Option("expression", tile->filterExpression));
136+
Stage *filterExpr = &manager->makeFilter( "filters.expression", filter_opts);
137+
for (Stage *s : last)
138+
filterExpr->setInput(*s);
139+
last.clear();
140+
last.push_back(filterExpr);
141+
}
142+
143+
if (readers.size() > 1)
144+
{
145+
// explicitly merge if there are multiple inputs, otherwise things don't work
146+
// (no pixels are written - not sure which downstream stage is to blame)
147+
Stage *merge = &manager->makeFilter("filters.merge");
148+
for (Stage *stage : last)
149+
merge->setInput(*stage);
150+
last.clear();
151+
last.push_back(merge);
110152
}
111153

154+
Stage &delaunay = manager->makeFilter("filters.delaunay");
155+
for (Stage *stage : last)
156+
delaunay.setInput(*stage);
157+
112158
pdal::Options faceRaster_opts;
113159
faceRaster_opts.add(pdal::Option("resolution", resolution));
114160

115-
if (tile->box.valid()) // if box is not provided, filters.faceraster will calculate it from data
161+
if (box.valid()) // if box is not provided, filters.faceraster will calculate it from data
116162
{
117-
faceRaster_opts.add(pdal::Option("origin_x", tile->box.minx));
118-
faceRaster_opts.add(pdal::Option("origin_y", tile->box.miny));
119-
faceRaster_opts.add(pdal::Option("width", (tile->box.maxx-tile->box.minx)/resolution));
120-
faceRaster_opts.add(pdal::Option("height", (tile->box.maxy-tile->box.miny)/resolution));
163+
faceRaster_opts.add(pdal::Option("origin_x", box.minx));
164+
faceRaster_opts.add(pdal::Option("origin_y", box.miny));
165+
faceRaster_opts.add(pdal::Option("width", ceil((box.maxx-box.minx)/resolution)));
166+
faceRaster_opts.add(pdal::Option("height", ceil((box.maxy-box.miny)/resolution)));
121167
}
122168

123169
Stage &faceRaster = manager->makeFilter("filters.faceraster", delaunay, faceRaster_opts);
@@ -177,13 +223,20 @@ void ToRasterTin::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>
177223
// to avoid empty areas in resulting rasters
178224
tileBox.clip(gridBounds);
179225

180-
ParallelJobInfo tile(ParallelJobInfo::Spatial, tileBox, filterExpression);
226+
if (!filterBounds.empty() && !intersectionBox2D(tileBox, parseBounds(filterBounds).to2d()).valid())
227+
{
228+
if (verbose)
229+
std::cout << "skipping tile " << iy << " " << ix << " -- " << tileBox.toBox() << std::endl;
230+
continue;
231+
}
232+
233+
ParallelJobInfo tile(ParallelJobInfo::Spatial, tileBox, filterExpression, filterBounds);
181234

182235
// add collar to avoid edge effects
183-
tile.boxWithCollar = tileBox;
184-
tile.boxWithCollar.grow(collarSize);
236+
BOX2D boxWithCollar = tileBox;
237+
boxWithCollar.grow(collarSize);
185238

186-
for (const VirtualPointCloud::File & f: vpc.overlappingBox2D(tile.boxWithCollar))
239+
for (const VirtualPointCloud::File & f: vpc.overlappingBox2D(boxWithCollar))
187240
{
188241
tile.inputFilenames.push_back(f.filename);
189242
totalPoints += f.count;
@@ -199,7 +252,7 @@ void ToRasterTin::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>
199252

200253
tileOutputFiles.push_back(tile.outputFilename);
201254

202-
pipelines.push_back(pipeline(&tile, resolution));
255+
pipelines.push_back(pipeline(&tile, resolution, collarSize));
203256
}
204257
}
205258
}
@@ -225,12 +278,19 @@ void ToRasterTin::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>
225278
{
226279
BOX2D tileBox = t.boxAt(ix, iy);
227280

228-
ParallelJobInfo tile(ParallelJobInfo::Spatial, tileBox, filterExpression);
281+
if (!filterBounds.empty() && !intersectionBox2D(tileBox, parseBounds(filterBounds).to2d()).valid())
282+
{
283+
if (verbose)
284+
std::cout << "skipping tile " << iy << " " << ix << " -- " << tileBox.toBox() << std::endl;
285+
continue;
286+
}
287+
288+
ParallelJobInfo tile(ParallelJobInfo::Spatial, tileBox, filterExpression, filterBounds);
229289
tile.inputFilenames.push_back(inputFile);
230290

231291
// add collar to avoid edge effects
232-
tile.boxWithCollar = tileBox;
233-
tile.boxWithCollar.grow(collarSize);
292+
BOX2D boxWithCollar = tileBox;
293+
boxWithCollar.grow(collarSize);
234294

235295
// create temp output file names
236296
// for tile (x=2,y=3) that goes to /tmp/hello.tif,
@@ -240,23 +300,28 @@ void ToRasterTin::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>
240300

241301
tileOutputFiles.push_back(tile.outputFilename);
242302

243-
pipelines.push_back(pipeline(&tile, resolution));
303+
pipelines.push_back(pipeline(&tile, resolution, collarSize));
244304
}
245305
}
246306
}
247307
else
248308
{
249-
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression);
309+
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression, filterBounds);
250310
tile.inputFilenames.push_back(inputFile);
251311
tile.outputFilename = outputFile;
252-
pipelines.push_back(pipeline(&tile, resolution));
312+
pipelines.push_back(pipeline(&tile, resolution, 0));
253313
}
254314
}
255315

256316
void ToRasterTin::finalize(std::vector<std::unique_ptr<PipelineManager>>& pipelines)
257317
{
258-
if (pipelines.size() > 1)
318+
if (!tileOutputFiles.empty())
259319
{
260320
rasterTilesToCog(tileOutputFiles, outputFile);
321+
322+
// clean up the temporary directory
323+
fs::path outputParentDir = fs::path(outputFile).parent_path();
324+
fs::path outputSubdir = outputParentDir / fs::path(outputFile).stem();
325+
fs::remove_all(outputSubdir);
261326
}
262327
}

‎external/pdal_wrench/to_vector.cpp

Lines changed: 28 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -55,17 +55,37 @@ static std::unique_ptr<PipelineManager> pipeline(ParallelJobInfo *tile, const st
5555

5656
Stage& r = manager->makeReader( tile->inputFilenames[0], "");
5757

58-
pdal::Options writer_opts;
59-
writer_opts.add(pdal::Option("ogrdriver", "GPKG"));
60-
if (!attributes.empty())
61-
writer_opts.add(pdal::Option("attr_dims", join_strings(attributes, ',')));
58+
Stage *last = &r;
6259

60+
// filtering
61+
if (!tile->filterBounds.empty())
62+
{
63+
Options filter_opts;
64+
filter_opts.add(pdal::Option("bounds", tile->filterBounds));
65+
66+
if (readerSupportsBounds(r))
67+
{
68+
// Reader of the format can do the filtering - use that whenever possible!
69+
r.addOptions(filter_opts);
70+
}
71+
else
72+
{
73+
// Reader can't do the filtering - do it with a filter
74+
last = &manager->makeFilter( "filters.crop", *last, filter_opts);
75+
}
76+
}
6377
if (!tile->filterExpression.empty())
6478
{
65-
writer_opts.add(pdal::Option("where", tile->filterExpression));
79+
Options filter_opts;
80+
filter_opts.add(pdal::Option("expression", tile->filterExpression));
81+
last = &manager->makeFilter( "filters.expression", *last, filter_opts);
6682
}
6783

68-
(void)manager->makeWriter( tile->outputFilename, "writers.ogr", r, writer_opts);
84+
pdal::Options writer_opts;
85+
writer_opts.add(pdal::Option("ogrdriver", "GPKG"));
86+
if (!attributes.empty())
87+
writer_opts.add(pdal::Option("attr_dims", join_strings(attributes, ',')));
88+
(void)manager->makeWriter( tile->outputFilename, "writers.ogr", *last, writer_opts);
6989

7090
return manager;
7191
}
@@ -87,7 +107,7 @@ void ToVector::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& p
87107

88108
for (const VirtualPointCloud::File& f : vpc.files)
89109
{
90-
ParallelJobInfo tile(ParallelJobInfo::FileBased, BOX2D(), filterExpression);
110+
ParallelJobInfo tile(ParallelJobInfo::FileBased, BOX2D(), filterExpression, filterBounds);
91111
tile.inputFilenames.push_back(f.filename);
92112

93113
// for input file /x/y/z.las that goes to /tmp/hello.vpc,
@@ -102,7 +122,7 @@ void ToVector::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& p
102122
}
103123
else
104124
{
105-
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression);
125+
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression, filterBounds);
106126
tile.inputFilenames.push_back(inputFile);
107127
tile.outputFilename = outputFile;
108128
pipelines.push_back(pipeline(&tile, attributes));

‎external/pdal_wrench/translate.cpp

Lines changed: 32 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -82,6 +82,32 @@ static std::unique_ptr<PipelineManager> pipeline(ParallelJobInfo *tile, std::str
8282

8383
Stage& r = manager->makeReader( tile->inputFilenames[0], "", reader_opts);
8484

85+
Stage *last = &r;
86+
87+
// filtering
88+
if (!tile->filterBounds.empty())
89+
{
90+
Options filter_opts;
91+
filter_opts.add(pdal::Option("bounds", tile->filterBounds));
92+
93+
if (readerSupportsBounds(r))
94+
{
95+
// Reader of the format can do the filtering - use that whenever possible!
96+
r.addOptions(filter_opts);
97+
}
98+
else
99+
{
100+
// Reader can't do the filtering - do it with a filter
101+
last = &manager->makeFilter( "filters.crop", *last, filter_opts);
102+
}
103+
}
104+
if (!tile->filterExpression.empty())
105+
{
106+
Options filter_opts;
107+
filter_opts.add(pdal::Option("expression", tile->filterExpression));
108+
last = &manager->makeFilter( "filters.expression", *last, filter_opts);
109+
}
110+
85111
// optional reprojection
86112
Stage* reproject = nullptr;
87113
if (!transformCrs.empty())
@@ -91,12 +117,13 @@ static std::unique_ptr<PipelineManager> pipeline(ParallelJobInfo *tile, std::str
91117
if (!transformCoordOp.empty())
92118
{
93119
transform_opts.add(pdal::Option("coord_op", transformCoordOp));
94-
reproject = &manager->makeFilter( "filters.projpipeline", r, transform_opts);
120+
reproject = &manager->makeFilter( "filters.projpipeline", *last, transform_opts);
95121
}
96122
else
97123
{
98-
reproject = &manager->makeFilter( "filters.reprojection", r, transform_opts);
124+
reproject = &manager->makeFilter( "filters.reprojection", *last, transform_opts);
99125
}
126+
last = reproject;
100127
}
101128

102129
pdal::Options writer_opts;
@@ -115,17 +142,7 @@ static std::unique_ptr<PipelineManager> pipeline(ParallelJobInfo *tile, std::str
115142
writer_opts.add(pdal::Option("offset_z", "auto"));
116143
}
117144

118-
Stage& writerInput = reproject ? *reproject : r;
119-
Stage& w = manager->makeWriter( tile->outputFilename, "", writerInput, writer_opts);
120-
121-
if (!tile->filterExpression.empty())
122-
{
123-
Options filter_opts;
124-
filter_opts.add(pdal::Option("where", tile->filterExpression));
125-
w.addOptions(filter_opts);
126-
if (reproject)
127-
reproject->addOptions(filter_opts);
128-
}
145+
Stage& w = manager->makeWriter( tile->outputFilename, "", *last, writer_opts);
129146

130147
return manager;
131148
}
@@ -153,7 +170,7 @@ void Translate::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>&
153170

154171
for (const VirtualPointCloud::File& f : vpc.files)
155172
{
156-
ParallelJobInfo tile(ParallelJobInfo::FileBased, BOX2D(), filterExpression);
173+
ParallelJobInfo tile(ParallelJobInfo::FileBased, BOX2D(), filterExpression, filterBounds);
157174
tile.inputFilenames.push_back(f.filename);
158175

159176
// for input file /x/y/z.las that goes to /tmp/hello.vpc,
@@ -168,7 +185,7 @@ void Translate::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>&
168185
}
169186
else
170187
{
171-
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression);
188+
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression, filterBounds);
172189
tile.inputFilenames.push_back(inputFile);
173190
tile.outputFilename = outputFile;
174191
pipelines.push_back(pipeline(&tile, assignCrs, transformCrs, transformCoordOp));

‎external/pdal_wrench/utils.cpp

Lines changed: 63 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@
1212

1313
#include "utils.hpp"
1414

15+
#include <filesystem>
1516
#include <iostream>
1617
#include <chrono>
1718

@@ -196,7 +197,68 @@ bool rasterTilesToCog(const std::vector<std::string> &inputFiles, const std::str
196197
rasterVrtToCog(ds, outputFile);
197198
GDALClose(ds);
198199

199-
// TODO: remove VRT + partial tifs after gdal_translate?
200+
std::filesystem::remove(outputVrt);
200201

201202
return true;
202203
}
204+
205+
bool readerSupportsBounds(Stage &reader)
206+
{
207+
// these readers support "bounds" option with a 2D/3D bounding box, and based
208+
// on it, they will do very efficient reading of data and only return what's
209+
// in the given bounding box
210+
return reader.getName() == "readers.copc" || reader.getName() == "readers.ept";
211+
}
212+
213+
bool allReadersSupportBounds(const std::vector<Stage *> &readers)
214+
{
215+
for (Stage *r : readers)
216+
{
217+
if (!readerSupportsBounds(*r))
218+
return false;
219+
}
220+
return true;
221+
}
222+
223+
pdal::Bounds parseBounds(const std::string &boundsStr)
224+
{
225+
// if the input string is not a correct 2D/3D PDAL bounds then parse()
226+
// will throw an exception
227+
pdal::Bounds b;
228+
std::string::size_type pos(0);
229+
b.parse(boundsStr, pos);
230+
return b;
231+
}
232+
233+
BOX2D intersectionBox2D(const BOX2D &b1, const BOX2D &b2)
234+
{
235+
BOX2D b;
236+
b.minx = b1.minx > b2.minx ? b1.minx : b2.minx;
237+
b.miny = b1.miny > b2.miny ? b1.miny : b2.miny;
238+
b.maxx = b1.maxx < b2.maxx ? b1.maxx : b2.maxx;
239+
b.maxy = b1.maxy < b2.maxy ? b1.maxy : b2.maxy;
240+
if (b.minx > b.maxx || b.miny > b.maxy)
241+
return BOX2D();
242+
return b;
243+
}
244+
245+
246+
BOX2D intersectTileBoxWithFilterBox(const BOX2D &tileBox, const BOX2D &filterBox)
247+
{
248+
if (tileBox.valid() && filterBox.valid())
249+
{
250+
return intersectionBox2D(tileBox, filterBox);
251+
}
252+
else if (tileBox.valid())
253+
{
254+
return tileBox;
255+
}
256+
else if (filterBox.valid())
257+
{
258+
return filterBox;
259+
}
260+
else
261+
{
262+
return BOX2D(); // invalid box
263+
}
264+
}

‎external/pdal_wrench/utils.hpp

Lines changed: 17 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -77,7 +77,8 @@ struct ParallelJobInfo
7777
} mode;
7878

7979
ParallelJobInfo(ParallelMode m = Single): mode(m) {}
80-
ParallelJobInfo(ParallelMode m, const BOX2D &b, const std::string fe): mode(m), box(b), filterExpression(fe) {}
80+
ParallelJobInfo(ParallelMode m, const BOX2D &b, const std::string fe, const std::string fb)
81+
: mode(m), box(b), filterExpression(fe), filterBounds(fb) {}
8182

8283
// what input point cloud files to read for a job
8384
std::vector<std::string> inputFilenames;
@@ -88,20 +89,20 @@ struct ParallelJobInfo
8889
// bounding box for this job (for input/output)
8990
BOX2D box;
9091

91-
// bounding box for the job with extra collar that some algs may use
92-
// in case they need access to neighboring points at the edges of tiles
93-
BOX2D boxWithCollar;
94-
9592
// PDAL filter expression to apply on all pipelines
9693
std::string filterExpression;
9794

95+
// PDAL filter on 2D or 3D bounds to apply on all pipelines
96+
// Format is "([xmin, xmax], [ymin, ymax])" or "([xmin, xmax], [ymin, ymax], [zmin, zmax])"
97+
std::string filterBounds;
98+
9899
// modes of operation:
99100
// A. multi input without box (LAS/LAZ) -- per file strategy
100101
// - all input files are processed, no filtering on bounding box
101102
// B. multi input with box (anything) -- tile strategy
102103
// - all input files are processed, but with filtering applied
103104
// - COPC: filtering inside readers.copc with "bounds" argument
104-
// - LAS/LAZ: filter either using CropFilter after reader -or- "where"
105+
// - LAS/LAZ: filter either using CropFilter after reader -or- "where"
105106

106107
// streaming algs:
107108
// - multi-las: if not overlapping: mode A
@@ -225,6 +226,16 @@ void runPipelineParallel(point_count_t totalPoints, bool isStreaming, std::vecto
225226

226227
std::string box_to_pdal_bounds(const BOX2D &box);
227228

229+
pdal::Bounds parseBounds(const std::string &boundsStr);
230+
231+
bool readerSupportsBounds(Stage &reader);
232+
233+
bool allReadersSupportBounds(const std::vector<Stage *> &readers);
234+
235+
BOX2D intersectionBox2D(const BOX2D &b1, const BOX2D &b2);
236+
237+
BOX2D intersectTileBoxWithFilterBox(const BOX2D &tileBox, const BOX2D &filterBox);
238+
228239

229240
inline bool ends_with(std::string const & value, std::string const & ending)
230241
{

‎external/pdal_wrench/vpc.cpp

Lines changed: 269 additions & 17 deletions
Large diffs are not rendered by default.

‎external/pdal_wrench/vpc.hpp

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,7 @@
1515
#include <vector>
1616
#include <string>
1717

18+
#include <pdal/Geometry.hpp>
1819
#include <pdal/pdal_types.hpp>
1920
#include <pdal/util/Bounds.hpp>
2021

@@ -26,6 +27,7 @@ void buildVpc(std::vector<std::string> args);
2627

2728
struct VirtualPointCloud
2829
{
30+
//! Schema for a single attribute
2931
struct SchemaItem
3032
{
3133
SchemaItem(const std::string &n, const std::string &t, int s): name(n), type(t), size(s) {}
@@ -35,14 +37,36 @@ struct VirtualPointCloud
3537
int size;
3638
};
3739

40+
//! Stats for a single attribute
41+
struct StatsItem
42+
{
43+
StatsItem(const std::string &n, uint32_t p, double a, point_count_t c, double max, double min, double st, double vr)
44+
: name(n), position(p), average(a), count(c), maximum(max), minimum(min), stddev(st), variance(vr) {}
45+
46+
std::string name;
47+
uint32_t position;
48+
double average;
49+
point_count_t count;
50+
double maximum;
51+
double minimum;
52+
double stddev;
53+
double variance;
54+
};
55+
3856
struct File
3957
{
4058
std::string filename;
4159
point_count_t count;
60+
std::string boundaryWkt; // not pdal::Geometry because of https://github.com/PDAL/PDAL/issues/4016
4261
BOX3D bbox;
4362
std::string crsWkt;
4463
std::string datetime; // RFC 3339 encoded date/time - e.g. 2023-01-01T12:00:00Z
4564
std::vector<SchemaItem> schema; // we're not using it, just for STAC export
65+
std::vector<StatsItem> stats;
66+
67+
// support for overview point clouds - currently we assume a file refers to at most a single overview file
68+
// (when building VPC with overviews, we create one overview file for all source data)
69+
std::string overviewFilename;
4670
};
4771

4872
std::vector<File> files;

0 commit comments

Comments
 (0)
Please sign in to comment.