initial implementation of OGR filtering with an 'ogr' json block for …

…readers.ept

io/EptReader.cpp

@@ -106,6 +106,8 @@ class EptBounds : public SrsBounds
     EptBounds() : SrsBounds(BOX3D(LOWEST, LOWEST, LOWEST,
         HIGHEST, HIGHEST, HIGHEST))
     {}
+    EptBounds(const SrsBounds& b) : SrsBounds( b )
+    {}
 };
 
 namespace Utils
@@ -140,6 +142,7 @@ struct EptReader::Args
 
     NL::json m_query;
     NL::json m_headers;
+    NL::json m_ogr;
 };
 
 EptReader::EptReader() : m_args(new EptReader::Args), m_currentIndex(0)
@@ -165,6 +168,8 @@ void EptReader::addArgs(ProgramArgs& args)
         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);
 }
 
 
@@ -243,6 +248,48 @@ void EptReader::initialize()
 
     setSpatialReference(m_info->srs());
 
+    if (!m_args->m_ogr.is_null())
+    {
+        // convert whatever we were given for an OGR source to
+        // polygons and override m_args->m_polys
+        m_args->m_polys.clear();
+        std::vector<OGRGeometry*> ogr_geoms = pdal::gdal::fetchOGRGeometries(m_args->m_ogr);
+        BOX3D layerBounds;
+
+        SpatialReference layerSrs;
+        for (auto g: ogr_geoms)
+        {
+
+            OGRSpatialReferenceH srs = (OGRSpatialReferenceH) g->getSpatialReference();
+            Polygon p = pdal::Polygon(g, srs);
+            m_args->m_polys.emplace_back(p);
+            if (layerSrs.empty())
+            {
+                SpatialReference gSrs = SpatialReference(srs);
+                if (!gSrs.empty())
+                    layerSrs = gSrs;
+            }
+
+            BOX3D box = p.bounds();
+            layerBounds.grow(box.minx, box.miny, EptBounds::LOWEST);
+
+            layerBounds.grow(box.maxx, box.maxy, EptBounds::HIGHEST);
+            std::cout << "layerBounds: " << layerBounds << std::endl;
+            std::cout << "box: " << box << std::endl;
+        }
+
+        // Apply our overrides from the fetched layer by taking
+        // the intersection of the layer's bounds with whatever the
+        // user provided before
+        BOX3D filtered = m_args->m_bounds.to3d();
+        std::cout << "filtered: " << filtered << std::endl;
+        filtered.clip(layerBounds);
+        std::cout << "filtered: " << filtered << std::endl;
+        SrsBounds newBounds = SrsBounds(filtered, layerSrs);
+    }
+
+
+
     // Transform query bounds to match point source SRS.
     const SpatialReference& boundsSrs = m_args->m_bounds.spatialReference();
     if (!m_info->srs().valid() && boundsSrs.valid())
@@ -253,6 +300,7 @@ void EptReader::initialize()
         gdal::reprojectBounds(m_queryBounds,
             boundsSrs.getWKT(), getSpatialReference().getWKT());
 
+
     // Transform polygons and bounds to point source SRS.
     std::vector <Polygon> exploded;
     for (Polygon& poly : m_args->m_polys)

io/EptReader.hpp

@@ -158,6 +158,8 @@ class PDAL_DLL EptReader : public Reader, public Streamable
     point_count_t m_currentIndex = -1;
     std::vector<char> m_temp_buffer;
     std::vector<std::unique_ptr<GridPnp>> m_queryGrids;
+
+
 };
 
 } // namespace pdal

pdal/GDALUtils.cpp

@@ -44,6 +44,9 @@
 
 #include <ogr_spatialref.h>
 #include <ogr_p.h>
+#include <ogr_api.h>
+#include <ogrsf_frmts.h>
+#include <nlohmann/json.hpp>
 
 #pragma warning(disable: 4127)  // conditional expression is constant
 
@@ -774,6 +777,86 @@ OGRGeometry *createFromGeoJson(const std::string& s, std::string& srs)
     return newGeom;
 }
 
+
+std::vector<OGRGeometry*> fetchOGRGeometries(const NL::json ogr)
+{
+//    const NL::json el = json.value("ogr", NL::json::null_value);
+    const NL::json& layer = ogr.at("layer");
+    const NL::json& datasource = ogr.at("datasource");
+
+    const NL::json options = ogr.at("options");//, nullptr);
+
+    registerDrivers();
+    std::string dsString = datasource.get<std::string>();
+    unsigned int openFlags = GDAL_OF_READONLY | GDAL_OF_VECTOR | GDAL_OF_VERBOSE_ERROR;
+    GDALDataset* ds = (GDALDataset*) GDALOpenEx( dsString.c_str(),
+                                                 openFlags,
+                                                 NULL, NULL, NULL );
+    if (!ds)
+        throw pdal_error("Unable to read datasource in fetchOGRGeometries " + datasource.dump() );
+
+
+    OGRLayer* poLayer(nullptr);
+
+    std::string lyrString = layer.get<std::string>();
+    poLayer = ds->GetLayerByName( lyrString.c_str() );
+
+    if (!poLayer)
+        throw pdal_error("Unable to read layer in fetchOGRGeometries " + layer.dump() );
+
+    std::vector<OGRGeometry*> output;
+
+    OGRFeature *poFeature (nullptr);
+    OGRGeometry* filterGeometry(nullptr);
+
+    poLayer->ResetReading();
+    if (ogr.count("sql"))
+    {
+        std::string dialect("OGRSQL");
+        if (options.count("dialect"))
+            dialect = options.at("dialect").get<std::string>();
+
+        if (options.count("geometry"))
+        {
+            std::string wkt_or_json = options.at("geometry").get<std::string>();
+            std::cout << "filter geometry!: " << wkt_or_json << std::endl;
+            bool isJson = (wkt_or_json.find("{") != wkt_or_json.npos) ||
+                          (wkt_or_json.find("}") != wkt_or_json.npos);
+
+            std::string srs;
+            if (isJson)
+            {
+                filterGeometry = createFromGeoJson(wkt_or_json, srs);
+                if (!filterGeometry)
+                    throw pdal_error("Unable to create filter geometry from input GeoJSON");
+            }
+            else
+            {
+                filterGeometry = createFromWkt(wkt_or_json, srs);
+                if (!filterGeometry)
+                    throw pdal_error("Unable to create filter geometry from input WKT");
+            }
+
+        }
+
+        std::string query = ogr.at("sql").get<std::string>();
+        std::cout << "query '" << query << "'" << std::endl;
+        poLayer = ds->ExecuteSQL( query.c_str(),
+                                  filterGeometry,
+                                  dialect.c_str());
+        if (!poLayer)
+            throw pdal_error("unable to execute sql query!");
+    }
+    while( (poFeature = poLayer->GetNextFeature()) != NULL )
+    {
+        OGRGeometry* poGeometry = poFeature->GetGeometryRef();
+        OGRGeometry* clone = poGeometry->clone();
+        output.emplace_back(clone);
+        OGRFeature::DestroyFeature( poFeature );
+    }
+    return output;
+}
+
 } // namespace gdal
 
 namespace oldgdalsupport

pdal/GDALUtils.hpp

@@ -45,6 +45,7 @@
 #include <pdal/Log.hpp>
 #include <pdal/SpatialReference.hpp>
 #include <pdal/util/Bounds.hpp>
+#include <pdal/JsonFwd.hpp>
 
 #include <cpl_conv.h>
 #include <gdal_priv.h>
@@ -866,6 +867,8 @@ OGRGeometry *createFromGeoJson(const char *s);
 OGRGeometry *createFromWkt(const std::string& s, std::string& srs);
 OGRGeometry *createFromGeoJson(const std::string& s, std::string& srs);
 
+std::vector<OGRGeometry*> fetchOGRGeometries(const NL::json json);
+
 inline OGRGeometry *fromHandle(OGRGeometryH geom)
 { return reinterpret_cast<OGRGeometry *>(geom); }
 

pdal/Geometry.cpp

@@ -74,6 +74,19 @@ Geometry::Geometry(OGRGeometryH g, const SpatialReference& srs) :
 }
 
 
+Geometry::Geometry(OGRGeometryH gh, OGRSpatialReferenceH srs) :
+    m_geom((reinterpret_cast<OGRGeometry *>(gh))->clone())
+{
+    OGRGeometry* g = reinterpret_cast<OGRGeometry *>(gh);
+    if (g)
+    {
+        OGRSpatialReferenceH ref = (OGRSpatialReferenceH) g->getSpatialReference();
+        OGR_G_AssignSpatialReference((OGRGeometry*)m_geom.get(), ref);
+        SpatialReference s(ref);
+        setSpatialReference(s);
+    }
+}
+
 Geometry::~Geometry()
 {}
 

pdal/Geometry.hpp

@@ -41,6 +41,7 @@
 
 class OGRGeometry;
 using OGRGeometryH = void *;
+using OGRSpatialReferenceH = void *;
 
 namespace pdal
 {
@@ -56,6 +57,7 @@ class PDAL_DLL Geometry
     Geometry(const std::string& wkt_or_json,
            SpatialReference ref = SpatialReference());
     Geometry(OGRGeometryH g, const SpatialReference& srs);
+    Geometry(OGRGeometryH g, OGRSpatialReferenceH srs);
 
 public:
     Geometry& operator=(const Geometry&);

pdal/Polygon.cpp

@@ -59,6 +59,28 @@ Polygon::Polygon(OGRGeometryH g, const SpatialReference& srs) : Geometry(g, srs)
     }
 }
 
+Polygon::Polygon(OGRGeometryH g, OGRSpatialReferenceH srs) : Geometry(g, srs)
+{
+    // If the handle was null, we need to create an empty polygon.
+    if (!m_geom)
+    {
+        m_geom.reset(new OGRPolygon());
+        return;
+    }
+
+    OGRwkbGeometryType t = m_geom->getGeometryType();
+
+    if (!(t == wkbPolygon ||
+        t == wkbMultiPolygon ||
+        t == wkbPolygon25D ||
+        t == wkbMultiPolygon25D))
+    {
+        throw pdal::pdal_error("pdal::Polygon() cannot construct geometry "
+            "because OGR geometry is not Polygon or MultiPolygon.");
+    }
+
+}
+
 Polygon::Polygon(const BOX2D& box)
 {
     OGRPolygon *poly = new OGRPolygon();

pdal/Polygon.hpp

@@ -59,6 +59,7 @@ class PDAL_DLL Polygon : public Geometry
     Polygon(const BOX2D&);
     Polygon(const BOX3D&);
     Polygon(OGRGeometryH g, const SpatialReference& srs);
+    Polygon(OGRGeometryH g, OGRSpatialReferenceH srs);
 
     OGRGeometryH getOGRHandle();
 

pdal/SpatialReference.cpp

@@ -81,6 +81,19 @@ SpatialReference::SpatialReference(const std::string& s)
     set(s);
 }
 
+SpatialReference::SpatialReference(OGRSpatialReferenceH ref)
+{
+
+    OGRSpatialReference* s = reinterpret_cast<OGRSpatialReference *>(ref);
+    if (s)
+    {
+        char *poWKT = 0;
+        s->exportToWkt(&poWKT);
+        m_wkt = poWKT;
+        CPLFree(poWKT);
+    }
+
+}
 
 bool SpatialReference::empty() const
 {

pdal/SpatialReference.hpp

@@ -36,6 +36,8 @@
 
 #include <pdal/pdal_internal.hpp>
 
+using OGRSpatialReferenceH = void *;
+
 namespace pdal
 {
 
@@ -70,6 +72,13 @@ class PDAL_DLL SpatialReference
     */
     SpatialReference(const std::string& wkt);
 
+    /**
+      Construct a spatial reference from an OGRSpatialReferenceH
+
+      \param srs  OGRSpatialReferenceH to clone via WKT translation.
+    */
+    SpatialReference(OGRSpatialReferenceH ref);
+
     /**
       Determine if this spatial reference is the same as another.