diff --git a/src/osgEarth/Math b/src/osgEarth/Math
index 6da88afa99..fc6adaa0c2 100644
--- a/src/osgEarth/Math
+++ b/src/osgEarth/Math
@@ -16,14 +16,13 @@
  * 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/>
  */
-
-#ifndef OSGEARTH_MATH_H
-#define OSGEARTH_MATH_H 1
+#pragma once
 
 #include <osgEarth/Common>
 #include <osgEarth/Geometry>
 #include <osg/Quat>
 #include <osg/Vec3>
+#include <osg/Vec4>
 #include <osg/BoundingBox>
 #include <osg/Array>
 #include <osg/Geometry>
@@ -512,6 +511,12 @@ namespace osgEarth
         return x + 1;
     }
 
+    template<typename T>
+    inline osg::Vec4f mult(const osg::Vec4f& in, T scalar)
+    {
+        return osg::Vec4f(in.x()*scalar, in.y()*scalar, in.z()*scalar, in.w()*scalar);
+    }
+
     // Adapted from Boost - see boost license
     // https://www.boost.org/users/license.html
     template <typename T> inline std::size_t hash_value_unsigned(T val) {
@@ -631,7 +636,4 @@ namespace osgEarth
             const osg::Matrix& m,
             double& L, double& R, double& B, double& T, double& N, double& F);
     };
-
 }
-#endif
-
diff --git a/src/osgEarth/MetaTile b/src/osgEarth/MetaTile
index 924c566148..6349360f47 100644
--- a/src/osgEarth/MetaTile
+++ b/src/osgEarth/MetaTile
@@ -27,6 +27,7 @@
 #include <osgEarth/GeoData>
 #include <osgEarth/TileKey>
 #include <osgEarth/Progress>
+#include <osgEarth/Math>
 #include <osg/Image>
 
 namespace osgEarth { namespace Util
@@ -135,6 +136,10 @@ namespace osgEarth { namespace Util
         inline const typename T::pixel_type* read(double u, double v);
         inline const typename T::pixel_type* read(int s, int t);
 
+        //! Read the interpolated value at the geospatial coordinate (x,y)
+        //! (Bilinear interpolation is used)
+        inline bool readAtCoord(typename T::pixel_type& output, double x, double y);
+
         //! The scale&bias of the tile relative to the key originally passed
         //! to setCenterTileKey
         inline const osg::Matrix& getScaleBias() const {
@@ -170,6 +175,9 @@ namespace osgEarth { namespace Util
         Grid _tiles;
         osg::Matrix _scale_bias; // scale/bias matrix of _centerKey
         unsigned _width, _height;
+
+        typename Tile* getTile(int tx, int ty);
+        bool getTileAndPixelCoords(double u, double v, typename Tile*& tile, double& s, double& t);
     };
 
     template<typename T>
@@ -402,6 +410,117 @@ namespace osgEarth { namespace Util
         return tile._data.read((unsigned)s, (unsigned)t);
     }
 
+    template<typename T>
+    typename MetaTile<T>::Tile* MetaTile<T>::getTile(int tx, int ty)
+    {
+        Tile& tile = _tiles(tx, ty);
+
+        // if we already tried to load this tile and failed, bail out
+        if (tile._failed)
+            return nullptr;
+
+        // if we still need to load this tile, do so
+        if (!tile._data.valid() && _createTile != nullptr)
+        {
+            TileKey key = _centerKey.createNeighborKey(tx, -ty);
+            tile._data = _createTile(key, nullptr);
+            if (!tile._data.valid())
+            {
+                tile._failed = true;
+            }
+        }
+
+        if (tile._failed)
+            return nullptr;
+
+        return &tile;
+    }
+
+    template<typename T>
+    bool MetaTile<T>::readAtCoord(typename T::pixel_type& output, double x, double y)
+    {
+        // geoextents of center tile:
+        double xmin, ymin, xmax, ymax;
+        _centerKey.getExtent().getBounds(xmin, ymin, xmax, ymax);
+
+        // (u, v) ndc relative to center tile:
+        double u = (x - xmin) / (xmax - xmin);
+        double v = (y - ymin) / (ymax - ymin);
+
+        int tx = (int)u;
+        int ty = (int)v;
+
+        // Anchor:
+        typename Tile* p00_tile = getTile(tx, ty);
+        if (!p00_tile)
+            return false;
+        double p00_s = fract(u) * (double)_width, p00_t = fract(v) * (double)_height;
+
+        // Neighbors:
+        typename Tile* p10_tile = p00_tile;
+        double p10_s = p00_s + 1, p10_t = p00_t;
+        if (p10_s >= (double)_width)
+        {
+            p10_tile = getTile(tx + 1, ty);
+            p10_s -= (double)_width;
+        }
+
+        typename Tile* p01_tile = p00_tile;
+        double p01_s = p00_s, p01_t = p00_t + 1;
+        if (p01_t >= _height)
+        {
+            p01_tile = getTile(tx, ty + 1);
+            p01_t -= (double)_height;
+        }
+
+        typename Tile* p11_tile = p00_tile;
+        double p11_s = p00_s + 1, p11_t = p00_t + 1;
+        if (p11_s >= (double)_width || p11_t >= (double)_height)
+        {
+            p11_tile = getTile(tx + 1, ty + 1);
+            if (p11_s > (double)_width) p11_s -= (double)_width;
+            if (p11_t > (double)_height) p11_t -= (double)_height;
+        }
+
+        // calculate weights:
+        unsigned s00 = (unsigned)floor(p00_s);
+        unsigned t00 = (unsigned)floor(p00_t);
+        double left_weight = 1.0 - (p00_s - (double)s00);
+        double bottom_weight = 1.0 - (p00_t - (double)t00);
+
+        if (!p01_tile || !p10_tile || !p11_tile)
+        {
+            return p00_tile->_data.read(output, s00, t00);
+        }
+        else
+        {
+            typename T::pixel_type p00, p10, p01, p11;
+
+            unsigned s01 = clamp((unsigned)floor(p01_s), 0u, _width-1);
+            unsigned t01 = clamp((unsigned)floor(p01_t), 0u, _height-1);
+            unsigned s10 = clamp((unsigned)floor(p10_s), 0u, _width-1);
+            unsigned t10 = clamp((unsigned)floor(p10_t), 0u, _height-1);
+            unsigned s11 = clamp((unsigned)floor(p11_s), 0u, _width-1);
+            unsigned t11 = clamp((unsigned)floor(p11_t), 0u, _height-1);
+
+            if (!p00_tile->_data.read(p00, s00, t00) ||
+                !p01_tile->_data.read(p01, s01, t01) ||
+                !p10_tile->_data.read(p10, s10, t10) ||
+                !p11_tile->_data.read(p11, s11, t11))
+            {
+                return false;
+            }
+
+            // bilinear:
+            typename T::pixel_type p0, p1;
+            p0 = mult(p00, left_weight) + mult(p10, 1.0 - left_weight);
+            p1 = mult(p01, left_weight) + mult(p11, 1.0 - left_weight);
+            output = mult(p0, bottom_weight) + mult(p1, 1.0 - bottom_weight);
+            return true;            
+        }
+    }
+       
+
 } }
 
 #endif // OSGEARTH_METATILE_H
diff --git a/src/osgEarth/XYZ b/src/osgEarth/XYZ
index c9fab2f744..26a70f2ba1 100644
--- a/src/osgEarth/XYZ
+++ b/src/osgEarth/XYZ
@@ -23,6 +23,7 @@
 #include <osgEarth/ElevationLayer>
 #include <osgEarth/URI>
 #include <osgEarth/Threading>
+#include <osgEarth/Containers>
 
 /**
  * XYZ layers. These are general purpose tiled layers that conform
@@ -67,8 +68,8 @@ namespace osgEarth { namespace XYZ
     public:
         META_LayerOptions(osgEarth, XYZImageLayerOptions, ImageLayer::Options);
         OE_OPTION(URI, url);
-        OE_OPTION(bool, invertY);
-        OE_OPTION(std::string, format);
+        OE_OPTION(bool, invertY, false);
+        OE_OPTION(std::string, format, {});
         static Config getMetadata();
         virtual Config getConfig() const;
     private:
@@ -81,9 +82,10 @@ namespace osgEarth { namespace XYZ
     public:
         META_LayerOptions(osgEarth, XYZElevationLayerOptions, ElevationLayer::Options);
         OE_OPTION(URI, url);
-        OE_OPTION(bool, invertY);
-        OE_OPTION(std::string, format);
-        OE_OPTION(std::string, elevationEncoding);
+        OE_OPTION(bool, invertY, false);
+        OE_OPTION(std::string, format, {});
+        OE_OPTION(std::string, elevationEncoding, {});
+        OE_OPTION(bool, stitchEdges, false);
         static Config getMetadata();
         virtual Config getConfig() const;
     private:
@@ -212,6 +214,7 @@ namespace osgEarth
 
     private:
         osg::ref_ptr<XYZImageLayer> _imageLayer;
+        mutable Util::LRUCache<TileKey, GeoImage> _stitchingCache{ true, 64 };
     };
 
 } // namespace osgEarth
diff --git a/src/osgEarth/XYZ.cpp b/src/osgEarth/XYZ.cpp
index 21999a7ddd..7b54635293 100644
--- a/src/osgEarth/XYZ.cpp
+++ b/src/osgEarth/XYZ.cpp
@@ -21,6 +21,7 @@
 #include <osgEarth/FileUtils>
 #include <osgEarth/XmlUtils>
 #include <osgEarth/ImageToHeightFieldConverter>
+#include "MetaTile"
 #include <osgDB/FileUtils>
 
 using namespace osgEarth;
@@ -156,6 +157,7 @@ XYZElevationLayerOptions::getConfig() const
     conf.set("format", _format);
     conf.set("invert_y", _invertY);
     conf.set("elevation_encoding", _elevationEncoding);
+    conf.set("stitch_edges", stitchEdges());
     return conf;
 }
 
@@ -167,6 +169,7 @@ XYZElevationLayerOptions::fromConfig(const Config& conf)
     conf.get("invert_y", _invertY);
     conf.get("elevation_encoding", _elevationEncoding);
     conf.get("interpretation", elevationEncoding()); // compat with QGIS
+    conf.get("stitch_edges", stitchEdges());
 }
 
 Config
@@ -218,6 +221,7 @@ XYZImageLayer::openImplementation()
         return parent;
 
     DataExtentList dataExtents;
+
     Status status = _driver.open(
         options().url().get(),
         options().format().get(),
@@ -233,6 +237,11 @@ XYZImageLayer::openImplementation()
         setProfile(Profile::create("spherical-mercator"));
     }
 
+    if (dataExtents.empty())
+    {
+        dataExtents.push_back(DataExtent(getProfile()->getExtent()));
+    }
+
     setDataExtents(dataExtents);
 
     return Status::NoError;
@@ -308,39 +317,96 @@ XYZElevationLayer::openImplementation()
 GeoHeightField
 XYZElevationLayer::createHeightFieldImplementation(const TileKey& key, ProgressCallback* progress) const
 {
-    // Make an image, then convert it to a heightfield
-    GeoImage geoImage = _imageLayer->createImageImplementation(key, progress);
-    if (geoImage.valid())
+    if (!_imageLayer->mayHaveData(key))
+    {
+        return GeoHeightField::INVALID;
+    }
+
+    std::function<float(const osg::Vec4f&)> decode;
+    
+    if (options().elevationEncoding() == "mapbox")
+    {
+        decode = [](const osg::Vec4f& p) -> float
+            {
+                return -10000.0f + ((p.r() * 255.0f * 256.0f * 256.0f + p.g() * 255.0f * 256.0f + p.b() * 255.0f) * 0.1f);
+            };
+    }
+    else if (options().elevationEncoding() == "terrarium")
+    {
+        decode = [](const osg::Vec4f& p) -> float
+            {
+                return (p.r() * 255.0f * 256.0f + p.g() * 255.0f + p.b() * 255.0f / 256.0f) - 32768.0f;
+            };
+    }
+    else
+    {
+        decode = [](const osg::Vec4f& p) -> float
+            {
+                return p.r();
+            };
+    }
+
+    if (options().stitchEdges() == true)
     {
-        const osg::Image* image = geoImage.getImage();
+        MetaTile<GeoImage> metaImage;
+        metaImage.setCreateTileFunction([this](const TileKey& key, ProgressCallback* progress)
+            {
+                Util::LRUCache<TileKey, GeoImage>::Record r;
+                if (_stitchingCache.get(key, r))
+                {
+                    return r.value();
+                }
+                else
+                {
+                    GeoImage image = _imageLayer->createImage(key, progress);
+                    if (image.valid())
+                    {
+                        _stitchingCache.insert(key, image);
+                    }
+                    return image;
+                }
+            });
+        metaImage.setCenterTileKey(key, progress);
 
-        if (options().elevationEncoding() == "mapbox")
+        if (!metaImage.getCenterTile().valid() || !metaImage.getScaleBias().isIdentity())
         {
-            // Allocate the heightfield.
-            osg::HeightField* hf = new osg::HeightField();
-            hf->allocate(image->s(), image->t());
+            return {};
+        }
 
-            ImageUtils::PixelReader reader(image);
-            osg::Vec4f pixel;
+        auto hf = new osg::HeightField();
+        hf->allocate(getTileSize(), getTileSize());
+        std::fill(hf->getHeightList().begin(), hf->getHeightList().end(), NO_DATA_VALUE);
 
-            for (int c = 0; c < image->s(); c++)
+        double xmin, ymin, xmax, ymax;
+        key.getExtent().getBounds(xmin, ymin, xmax, ymax);
+
+        osg::Vec4f pixel;
+        for (int c = 0; c < getTileSize(); c++)
+        {
+            double u = (double)c / (double)(getTileSize() - 1);
+            double x = xmin + u * (xmax - xmin);
+            for (int r = 0; r < getTileSize(); r++)
             {
-                for (int r = 0; r < image->t(); r++)
+                double v = (double)r / (double)(getTileSize() - 1);
+                double y = ymin + v * (ymax - ymin);
+                if (metaImage.readAtCoord(pixel, x, y))
                 {
-                    reader(pixel, c, r);
-                    pixel.r() *= 255.0;
-                    pixel.g() *= 255.0;
-                    pixel.b() *= 255.0;
-                    float h = -10000.0f + ((pixel.r() * 256.0f * 256.0f + pixel.g() * 256.0f + pixel.b()) * 0.1f);
-                    hf->setHeight(c, r, h);
+                    hf->setHeight(c, r, decode(pixel));
                 }
             }
-
-            return GeoHeightField(hf, key.getExtent());
         }
 
-        else if (options().elevationEncoding() == "terrarium")
+        return GeoHeightField(hf, key.getExtent());
+    }
+
+    else
+    {
+        // Make an image, then convert it to a heightfield
+        GeoImage geoImage = _imageLayer->createImageImplementation(key, progress);
+        if (geoImage.valid())
         {
+            const osg::Image* image = geoImage.getImage();
+
             // Allocate the heightfield.
             osg::HeightField* hf = new osg::HeightField();
             hf->allocate(image->s(), image->t());
@@ -353,25 +419,13 @@ XYZElevationLayer::createHeightFieldImplementation(const TileKey& key, ProgressC
                 for (int r = 0; r < image->t(); r++)
                 {
                     reader(pixel, c, r);
-                    pixel.r() *= 255.0;
-                    pixel.g() *= 255.0;
-                    pixel.b() *= 255.0;
-                    float h = (pixel.r() * 256.0f + pixel.g() + pixel.b() / 256.0f) - 32768.0f;
-                    hf->setHeight(c, r, h);
+                    hf->setHeight(c, r, decode(pixel));
                 }
             }
 
             return GeoHeightField(hf, key.getExtent());
         }
 
-        else
-        {
-            ImageToHeightFieldConverter conv;
-            osg::HeightField* hf = conv.convert(image);
-            return GeoHeightField(hf, key.getExtent());
-        }
+        return GeoHeightField::INVALID;
     }
-    
-
-    return GeoHeightField::INVALID;
 }
diff --git a/tests/mapzen.earth b/tests/mapzen.earth
new file mode 100644
index 0000000000..190f2de93e
--- /dev/null
+++ b/tests/mapzen.earth
@@ -0,0 +1,28 @@
+<Map name="MapZen Elevation">
+    <options>
+        <terrain skirt_ratio="0.05"/>
+    </options>
+    
+    <XYZElevation name="MapZen global elevation (S3)">
+        <url>https://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png</url>
+        <elevation_encoding>terrarium</elevation_encoding>
+        <stitch_edges>true</stitch_edges>
+    </XYZElevation>
+
+    <!-- The stops are optional. Leave them out for a default. -->
+    <ContourMap>
+    <!--
+        <stops>
+            <stop elevation="-3000" color="#00002f"/>
+            <stop elevation="-2500" color="#00007f"/>
+            <stop elevation="-625" color="#0000ff"/>
+            <stop elevation="0" color="#007fff"/>
+            <stop elevation="3" color="#c2b280"/>
+            <stop elevation="160" color="#1f8f00"/>
+            <stop elevation="850" color="#bfaf72"/>
+            <stop elevation="1875" color="#7f7f7f"/>
+            <stop elevation="2500" color="#ffffff"/>
+        </stops>
+        -->
+    </ContourMap>
+</Map>