Skip to content

Commit

Permalink
reproject points to clip rather than the other way around
Browse files Browse the repository at this point in the history
  • Loading branch information
@abellgithub
abellgithub committed Apr 30, 2021
1 parent 2879f6a commit 955f333
Show file tree
Hide file tree
Showing 4 changed files with 139 additions and 72 deletions.
110 changes: 73 additions & 37 deletions io/EptReader.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@
#include <pdal/pdal_features.hpp>
#include <pdal/util/ThreadPool.hpp>
#include <pdal/private/gdal/GDALUtils.hpp>
#include <pdal/private/SrsTransform.hpp>

#include "private/ept/Connector.hpp"
#include "private/ept/EptArtifact.hpp"
Expand Down Expand Up @@ -101,6 +102,12 @@ namespace Utils
}
}

struct PolySrs
{
Polygon poly;
SrsTransform xform;
};

struct EptReader::Args
{
public:
Expand Down Expand Up @@ -129,6 +136,8 @@ struct EptReader::Private
AddonList addons;
std::mutex mutex;
std::condition_variable contentsCv;
std::vector<PolySrs> polys;
SrsTransform boundsXform;
};

EptReader::EptReader() : m_args(new EptReader::Args), m_p(new EptReader::Private),
Expand All @@ -146,17 +155,13 @@ void EptReader::addArgs(ProgramArgs& args)
args.add("origin", "Origin of source file to fetch", m_args->m_origin);
args.add("threads", "Number of worker threads", m_args->m_threads);
args.add("resolution", "Resolution limit", m_args->m_resolution);
args.add("addons", "Mapping of addon dimensions to their output directory",
m_args->m_addons);
args.add("addons", "Mapping of addon dimensions to their output directory", m_args->m_addons);
args.add("polygon", "Bounding polygon(s) to crop requests",
m_args->m_polys).setErrorText("Invalid polygon specification. "
"Must be valid GeoJSON/WKT");
args.add("header", "Header fields to forward with HTTP requests",
m_args->m_headers);
args.add("query", "Query parameters to forward with HTTP requests",
m_args->m_query);
args.add("ogr", "OGR filter geometries",
m_args->m_ogr);
args.add("header", "Header fields to forward with HTTP requests", m_args->m_headers);
args.add("query", "Query parameters to forward with HTTP requests", m_args->m_query);
args.add("ogr", "OGR filter geometries", m_args->m_ogr);
}


Expand Down Expand Up @@ -218,34 +223,32 @@ void EptReader::initialize()
plist.insert(plist.end(), ogrPolys.begin(), ogrPolys.end());
}

// Transform query bounds to match point source SRS.
// Create transformations from our source data to the bounds SRS.
const SpatialReference& boundsSrs = m_args->m_bounds.spatialReference();
if (!m_p->info->srs().valid() && boundsSrs.valid())
throwError("Can't use bounds with SRS with data source that has "
"no SRS.");
m_queryBounds = m_args->m_bounds.to3d();
throwError("Can't use bounds with SRS with data source that has no SRS.");
if (boundsSrs.valid())
gdal::reprojectBounds(m_queryBounds,
boundsSrs.getWKT(), getSpatialReference().getWKT());
m_p->boundsXform = SrsTransform(m_p->info->srs(), boundsSrs);
m_queryBounds = m_args->m_bounds.to3d();

// Transform polygons and bounds to point source SRS.
std::vector <Polygon> exploded;
for (Polygon& poly : m_args->m_polys)
{
if (!poly.valid())
throwError("Geometrically invalid polygon in option 'polygon'.");

// Get the sub-polygons from a multi-polygon.
std::vector<Polygon> exploded = poly.polygons();
SrsTransform xform;
if (poly.srsValid())
xform.set(m_p->info->srs(), poly.getSpatialReference());
for (Polygon& p : exploded)
{
auto ok = poly.transform(getSpatialReference());
if (!ok)
throwError(ok.what());
PolySrs ps { std::move(p), xform };
m_p->polys.push_back(ps);
}
std::vector<Polygon> polys = poly.polygons();
exploded.insert(exploded.end(),
std::make_move_iterator(polys.begin()),
std::make_move_iterator(polys.end()));
}
m_args->m_polys = std::move(exploded);

try
{
Expand Down Expand Up @@ -549,20 +552,44 @@ void EptReader::overlaps()

void EptReader::overlaps(Hierarchy& target, const NL::json& hier, const Key& key)
{
// If this key doesn't overlap our query
// we can skip
if (!key.b.overlaps(m_queryBounds))
return;
auto reproject = [](BOX3D src, SrsTransform& xform) -> BOX3D
{
if (!xform.valid())
return src;

// Check the box of the key against our
// query polygon(s). If it doesn't overlap,
BOX3D b;
auto reprogrow = [&b, &xform](double x, double y, double z)
{
xform.transform(x, y, z);
b.grow(x, y, z);
};

reprogrow(src.minx, src.miny, src.minz);
reprogrow(src.maxx, src.miny, src.minz);
reprogrow(src.minx, src.maxy, src.minz);
reprogrow(src.maxx, src.maxy, src.minz);
reprogrow(src.minx, src.miny, src.maxz);
reprogrow(src.maxx, src.miny, src.maxz);
reprogrow(src.minx, src.maxy, src.maxz);
reprogrow(src.maxx, src.maxy, src.maxz);
return b;
};

{
std::lock_guard<std::mutex> lock(m_p->mutex);
// If the reprojected source bounds doesn't overlap our query bounds, we're done.
if (!reproject(key.b, m_p->boundsXform).overlaps(m_queryBounds))
return;
}

// Check the box of the key against our query polygon(s). If it doesn't overlap,
// we can skip
auto polysOverlap = [this, &key]()
auto polysOverlap = [this, &reproject, &key]()
{
if (m_args->m_polys.empty())
if (m_p->polys.empty())
return true;
for (auto& p: m_args->m_polys)
if (!p.disjoint(key.b))
for (auto& ps : m_p->polys)
if (!ps.poly.disjoint(reproject(key.b, ps.xform)))
return true;
return false;
};
Expand Down Expand Up @@ -652,22 +679,31 @@ bool EptReader::processPoint(PointRef& dst, const TileContents& tile)
if (m_queryOriginId != -1 && originId != m_queryOriginId)
return false;

auto passesPolyFilter = [this](double x, double y)
auto passesBoundsFilter = [this](double x, double y, double z)
{
if (m_args->m_polys.empty())
m_p->boundsXform.transform(x, y, z);
return m_queryBounds.contains(x, y, z);
};

auto passesPolyFilter = [this](double x, double y, double z)
{
if (m_p->polys.empty())
return true;

for (Polygon& poly : m_args->m_polys)
if (poly.contains(x, y))
for (PolySrs& ps : m_p->polys)
{
ps.xform.transform(x, y, z);
if (ps.poly.contains(x, y))
return true;
}
return false;
};

double x = p.getFieldAs<double>(Id::X);
double y = p.getFieldAs<double>(Id::Y);
double z = p.getFieldAs<double>(Id::Z);

if (!m_queryBounds.contains(x, y, z) || !passesPolyFilter(x, y))
if (!passesBoundsFilter(x, y, z) || !passesPolyFilter(x, y, z))
return false;

for (auto& el : m_p->info->dims())
Expand Down
63 changes: 46 additions & 17 deletions pdal/private/SrsTransform.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,25 +36,36 @@
namespace pdal
{

SrsTransform::SrsTransform(const SpatialReference& src, const SpatialReference& dst) :
SrsTransform(OGRSpatialReference(src.getWKT().data()),
OGRSpatialReference(dst.getWKT().data()))
SrsTransform::SrsTransform()
{}

SrsTransform::SrsTransform(const SrsTransform& src) : m_transform(src.m_transform->Clone())
{}


SrsTransform::SrsTransform(SrsTransform&& src) : m_transform(std::move(src.m_transform))
{}


SrsTransform::SrsTransform(OGRSpatialReference srcRef, OGRSpatialReference dstRef)
SrsTransform::~SrsTransform()
{}


SrsTransform& SrsTransform::operator=(SrsTransform&& src)
{
// Starting with version 3 of GDAL, the axes (X, Y, Z or lon, lat, h or whatever)
// are mapped according to the WKT definition. In particular, this means
// that for EPSG:4326 the mapping is X -> lat, Y -> lon, rather than the
// more conventional X -> lon, Y -> lat. Setting this flag reverses things
// such that the traditional ordering is maintained. There are other
// SRSes where this comes up. See "axis order issues" in the GDAL WKT2
// discussion for more info.
//
srcRef.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
dstRef.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
m_transform.reset(OGRCreateCoordinateTransformation(&srcRef, &dstRef));
m_transform = std::move(src.m_transform);
return *this;
}

SrsTransform::SrsTransform(const SpatialReference& src, const SpatialReference& dst)
{
set(src, dst);
}


SrsTransform::SrsTransform(const OGRSpatialReference& srcRef, const OGRSpatialReference& dstRef)
{
set(srcRef, dstRef);
}


Expand All @@ -81,8 +92,26 @@ SrsTransform::SrsTransform(const SpatialReference& src,
m_transform.reset(OGRCreateCoordinateTransformation(&srcRef, &dstRef));
}

SrsTransform::~SrsTransform()
{}
void SrsTransform::set(const SpatialReference& src, const SpatialReference& dst)
{
set(OGRSpatialReference(src.getWKT().data()), OGRSpatialReference(dst.getWKT().data()));
}


void SrsTransform::set(OGRSpatialReference src, OGRSpatialReference dst)
{
// Starting with version 3 of GDAL, the axes (X, Y, Z or lon, lat, h or whatever)
// are mapped according to the WKT definition. In particular, this means
// that for EPSG:4326 the mapping is X -> lat, Y -> lon, rather than the
// more conventional X -> lon, Y -> lat. Setting this flag reverses things
// such that the traditional ordering is maintained. There are other
// SRSes where this comes up. See "axis order issues" in the GDAL WKT2
// discussion for more info.
//
src.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
dst.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
m_transform.reset(OGRCreateCoordinateTransformation(&src, &dst));
}


OGRCoordinateTransformation *SrsTransform::get() const
Expand Down
20 changes: 15 additions & 5 deletions pdal/private/SrsTransform.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,14 +44,19 @@ class PDAL_DLL SrsTransform
public:
/// Object that performs transformation from a \src spatial reference
/// to a \dest spatial reference.
SrsTransform(OGRSpatialReference src, OGRSpatialReference dst);
SrsTransform();
SrsTransform(const SrsTransform& src);
SrsTransform(SrsTransform&& src);
SrsTransform& operator=(SrsTransform&& src);
SrsTransform(const OGRSpatialReference& src, const OGRSpatialReference& dst);
SrsTransform(const SpatialReference& src, const SpatialReference& dst);
SrsTransform(const SpatialReference& src,
std::vector<int> srcOrder,
const SpatialReference& dst,
std::vector<int> dstOrder);
SrsTransform(const SpatialReference& src, std::vector<int> srcOrder,
const SpatialReference& dst, std::vector<int> dstOrder);
~SrsTransform();

void set(OGRSpatialReference src, OGRSpatialReference dst);
void set(const SpatialReference& src, const SpatialReference& dst);

/// Get the underlying transformation.
/// \return Pointer to the underlying coordinate transform.
OGRCoordinateTransformation *get() const;
Expand All @@ -71,6 +76,11 @@ class PDAL_DLL SrsTransform
bool transform(std::vector<double>& x, std::vector<double>& y,
std::vector<double>& z) const;

/// Determine if this represents a valid transform.
/// \return Whether the transform is valid or not.
bool valid() const
{ return (bool)m_transform.get(); }

private:
std::unique_ptr<OGRCoordinateTransformation> m_transform;
};
Expand Down
18 changes: 5 additions & 13 deletions test/unit/io/EptReaderTest.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -292,8 +292,6 @@ TEST(EptReaderTest, bounds2dXform)
{
BOX2D b(515380, 4918360, 515390, 4918370);
SrsBounds eptBounds(b);
gdal::reprojectBounds(b, "EPSG:26912", "EPSG:4326");
SrsBounds boxBounds(b, "EPSG:4326");

PointViewPtr v1;
PointViewPtr v2;
Expand All @@ -308,6 +306,9 @@ TEST(EptReaderTest, bounds2dXform)
auto vset = reader.execute(eptTable);
v1 = *vset.begin();
}

gdal::reprojectBounds(b, "EPSG:26912", "EPSG:4326");
SrsBounds boxBounds(b, "EPSG:4326");
{
EptReader reader;
Options options;
Expand Down Expand Up @@ -622,8 +623,7 @@ TEST(EptReaderTest, boundedCrop)
{
Options options;
options.add("filename", eptAutzenPath);
std::string overrides(wkt + "/ EPSG:3644");
Option polygon("polygon", overrides);
Option polygon("polygon", wkt + "/ EPSG:3644");
options.add(polygon);
reader.setOptions(options);
}
Expand All @@ -647,9 +647,8 @@ TEST(EptReaderTest, boundedCrop)
}
CropFilter crop;
{
std::string overrides(wkt + "/ EPSG:3644");
Options options;
Option polygon("polygon", overrides);
Option polygon("polygon", wkt + "/ EPSG:3644");
options.add(polygon);
crop.setOptions(options);
crop.setInput(source);
Expand Down Expand Up @@ -734,7 +733,6 @@ TEST(EptReaderTest, boundedCropReprojection)

TEST(EptReaderTest, ogrCrop)
{

EptReader reader;
{
Options options;
Expand All @@ -754,9 +752,7 @@ TEST(EptReaderTest, ogrCrop)

uint64_t eptNp(0);
for (const PointViewPtr& view : reader.execute(eptTable))
{
eptNp += view->size();
}

// Now we'll check the result against a crop filter of the source file with
// the same bounds.
Expand All @@ -770,15 +766,11 @@ TEST(EptReaderTest, ogrCrop)
source.prepare(sourceTable);
uint64_t sourceNp(0);
for (const PointViewPtr& view : source.execute(sourceTable))
{
sourceNp += view->size();
}

EXPECT_EQ(eptNp, sourceNp);
EXPECT_EQ(eptNp, 86u);
EXPECT_EQ(sourceNp, 86u);
}



} // namespace pdal

0 comments on commit 955f333

Please sign in to comment.