support fetching bounds information for readers.gdal QuickInfo (#2177)

* support fetching bounds information for readers.gdal QuickInfo

* rearrange into Raster::bounds method

* refactor to return a BOX2D or a BOX3D if given a bandID

* add 'header' option to readers.gdal to allow mapping of band->dimension

* test 'header' option

doc/stages/readers.gdal.rst

@@ -85,3 +85,8 @@ filename
 
 count
     Maximum number of points to read [Optional]
+
+header
+    A comma-separated list of :ref:`dimensions` IDs to map
+    bands to. The length of the list must match the number
+    of bands in the raster.

io/GDALReader.cpp

@@ -47,7 +47,7 @@ static StaticPluginInfo const s_info
     "readers.gdal",
     "Read GDAL rasters as point clouds.",
     "http://pdal.io/stages/reader.gdal.html",
-    { "tif", "tiff", "jpeg", "jpg" }
+    { "tif", "tiff", "jpeg", "jpg", "png" }
 
 };
 
@@ -92,15 +92,24 @@ QuickInfo GDALReader::inspect()
     QuickInfo qi;
     std::unique_ptr<PointLayout> layout(new PointLayout());
 
-    addDimensions(layout.get());
     initialize();
+    addDimensions(layout.get());
 
     m_raster = std::unique_ptr<gdal::Raster>(new gdal::Raster(m_filename));
     if (m_raster->open() == gdal::GDALError::CantOpen)
         throwError("Couldn't open raster file '" + m_filename + "'.");
 
     qi.m_pointCount = m_raster->width() * m_raster->height();
-    // qi.m_bounds = ???;
+
+    auto p = std::find(m_bandIds.begin(), m_bandIds.end(), Dimension::Id::Z);
+
+    int nBand(1);
+    if (p != m_bandIds.end())
+    {
+        nBand = (int) std::distance(m_bandIds.begin(), p);
+    }
+
+    qi.m_bounds = m_raster->bounds(nBand);
     qi.m_srs = m_raster->getSpatialRef();
     qi.m_valid = true;
 
@@ -112,15 +121,41 @@ void GDALReader::addDimensions(PointLayoutPtr layout)
 {
     layout->registerDim(pdal::Dimension::Id::X);
     layout->registerDim(pdal::Dimension::Id::Y);
-    for (int i = 0; i < m_raster->bandCount(); ++i)
+
+    if (m_header.size())
+    {
+        std::vector<std::string> dimNames = Utils::split(m_header, ',');
+
+        if (dimNames.size() != (size_t)m_raster->bandCount())
+            throwError("Dimension names are not the same count as raster bands!");
+
+        for (int i=0; i < m_raster->bandCount(); ++i)
+        {
+
+            std::cerr << "Adding dimension '" << dimNames[i] << "'" << std::endl;
+            Dimension::Id id = layout->registerOrAssignDim(dimNames[i], m_bandTypes[i]);
+            m_bandIds.push_back(id);
+        }
+
+    }
+    else
     {
-        std::ostringstream oss;
-        oss << "band-" << (i + 1);
-        Dimension::Id id = layout->registerOrAssignDim(oss.str(), m_bandTypes[i]);
-        m_bandIds.push_back(id);
+        for (int i = 0; i < m_raster->bandCount(); ++i)
+        {
+            std::ostringstream oss;
+            oss << "band-" << (i + 1);
+            Dimension::Id id = layout->registerOrAssignDim(oss.str(), m_bandTypes[i]);
+            m_bandIds.push_back(id);
+        }
     }
 }
 
+void GDALReader::addArgs(ProgramArgs& args)
+{
+    args.add("header", "A comma-separated list of dimension IDs to map raster bands to dimension id", m_header);
+}
+
+
 
 void GDALReader::ready(PointTableRef table)
 {

io/GDALReader.hpp

@@ -64,10 +64,12 @@ class PDAL_DLL GDALReader : public Reader , public Streamable
         { m_raster->close(); }
     virtual bool processOne(PointRef& point);
     virtual QuickInfo inspect();
+    virtual void addArgs(ProgramArgs& args);
 
     std::unique_ptr<gdal::Raster> m_raster;
     std::vector<Dimension::Type> m_bandTypes;
     std::vector<Dimension::Id> m_bandIds;
+    std::string m_header;
     point_count_t m_index;
     int m_row;
     int m_col;

pdal/GDALUtils.hpp

@@ -511,8 +511,19 @@ friend class Raster;
         if (m_band->WriteBlock(x, y, m_buf.data()) != CPLE_None)
             throw CantWriteBlock();
     }
+
+    void statistics(double* minimum, double* maximum,
+                    double* mean, double* stddev,
+                    int bApprox, int bForce) const
+    {
+        m_band->GetStatistics(bApprox, bForce, minimum, maximum, mean, stddev);
+    }
+
+
+
 };
 
+
 class PDAL_DLL Raster
 {
 public:
@@ -733,6 +744,51 @@ class PDAL_DLL Raster
 
     std::string const& filename() { return m_filename; }
 
+    void statistics(int nBand,
+                    double* minimum,
+                    double* maximum,
+                    double* mean,
+                    double* stddev,
+                    int bApprox=TRUE,
+                    int bForce=TRUE) const
+    {
+        Band<double>(m_ds, nBand).statistics(minimum, maximum, mean, stddev, bApprox, bForce);
+    }
+
+    BOX2D bounds() const
+    {
+
+        std::array<double, 2> coords;
+
+        pixelToCoord(height(),
+                     width(),
+                     coords);
+        double maxx = coords[0];
+        double maxy = coords[1];
+
+        pixelToCoord(0,
+                     0,
+                     coords);
+        double minx = coords[0];
+        double miny = coords[1];
+        BOX2D output(minx, miny, maxx, maxy);
+        return output;
+    }
+
+    BOX3D bounds(int nBand) const
+    {
+
+        BOX2D box2 = bounds();
+
+        double minimum; double maximum;
+        double mean; double stddev;
+        statistics(nBand, &minimum, &maximum, &mean, &stddev);
+
+        BOX3D output(box2.minx, box2.miny, minimum,
+                     box2.maxx, box2.maxy, maximum);
+        return output;
+    }
+
 private:
     std::string m_filename;
 

test/unit/io/GDALReaderTest.cpp

@@ -60,6 +60,7 @@ TEST(GDALReaderTest, simple)
 {
     Options ro;
     ro.add("filename", Support::datapath("png/autzen-height.png"));
+    ro.add("header", "Intensity,Userdata,Z");
 
     GDALReader gr;
     gr.setOptions(ro);
@@ -69,9 +70,9 @@ TEST(GDALReaderTest, simple)
     PointViewSet s = gr.execute(t);
     PointViewPtr v = *s.begin();
     PointLayoutPtr l = t.layout();
-    Dimension::Id id1 = l->findDim("band-1");
-    Dimension::Id id2 = l->findDim("band-2");
-    Dimension::Id id3 = l->findDim("band-3");
+    Dimension::Id id1 = l->findDim("Intensity");
+    Dimension::Id id2 = l->findDim("Userdata");
+    Dimension::Id id3 = l->findDim("Z");
     EXPECT_EQ(v->size(), (size_t)(735 * 973));
 
     auto verify = [v, id1, id2, id3]