/* -*-c++-*- */
/* osgEarth - Geospatial SDK for OpenSceneGraph
 * Copyright 2020 Pelican Mapping
 * http://osgearth.org
 *
 * osgEarth is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>
 */
#pragma once

#include <osgEarth/Common>
#include <osgEarth/GeoCommon>
#include <osgEarth/Bounds>
#include <osgEarth/SpatialReference>
#include <osgEarth/ImageUtils>
#include <osgEarth/Status>
#include <osgEarth/Threading>

#include <osg/Referenced>
#include <osg/Image>
#include <osg/Shape>

#include <unordered_map>

namespace osg
{
    class Camera;
    class Polytope;
}

namespace osgEarth
{
    class TerrainResolver;
    class GeoExtent;

    /**
     * A georeferenced 3D point.
     */
    class OSGEARTH_EXPORT GeoPoint
    {
    public:

        /**
         * Constructs a GeoPoint.
         */
        GeoPoint(
            const SpatialReference* srs,
            double x,
            double y,
            double z,
            const AltitudeMode& mode );

        /**
         * Constructs a GeoPoint with an absolute Z.
         */
        GeoPoint(
            const SpatialReference* srs,
            double x,
            double y,
            double z );

        /**
         * Constructs a GeoPoint with X and Y coordinates. The Z defaults
         * to zero with an ALTMODE_RELATIVE altitude mode (i.e., 0 meters
         * above the terrain).
         */
        GeoPoint(
            const SpatialReference* srs,
            double x,
            double y );

        /**
         * Constructs a GeoPoint from a vec3.
         */
        GeoPoint(
            const SpatialReference* srs,
            const osg::Vec3d&       xyz,
            const AltitudeMode&     mode );

        /**
         * Constructs a GeoPoint from a vec3 with absolute Z.
         */
        GeoPoint(
            const SpatialReference* srs,
            const osg::Vec3d&       xyz);

        /**
         * Constructs a new GeoPoint by transforming an existing GeoPoint into
         * the specified spatial reference.
         */
        GeoPoint(
            const SpatialReference* srs,
            const GeoPoint&         rhs );

        /**
         * Copy constructor
         */
        GeoPoint(const GeoPoint& rhs);

        /**
         * Constructs an empty (and invalid) geopoint.
         */
        GeoPoint();
        
        /**
         * Constructs a GeoPoint at 0,0,0 absolute.
         */
        GeoPoint(const SpatialReference* srs);

        /**
         * Constructs a geopoint from serialization
         */
        GeoPoint( const Config& conf, const SpatialReference* srs =0L );

        /** dtor */
        virtual ~GeoPoint() { }

        /**
         * Sets the SRS and coords 
         */
        void set(
            const SpatialReference* srs,
            const osg::Vec3d&       xyz,
            const AltitudeMode&     mode );

        void set(
            const SpatialReference* srs,
            double                  x,
            double                  y,
            double                  z,
            const AltitudeMode&     mode );

        // component getter/setters

        double& x() { return _p.x(); }
        double  x() const { return _p.x(); }

        double& y() { return _p.y(); }
        double  y() const { return _p.y(); }

        double& z() { return _p.z(); }
        double  z() const { return _p.z(); }

        double& alt() { return _p.z(); }
        double  alt() const { return _p.z(); }

        osg::Vec3d& vec3d() { return _p; }
        const osg::Vec3d& vec3d() const { return _p; }

        const SpatialReference* getSRS() const { return _srs.get(); }

        /**
         * AltitudeMode reflects whether the Z coordinate is absolute with
         * respect to MSL or relative to the terrain elevation at that
         * point. When using relative points, GeoPoint usually requires
         * access to a TerrainProvider in order to resolve the altitude.
         */
        AltitudeMode& altitudeMode() { return _altMode; }
        const AltitudeMode& altitudeMode() const { return _altMode; }
        bool isRelative() const { return _altMode == ALTMODE_RELATIVE; }
        bool isAbsolute() const { return _altMode == ALTMODE_ABSOLUTE; }

        /**
         * Returns a copy of this geopoint transformed into another SRS.
         */
        GeoPoint transform(const SpatialReference* outSRS) const;

        /**
         * Transforms this geopoint into another SRS.
         */
        bool transform(const SpatialReference* outSRS, GeoPoint& output) const;

        /**
         * Transforms this point in place to another SRS.
         */
        bool transformInPlace(const SpatialReference* srs);

        /**
         * Transforms this geopoint's Z coordinate (in place)
         */
        bool transformZ(const AltitudeMode& altMode, const TerrainResolver* t );

        /**
         * Transforms and returns the geopoints Z coordinate.
         */
        bool transformZ(const AltitudeMode& altMode, const TerrainResolver* t, double& out_z) const;

        //! Transforms a resolution distance to a cartesian value.
        Distance transformResolution(const Distance& d, const UnitsType& outUnits) const;

        /**
         * Transforms this geopoint's Z to be absolute w.r.t. the vertical datum
         */
        bool makeAbsolute(const TerrainResolver* t) { return transformZ(ALTMODE_ABSOLUTE, t); }

        /**
         * Transforms this geopoint's Z to be terrain-relative.
         */
        bool makeRelative(const TerrainResolver* t) { return transformZ(ALTMODE_RELATIVE, t); }

        /**
         * Transforms this GeoPoint to geographic (lat/long) coords in place.
         */
        bool makeGeographic();

        /**
         * Outputs world coordinates corresponding to this point. If the point
         * is ALTMODE_RELATIVE, this will fail because there's no way to resolve
         * the actual Z coordinate. Use the variant of toWorld that takes a
         * Terrain* instead.
         */
        bool toWorld( osg::Vec3d& out_world ) const;

        /**
         * Outputs world coordinates corresponding to this point, passing in a Terrain
         * object that will be used if the point needs to be converted to absolute
         * altitude
         */
        bool toWorld( osg::Vec3d& out_world, const TerrainResolver* terrain ) const;

        /**
         * Converts world coordinates into a geopoint
         */
        bool fromWorld(const SpatialReference* srs, const osg::Vec3d& world);

        /**
         * geopoint into absolute world coords.
         */
        bool createLocalToWorld( osg::Matrixd& out_local2world ) const;

        /**
         * Outputs a matrix that will transform absolute world coordiantes so they are
         * localized into a local tangent place around this geopoint.
         */
        bool createWorldToLocal( osg::Matrixd& out_world2local ) const;

        /**
         * Converts this point to the same point in a local tangent plane.
         */
        GeoPoint toLocalTangentPlane() const;

        /**
         * Outputs an "up" vector corresponding to the given point. The up vector
         * is orthogonal to a local tangent plane at that point on the map.
         */
        bool createWorldUpVector( osg::Vec3d& out_up ) const;

        //! Geodesic distance from this point to another.
        //! This is the distance along the real-world ellipsoidal surface
        //! of the Earth (or other body), regardless of map projection.
        //! It does not account for Z/altitude.
        //! @param rhs Other point
        //! @return Geodesic distance between the two points
        Distance geodesicDistanceTo(const GeoPoint& rhs) const;

        /**
         * @deprecated - ambiguous, will be removed. Use geodesicDistanceTo() or toWorld()/length instead.
         * Calculates the distance in meters from this geopoint to another.
         */
        double distanceTo(const GeoPoint& rhs) const;

        /**
         * Interpolates a point between this point and another point
         * using the parameter t [0..1].
         */
        GeoPoint interpolate(const GeoPoint& rhs, double t) const;

        //! Convert this point to screen coordinates.
        osg::Vec2d toScreen(const osg::Camera* camera) const;

        //! Convenience function to return xy units
        const UnitsType& getXYUnits() const;

        bool operator == (const GeoPoint& rhs) const;
        bool operator != (const GeoPoint& rhs) const { return !operator==(rhs); }
        bool isValid() const { return _srs.valid(); }

        Config getConfig() const;

        /**
         * Represent this point as a string
         */
        std::string toString() const;

    public:
        static GeoPoint INVALID;

    protected:
        osg::Vec3d _p;
        osg::ref_ptr<const SpatialReference> _srs;
        AltitudeMode _altMode;
    };


    /**
     * A simple circular bounding area consiting of a GeoPoint and a linear radius.
     */
     class OSGEARTH_EXPORT GeoCircle
     {
     public:
         /** Construct an INVALID GeoCircle */
        GeoCircle();

        /** Copy another GoeCircle */
        GeoCircle(const GeoCircle& rhs);

        /** Construct a new GeoCircle */
        GeoCircle(
             const GeoPoint& center,
             double          radius );

        virtual ~GeoCircle() { }

        /** The center point of the circle */
        const GeoPoint& getCenter() const { return _center; }
        void setCenter( const GeoPoint& value ) { _center = value; }

        /** Circle's radius, in linear map units (or meters for a geographic SRS) */
        double getRadius() const { return _radius; }
        void setRadius( double value ) { _radius = value; }

        /** SRS of the center point */
        const SpatialReference* getSRS() const { return _center.getSRS(); }

        /** equality test */
        bool operator == ( const GeoCircle& rhs ) const;

        /** inequality test */
        bool operator != ( const GeoCircle& rhs ) const { return !operator==(rhs); }

        /** validity test */
        bool isValid() const { return _center.isValid() && _radius > 0.0; }

        /** transform the GeoCircle to another SRS */
        GeoCircle transform( const SpatialReference* srs ) const;

        /** transform the GeoCircle to another SRS */
        bool transform( const SpatialReference* srs, GeoCircle& out_circle ) const;

        /** does this GeoCircle intersect another? */
        bool intersects( const GeoCircle& rhs ) const;

     public:

         static GeoCircle INVALID;

     protected:
         GeoPoint _center;
         double   _radius;
     };
}
OSGEARTH_SPECIALIZE_CONFIG(osgEarth::GeoPoint);

namespace osgEarth
{
    /**
     * An axis-aligned geospatial extent. A bounding box that is aligned with a
     * spatial reference's coordinate system.
     */
    class OSGEARTH_EXPORT GeoExtent
    {
    public:
        /** Default ctor creates an "invalid" extent */
        GeoExtent();

        /** Contructs a valid extent */
        GeoExtent(
            const SpatialReference* srs,
            double west, double south,
            double east, double north);

        /** Contructs an invalid extent that you can grow with the expandToInclude method */
        GeoExtent(const SpatialReference* srs);

        /** Copy ctor */
        GeoExtent(const GeoExtent& rhs);

        /** create from Bounds object */
        GeoExtent(const SpatialReference* srs, const Bounds& bounds);

        /** dtor */
        virtual ~GeoExtent() { }

        //! Set from the SW and NE corners.
        void set(double west, double south, double east, double north);

        bool operator == (const GeoExtent& rhs) const;
        bool operator != (const GeoExtent& rhs) const;

        /** Gets the spatial reference system underlying this extent. */
        const SpatialReference* getSRS() const { return _srs.get(); }

        //! Coordinates of the bounding edges, normalized for the lat/long frame if necessary
        inline double west() const;
        inline double east() const;
        inline double south() const;
        inline double north() const;

        //! Coordinates of the bounds, NOT normalized to the lat/long frame.
        inline double xMin() const { return _west; }
        inline double xMax() const { return _west + _width; }
        inline double yMin() const { return _south; }
        inline double yMax() const { return _south + _height; }

        //! East-to-west span of the extent
        inline double width() const { return _width; }

        //! East-to-est span in specified units
        double width(const UnitsType& units) const;

        //! North-to-south span of the extent
        inline double height() const { return _height; }

        //! North-to-south span in specified units
        double height(const UnitsType& units) const;

        //! Gets the centroid of the bounds
        GeoPoint getCentroid() const;

        //! Legacy @deprecated
        bool getCentroid(double& out_x, double& out_y) const;

        //! True if the extent is geographic and crosses the 180 degree meridian.
        bool crossesAntimeridian() const;

        //! Raw bounds of the extent (unnormalized)
        void getBounds(double& xmin, double& ymin, double& xmax, double& ymax) const;

        /** True if this object defines a real, valid extent with positive area */
        inline bool isValid() const
        {
            return _srs.valid() && _width >= 0.0 && _height >= 0.0;
        }

        inline bool isInvalid() const { return !isValid(); }

        /**
         * If this extent crosses the international date line, populates two extents, one for
         * each side, and returns true. Otherwise returns false and leaves the reference
         * parameters untouched.
         */
        bool splitAcrossAntimeridian(GeoExtent& first, GeoExtent& second) const;

        /**
         * Returns this extent transformed into another spatial reference.
         *
         * NOTE! It is possible that the target SRS will not be able to accomadate the
         * extents of the source SRS. (For example, transforming a full WGS84 extent
         * to Mercator will resultin an error since Mercator does not cover the entire
         * globe.) Consider using Profile:clampAndTransformExtent() instead of using
         * this method directly.
         */
        GeoExtent transform(const SpatialReference* to_srs) const;

        /**
         * Same as transform(srs) but puts the result in the output extent
         */
        bool transform(const SpatialReference* to_srs, GeoExtent& output) const;

        /**
         * Returns true if the specified point falls within the bounds of the extent.
         *
         * @param x, y
         *      Coordinates to test
         * @param xy_srs
         *      SRS of input x and y coordinates; if null, the method assumes x and y
         *      are in the same SRS as this object.
         */
        bool contains(double x, double y, const SpatialReference* srs = 0L) const;

        template<typename VEC>
        bool contains(const VEC& a, const SpatialReference* srs = nullptr) const { 
            return contains(a.x(), a.y(), srs);
        }

        /**
         * Returns true if the point falls within this extent.
         */
        bool contains(const GeoPoint& rhs) const;

        /**
         * Returns true if this extent fully contains another extent.
         */
        bool contains(const GeoExtent& rhs) const;

        /**
         * Returns true if this extent fully contains the target bounds.
         */
        bool contains(const Bounds& rhs) const;

        /**
         * Returns TRUE if this extent intersects another extent.
         * @param[in ] rhs      Extent against which to perform intersect test
         * @param[in ] checkSRS If false, assume the extents have the same SRS (don't check).
         */
        bool intersects(const GeoExtent& rhs, bool checkSRS = true) const;

        /**
         * Copy of the anonymous bounding box
         */
        Bounds bounds() const;

        /**
         * Gets a geo circle bounding this extent.
         */
        GeoCircle computeBoundingGeoCircle() const;

        /**
         * Grow this extent to include the specified point (which is assumed to be
         * in the extent's SRS.
         */
        void expandToInclude(double x, double y);
        void expandToInclude(const osg::Vec3d& v) { expandToInclude(v.x(), v.y()); }
        //void expandToInclude( const Bounds& bounds );

        /**
         * Expand this extent to include the bounds of another extent.
         */
        bool expandToInclude(const GeoExtent& rhs);

        /**
         * Intersect this extent with another extent in the same SRS and return the
         * result.
         */
        GeoExtent intersectionSameSRS(const GeoExtent& rhs) const; //const Bounds& rhs ) const;

        /**
         * Returns a human-readable string containing the extent data (without the SRS)
         */
        std::string toString() const;

        /**
         *Inflates this GeoExtent by the given ratios
         */
        void scale(double x_scale, double y_scale);

        /**
         * Expands the extent by x and y.
         */
        void expand(double x, double y);

        /**
        * Expands the extent by x and y.
        */
        void expand(const Distance& x, const Distance& y);

        /**
         *Gets the area of this GeoExtent
         */
        double area() const;

        /**
         * Generate a polytope in world coordinates that bounds the extent.
         * Return false if the extent it invalid.
         */
        bool createPolytope(osg::Polytope& out) const;

        /**
         * Computes a scale/bias matrix that transforms parametric coordinates [0..1]
         * from this extent into the target extent. Return false if the extents are
         * incompatible (different SRS, etc.)
         *
         * For example, if this extent is 100m wide, and the target extent is
         * 200m wide, the output matrix will have an x_scale = 0.5.
         *
         * Note! For the sake of efficiency, this method does NOT check for
         * validity nor for SRS equivalence -- so be sure to validate those beforehand.
         * It also assumes the output matrix is preset to the identity.
         */
        bool createScaleBias(const GeoExtent& target, osg::Matrixd& output) const;
        bool createScaleBias(const GeoExtent& target, osg::Matrixf& output) const;

        /**
         * Generates a BoundingSphere encompassing the extent and a vertical
         * volume, in world coordinates.
         */
        osg::BoundingSphered createWorldBoundingSphere(double minElev, double maxElev) const;

        /**
         * Returns true if the extent is the entire Earth, for geographic
         * extents. Otherwise returns false.
         */
        bool isWholeEarth() const;

    public:
        static GeoExtent INVALID;

    public: // config
        Config getConfig() const;
        void fromConfig(const Config& conf);

    private:
        double _west, _width, _south, _height;
        osg::ref_ptr<const SpatialReference> _srs;

        double normalizeX(double longitude) const;

        void clamp();

        bool isGeographic() const;

        void setOriginAndSize(double west, double south, double width, double height);
    };

    inline double GeoExtent::west() const { return _west; }

    // if east is exactly on the antimeridian, return 180.0
    inline double GeoExtent::east() const {
        double a = normalizeX(_west + _width);
        return (isValid() && _srs->isGeographic() && a == -180.0) ? 180.0 : a;
    }

    inline double GeoExtent::south() const { return _south; }

    inline double GeoExtent::north() const { return _south + _height; }


    /**
     * A geospatial area with tile data LOD extents
     */
    class OSGEARTH_EXPORT DataExtent : public GeoExtent
    {
    public:
        DataExtent();
        DataExtent(const GeoExtent& extent);
        DataExtent(const GeoExtent& extent, const std::string &description);
        DataExtent(const GeoExtent& extent, unsigned minLevel);
        DataExtent(const GeoExtent& extent, unsigned minLevel, const std::string &description);
        DataExtent(const GeoExtent& extent, unsigned minLevel, unsigned maxLevel);
        DataExtent(const GeoExtent& extent, unsigned minLevel, unsigned maxLevel, const std::string &description);

        /** dtor */
        virtual ~DataExtent() { }

        /** The minimum LOD of the extent */
        optional<unsigned>& minLevel() { return _minLevel; }
        const optional<unsigned>& minLevel() const { return _minLevel; }

        /** The maximum LOD of the extent */
        optional<unsigned>& maxLevel() { return _maxLevel; }
        const optional<unsigned>& maxLevel() const { return _maxLevel; }

        /** description for the data extents */
        const osgEarth::optional<std::string>& description() const { return _description; }

    private:
        optional<unsigned> _minLevel;
        optional<unsigned> _maxLevel;
        optional<std::string> _description;
    };

    typedef std::vector< DataExtent > DataExtentList;


    /**
     * A georeferenced image; i.e. an osg::Image and an associated GeoExtent with SRS.
     */
    class OSGEARTH_EXPORT GeoImage
    {
    public:
        using pixel_type = osg::Vec4f;

        /** Construct an empty (invalid) geoimage. */
        GeoImage();

        //! Construct an image with an error status
        GeoImage(const Status&);

        //! Constructs a new goereferenced image.
        GeoImage(const osg::Image* image, const GeoExtent& extent);

        //! Constructs a new goereferenced image from a future.
        GeoImage(jobs::future<osg::ref_ptr<osg::Image>> image, const GeoExtent& extent);

        /** dtor */
        virtual ~GeoImage() { }

        static GeoImage INVALID;

    public:
        /**
         * True if this is a valid geo image.
         */
        bool valid() const;

        //! Error status if set
        const Status& getStatus() const { return _status; }

        /**
         * Gets a pointer to the underlying OSG image.
         */
        const osg::Image* getImage() const;

        /**
         * Gets the geospatial extent of the image.
         */
        const GeoExtent& getExtent() const;

        /**
         * Shortcut to get the spatial reference system describing
         * the projection of the image.
         */
        const SpatialReference* getSRS() const;

        /**
         * Crops the image to a new geospatial extent.
         *
         * @param extent
         *      New extent to which to crop the image.
         * @param exact
         *      If "exact" is true, the output image will have exactly the extents requested;
         *      this process may require resampling and will therefore be more expensive. If
         *      "exact" is false, we do a simple crop of the image that is rounded to the nearest
         *      pixel. The resulting extent will be close but usually not exactly what was
         *      requested - however, this method is faster.
         * @param width, height
         *      New pixel size for the output image. By default, the method will automatically
         *      calculate a new pixel size.
         */
        GeoImage crop(
            const GeoExtent& extent,
            bool exact = false,
            unsigned int width = 0,
            unsigned int height = 0,
            bool useBilinearInterpolation = true) const;

        /**
         * Warps the image into a new spatial reference system.
         *
         * @param to_srs
         *      SRS into which to warp the image.
         * @param to_extent
         *      Supply this extent if you wish to warp AND crop the image in one step. This is
         *      faster than calling reproject() and then crop().
         * @param width, height
         *      New pixel size for the output image. Be default, the method will automatically
         *      calculate a new pixel size.
         */
        GeoImage reproject(
            const SpatialReference* to_srs,
            const GeoExtent* to_extent = 0,
            unsigned int width = 0,
            unsigned int height = 0,
            bool useBilinearInterpolation = true) const;

        /**
         * Returns the underlying OSG image and releases the reference pointer.
         */
        osg::ref_ptr<osg::Image> takeImage();

        /**
         * Gets the units per pixel of this geoimage
         */
        double getUnitsPerPixel() const;

        //! Gets the coordinate at the image's s,t
        bool getCoord(int s, int t, double& out_x, double& out_y) const;

        //! Sets a token object in this GeoImage. You can use this to
        //! keep a weak reference to the object after creating it,
        //! for example to detect destruction.
        void setTrackingToken(osg::Object* token);
        osg::Object* getTrackingToken() const;


        unsigned s() const {
            return _myimage.valid() ? _myimage->s() : 0;
        }
        unsigned t() const {
            return _myimage.valid() ? _myimage->t() : 0;
        }
        bool read(pixel_type& output, unsigned s, unsigned t) const {
            _read(output, s, t);
            return output.r() != NO_DATA_VALUE;
        }
        bool read(pixel_type& output, double u, double v) const {
            _read(output, u, v);
            return output.r() != NO_DATA_VALUE;
        }

        ImageUtils::PixelReader& getReader() { return _read; };
        const ImageUtils::PixelReader& getReader() const { return _read; }

        //! Read the value of a pixel at a geopoint.
        bool read(
            pixel_type& output,
            const GeoPoint& p) const;

    private:
        GeoExtent _extent;
        Status _status;
        osg::ref_ptr<const osg::Image> _myimage;
        mutable optional<jobs::future<osg::ref_ptr<osg::Image>>> _future;
        osg::ref_ptr<osg::Object> _token;
        ImageUtils::PixelReader _read;
    };

    typedef std::vector<GeoImage> GeoImageVector;

    struct GeoImageIterator : public ImageUtils::ImageIteratorWithExtent<GeoExtent>
    {
        GeoImageIterator(const GeoImage& im) :
            ImageUtils::ImageIteratorWithExtent<GeoExtent>(
                im.getImage(), im.getExtent()) { }
    };

    struct GeoImagePixelReader : public ImageUtils::PixelReaderWithExtent<GeoExtent>
    {
        GeoImagePixelReader(const GeoImage& im) :
            ImageUtils::PixelReaderWithExtent<GeoExtent>(
                im.getImage(), im.getExtent()) { }
    };

    /**
     * A georeferenced heightfield and associated normal/curvature map.
     */
    class OSGEARTH_EXPORT GeoHeightField
    {
    public:
        /** Constructs an empty (invalid) heightfield. */
        GeoHeightField();

        //! Construct a heightfield with an error status
        GeoHeightField(const Status&);

        /**
         * Constructs a new georeferenced heightfield.
         */
        GeoHeightField(
            const osg::HeightField* heightField,
            const GeoExtent&  extent);

        /** dtor */
        virtual ~GeoHeightField() { }

        static GeoHeightField INVALID;

        /**
         * True if this is a valid heightfield.
         */
        bool valid() const;

        //! Status indicator
        const Status& getStatus() const { return _status; }

        /**
         * Gets the elevation value at a specified point.
         *
         * @param srs
         *      Spatial reference of the query coordinates. (If you pass in NULL, the method
         *      will assume that the SRS is equivalent to that of the GeoHeightField. Be sure
         *      this is case of you will get incorrect results.)
         * @param x, y
         *      Coordinates at which to query the elevation value.
         * @param interp
         *      Interpolation method for the elevation query.
         * @param srsWithOutputVerticalDatum
         *      Transform the output elevation value according to the vertical datum
         *      associated with this SRS. If the SRS is NULL, assume a geodetic vertical datum
         *      relative to this object's reference ellipsoid.
         * @param out_elevation
         *      Output: the elevation value
         * @return
         *      True if the elevation query was succesful; false if not (e.g. if the query
         *      fell outside the geospatial extent of the heightfield)
         */
        bool getElevation(
            const SpatialReference* inputSRS,
            double                  x,
            double                  y,
            RasterInterpolation     interp,
            const SpatialReference* srsWithOutputVerticalDatum,
            float&                  out_elevation) const;

        //! Gets the elevation at a point (must be in the same SRS; bilinear interpolation)
        float getElevation(double x, double y, RasterInterpolation interp = INTERP_BILINEAR) const;

        /**
         * Subsamples the heightfield, returning a new heightfield corresponding to
         * the destEx extent. The destEx must be a smaller, inset area of sourceEx.
         */
        GeoHeightField createSubSample(const GeoExtent& destEx, unsigned int width, unsigned int height, RasterInterpolation interpolation) const;

        /**
         * Gets the geospatial extent of the heightfield.
         */
        const GeoExtent& getExtent() const;

        /**
         * The minimum height in the heightfield
         */
        float getMinHeight() const { return _minHeight; }

        /**
         * The maximum height in the heightfield
         */
        float getMaxHeight() const { return _maxHeight; }

        /**
         * Gets a pointer to the underlying OSG heightfield.
         */
        const osg::HeightField* getHeightField() const;

        /**
         * Gets the X interval of this GeoHeightField
         */
        double getXInterval() const;

        /**
         * Gets the Y interval of this GeoHeightField
         */
        double getYInterval() const;


        //Sort GeoHeightField's by resolution
        struct SortByResolutionFunctor
        {
            inline bool operator() (const GeoHeightField& lhs, const GeoHeightField& rhs) const
            {
                return lhs.getXInterval() < rhs.getXInterval();
            }
        };

    protected:
        GeoExtent _extent;
        Status _status;
        float _minHeight, _maxHeight;
        osg::ref_ptr<const osg::HeightField> _heightField;

        void init();
    };

    typedef std::vector<GeoHeightField> GeoHeightFieldVector;


    /**
    * A georeferenced node.
    */
    class OSGEARTH_EXPORT GeoNode
    {
    public:
        //! Construct an empty (invalid) GeoNode
        GeoNode();

        //! Construct an image with an error status
        GeoNode(const Status&);

        //! Constructs a new goereferenced node
        GeoNode(const osg::Node* node, const GeoExtent& extent);

        //! Global invalid node
        static GeoNode INVALID;

    public:
        //! True if this is a valid object
        bool valid() const;

        //! Error status if set
        const Status& getStatus() const { return _status; }

        //! The underlying OSG node
        const osg::Node* getNode() const { return _node.get(); }

        //! The node's extent
        const GeoExtent& getExtent() const { return _extent; }

    private:
        GeoExtent _extent;
        Status _status;
        osg::ref_ptr<const osg::Node> _node;
    };

    typedef std::vector<GeoNode> GeoNodeVector;
}