Merge remote-tracking branch 'origin/master' into filters-colorinterp…

…olate

doc/stages/writers.gdal.rst

@@ -0,0 +1,127 @@
+.. _writers.gdal:
+
+writers.gdal
+================================================================================
+
+The `GDAL`_ writer creates a raster from a point cloud using an interpolation
+algorithm.  Output is produced using GDAL and can therefore use any `driver
+that supports creation of rasters`_.
+
+.. _`GDAL`: http://gdal.org
+.. _`driver that supports creation of rasters`: http://www.gdal.org/formats_list.html
+
+The technique used to create the raster is a simple interpolation where
+each point that falls within a given radius of a raster cell center
+potentially contributes to the raster's value.
+
+.. note::
+
+    If a circle with the provided radius doesn't encompass the entire cell,
+    it is possible that some points will not be considered at all, including
+    those that may be within the bounds of the raster cell.
+
+The GDAL writer creates rasters using the data in the 'Z' dimension of the 
+points.  In a :ref:`pipeline_command` one can precede the GDAL writer with
+the ferry filter (:ref:`filters.ferry`) in order to copy data from some
+other dimension to the 'Z' dimension if desired.
+
+The writer will creates up to six rasters based on different statistics in
+the output dataset.  The order of the layers in the dataset is as follows:
+
+min
+    Give the cell the minimum value of all points within the given radius.
+
+max
+    Give the cell the maximum value of all points within the given radius.
+
+mean
+    Give the cell the mean value of all points within the given radius.
+
+idw
+    Cells are assigned a value based on `Shepard's inverse distance weighting`_
+    algorithm, considering all points within the given radius.
+
+count
+    Give the cell the number of points that lie within the given radius.
+
+stdev
+    Give the cell the population standard deviation of the points that lie
+    within the given radius.
+
+.. _`Shepard's inverse distance weighting`: https://en.wikipedia.org/wiki/Inverse_distance_weighting
+
+If no points fall within the circle about a raster cell, a secondary
+algorithm can be used to attempt to provide a value after the standard
+interpolation is complete.  If the window_size_ option is set to a non-zero
+value, a square of rasters surrounding an empty cell, and the value of each
+non-empty surrounding is averaged using inverse distance weighting to provide
+a value for the subject cell.  The value provided for window_size is the
+maximum horizontal or vertical distance that a donor cell may be in order to
+contribute to the subject cell (A window_size of 1 essentially creates a 3x3
+array around the subject cell.  A window_size of 2 creates a 5x5 array, and
+so on.)
+
+Cells that have no value after interpolation are given the empty value of -9999.
+
+Basic Example
+--------------------------------------------------------------------------------
+
+This  pipeline reads the file autzen_trim.las and creates a Geotiff dataset
+called outputfile.tif.  Since output_type isn't specified, it creates six
+raster bands ("min", "max", "mean", "idx", "count" and "stdev") in the output
+dataset.  The raster cells are 10x10 and the radius used to locate points
+whose values contribute to the cell value is 14.14.
+
+.. code-block:: json
+
+    {
+      "pipeline":[
+        "pdal/test/data/las/autzen_trim.las",
+        {
+          "edge_length": 10,
+          "radius": 14.14,
+          "filename":"outputfile.tif"
+        }
+      ]
+    }
+
+
+Options
+--------------------------------------------------------------------------------
+
+filename
+    Name of output file. [Required]
+
+edge_length
+    Length of raster cell edges.  [Required]
+
+radius
+    Radius about cell center bounding points to use to calculate a cell value.
+    [Required]
+
+gdaldriver
+    Name of the GDAL driver to use to write the output. [Default: "GTiff"]
+
+gdalopts
+    A list of key/value options to pass directly to the GDAL driver.  The
+    format is name=value,name=value,...  The option may be specified
+    any number of times.
+
+    .. note::
+        The INTERLEAVE GDAL driver option is not supported.  writers.gdal
+        always uses BAND interleaving.
+
+.. _output_type:
+
+output_type
+    A comma separated list of statistics for which to produce raster layers.
+    The supported values are "min", "max", "mean", "idw", "count", "stdev"
+    and "all".  The option may be specified more than once. [Default: "all"]
+
+.. _window_size:
+
+window_size
+    The maximum distance from a donor cell to a target cell when applying
+    the fallback interpolation method.  See the stage description for more
+    information. [Default: 0]
+

include/pdal/GDALUtils.hpp

@@ -36,6 +36,7 @@
 
 #include <pdal/pdal_internal.hpp>
 #include <pdal/Dimension.hpp>
+#include <pdal/SpatialReference.hpp>
 #include <pdal/util/Bounds.hpp>
 
 #include <pdal/Log.hpp>
@@ -46,7 +47,7 @@
 #include <vector>
 
 #include <cpl_port.h>
-#include <gdal.h>
+#include <gdal_priv.h>
 #include <cpl_vsi.h>
 #include <cpl_conv.h>
 #include <ogr_api.h>
@@ -137,7 +138,7 @@ class Geometry
         }
         OGRErr err = OGR_G_CreateFromWkt(&p_wkt, ref, &geom);
         if (err != OGRERR_NONE)
-            throw pdal::pdal_error("unable to construct OGR geometry from wkt!");
+            throw pdal::pdal_error("Unable to construct OGR geometry from wkt");
         newRef(geom);
     }
 
@@ -236,7 +237,6 @@ class PDAL_DLL ErrorHandler
 };
 
 
-
 enum class GDALError
 {
     None,
@@ -246,51 +246,162 @@ enum class GDALError
     InvalidBand,
     NoTransform,
     NotInvertible,
-    CantReadBlock
+    CantReadBlock,
+    InvalidDriver,
+    DriverNotFound,
+    CantCreate,
+    InvalidOption,
+    CantWriteBlock
 };
 
 class PDAL_DLL Raster
 {
 
 public:
-    Raster(const std::string& filename);
+    /**
+      Constructor.
+
+      \param filename  Filename of raster file.
+      \param drivername  Optional name of driver to use to open raster file.
+    */
+    Raster(const std::string& filename, const std::string& drivername = "");
+
+    /**
+      Constructor.
+
+      \param filename  Filename of raster file.
+      \param drivername  Optional name of driver to use to open raster file.
+      \param srs  SpatialReference of the raster.
+      \param pixelToPos  Transformation matrix to convert raster positions to
+        geolocations.
+    */
+    Raster(const std::string& filename, const std::string& drivername,
+        const SpatialReference& srs, const std::array<double, 6> pixelToPos);
+
+
+    /**
+      Destructor.  Closes an open raster.
+    */
     ~Raster();
+
+    /**
+      Open raster file for reading.
+    */
     GDALError open();
+
+    /**
+      Open a raster for writing.
+
+      \param width  Width of the raster in cells (X direction)
+      \param height  Height of the raster in cells (Y direction)
+      \param numBands  Number of bands in the raster.
+      \param type  Datatype (int, float, etc.) of the raster data.
+      \param noData  Value that indiciates no data in a raster cell.
+      \param options  GDAL driver options.
+    */
+    GDALError open(int width, int height, int numBands, Dimension::Type type,
+        double noData, StringList options = StringList());
+
+    /**
+      Close the raster and deallocate the underlying dataset.
+    */
     void close();
 
-    GDALError read(double x, double y, std::vector<double>& data);
-    std::vector<pdal::Dimension::Type> getPDALDimensionTypes() const
-       { return m_types; }
     /**
-      Read a raster band (layer) into a vector.
+      Read an entire raster band (layer) into a vector.
 
       \param band  Vector into which data will be read.  The vector will
         be resized appropriately to hold the data.
       \param nBand  Band number to read.  Band numbers start at 1.
+      \return Error code or GDALError::None.
     */
     GDALError readBand(std::vector<uint8_t>& band, int nBand);
 
+    /**
+      Write an entire raster band (layer) into raster to be written with GDAL.
+
+      \param data  Linearized raster data to be written.
+      \param nBand  Band number to write.
+      \param name  Name of the raster band.
+    */
+    GDALError writeBand(const uint8_t *data, int nBand,
+        const std::string& name = "");
+
+    /**
+      Read the data for each band at x/y into a vector of doubles.  x and y
+      are transformed to the basis of the raster before the data is fetched.
+
+      \param x  X position to read
+      \param y  Y position to read
+      \param data  Vector in which to store data.
+    */
+    GDALError read(double x, double y, std::vector<double>& data);
+
+    /**
+      Get a vector of dimensions that map to the bands of a raster. 
+    */
+    std::vector<pdal::Dimension::Type> getPDALDimensionTypes() const
+       { return m_types; }
+
+    /**
+      Convert an X/Y raster position into geo-located position using the
+      raster's transformation matrix.
+
+      \param column  raster column whose position should be calculated
+      \param row  raster row whose position should be calculated
+      \param[out]  Array containing the geo-located position of the pixel.
+    */
     void pixelToCoord(int column, int row, std::array<double, 2>& output) const;
+
+    /**
+      Get the spatial reference associated with the raster.
+
+      \return  The associated spatial reference.
+    */
     SpatialReference getSpatialRef() const;
+
+    /**
+      Get the most recent error message.
+    */
     std::string errorMsg() const
         { return m_errorMsg; }
 
-    std::string m_filename;
+    /**
+      Get the number of bands in the raster.
 
-    std::array<double, 6> m_forward_transform;
-    std::array<double, 6> m_inverse_transform;
+      \return  The number of bands in the raster.
+    */
+    int bandCount() const
+        { return m_numBands; }
 
-    int m_raster_x_size;
-    int m_raster_y_size;
+    /**
+      Get the width of the raster (X direction)
+    */
+    int width() const
+        { return m_width; }
 
-    int m_band_count;
-    mutable std::vector<pdal::Dimension::Type> m_types;
-    std::vector<std::array<double, 2>> m_block_sizes;
+    /**
+      Get the height of the raster (Y direction)
+    */
+    int height() const
+        { return m_height; }
+
+private:
+    std::string m_filename;
+
+    int m_width;
+    int m_height;
+    int m_numBands;
+    std::string m_drivername;
+    std::array<double, 6> m_forwardTransform;
+    std::array<double, 6> m_inverseTransform;
+    SpatialReference m_srs;
+    GDALDataset *m_ds;
 
-    GDALDatasetH m_ds;
     std::string m_errorMsg;
+    mutable std::vector<pdal::Dimension::Type> m_types;
+    std::vector<std::array<double, 2>> m_block_sizes;
 
-private:
     bool getPixelAndLinePosition(double x, double y,
         int32_t& pixel, int32_t& line);
     GDALError computePDALDimensionTypes();
@@ -302,6 +413,5 @@ class PDAL_DLL Raster
 PDAL_DLL std::string transformWkt(std::string wkt, const SpatialReference& from,
     const SpatialReference& to);
 
-
 } // namespace pdal
 

include/pdal/SpatialReference.hpp

@@ -118,10 +118,7 @@ class PDAL_DLL SpatialReference
     // Returns true of OSR can validate the SRS
     bool valid() const;
 
-    // (this is a cleaner way of saying "getWKT() == "")
     /// Returns the OGC WKT describing Spatial Reference System.
-    /// If GDAL is linked, it uses GDAL's operations and methods to determine
-    /// the WKT.  If GDAL is not linked, no WKT is returned.
     /// \param mode_flag May be eHorizontalOnly indicating the WKT will not
     /// include vertical coordinate system info (the default), or
     /// eCompoundOK indicating the the returned WKT may be a compound
@@ -130,8 +127,7 @@ class PDAL_DLL SpatialReference
     std::string getWKT(WKTModeFlag mode_flag = eHorizontalOnly) const;
     std::string getWKT(WKTModeFlag mode_flag, bool pretty) const;
 
-    /// Sets the SRS using GDAL's OGC WKT. If GDAL is not linked, this
-    /// operation has no effect.
+    /// Sets the SRS using GDAL's OGC WKT.
     /// \param v - a string containing the WKT string.
     void setWKT(std::string const& v)
         { m_wkt = v; }

include/pdal/util/Utils.hpp

@@ -906,7 +906,7 @@ namespace Utils
     {
         if (s == "nan" || s == "NaN")
         {
-            d = std::numeric_limits<double>::quiet_NaN();   
+            d = std::numeric_limits<double>::quiet_NaN();
             return true;
         }
 

io/gdal/CMakeLists.txt

@@ -1,10 +1,8 @@
 set(src
     GDALReader.cpp
+    GDALWriter.cpp
+    GDALGrid.cpp
 )
 
-set(inc
-    GDALReader.hpp
-)
-
-PDAL_ADD_DRIVER(reader gdal "${src}" "${inc}" objs)
+PDAL_ADD_DRIVER(reader gdal "${src}" "${incs}" objs)
 set(PDAL_TARGET_OBJECTS ${PDAL_TARGET_OBJECTS} ${objs} PARENT_SCOPE)