Update contourmap to use MSL heights by default

src/osgEarth/ContourMap

@@ -19,8 +19,7 @@
 * 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 OSGEARTHUTIL_CONTOUR_MAP_H
-#define OSGEARTHUTIL_CONTOUR_MAP_H
+#pragma once
 
 #include <osgEarth/VisibleLayer>
 #include <osgEarth/Color>
@@ -44,6 +43,7 @@ namespace osgEarth
                 Color color;
             };
             OE_OPTION_VECTOR(Stop, stops);
+            OE_OPTION(std::string, vdatum, "egm96");
             virtual Config getConfig() const;
             static Config getMetadata();
         private:
@@ -87,16 +87,14 @@ namespace osgEarth
         typedef osg::Texture1D TextureType;
 #endif
 
-        TextureImageUnitReservation           _reservation;
+        TextureImageUnitReservation _reservationColorRamp;
+        TextureImageUnitReservation _reservationGeoid;
         osg::ref_ptr<osg::TransferFunction1D> _xfer;
-        osg::ref_ptr<TextureType>             _xferTexture;
-        osg::ref_ptr<osg::Uniform>            _xferSampler;
-        osg::ref_ptr<osg::Uniform>            _xferMin;
-        osg::ref_ptr<osg::Uniform>            _xferRange;
+        osg::ref_ptr<TextureType> _xferTexture;
+        osg::ref_ptr<osg::Uniform> _xferSampler;
+        osg::ref_ptr<osg::Uniform> _xferMin;
+        osg::ref_ptr<osg::Uniform> _xferRange;
     };
 
 }
 OSGEARTH_SPECIALIZE_CONFIG(osgEarth::ContourMapLayer::Options);
-
-
-#endif // OSGEARTHUTIL_CONTOUR_MAP_H

src/osgEarth/ContourMap.cpp

@@ -24,6 +24,8 @@
 #include <osgEarth/Registry>
 #include <osgEarth/TerrainEngineNode>
 #include <osgEarth/MapNode>
+#include <osgEarth/VerticalDatum>
+#include <osgDB/WriteFile>
 
 #define LC "[ContourMap] "
 
@@ -56,6 +58,7 @@ ContourMapLayer::Options::getConfig() const
     }
     if (!stopsConf.empty())
         conf.add(stopsConf);
+    conf.set("vdatum", vdatum());
     return conf;
 }
 
@@ -70,6 +73,7 @@ ContourMapLayer::Options::fromConfig(const Config& conf)
         if (stop.get("elevation", s.elevation) && stop.get("color", s.color))
             stops().emplace_back(s);
     }
+    conf.get("vdatum", vdatum());
 }
 
 //........................................................................
@@ -133,6 +137,7 @@ ContourMapLayer::init()
 
     // build a transfer function.
     osg::TransferFunction1D* xfer = new osg::TransferFunction1D();
+    xfer->allocate(1024);
 
     if (options().stops().empty() == false)
     {
@@ -146,8 +151,8 @@ ContourMapLayer::init()
         xfer->setColor(-7500.0f, Color("#000016"));
         xfer->setColor(-2500.0f, Color("#00007f"));
         xfer->setColor(-625.0f, Color("#0000ff"));
-        xfer->setColor(0.0f, Color("#007fff"));
-        xfer->setColor(2.5f, Color("#c2b280"));
+        xfer->setColor(-1.0f, Color("#007fff"));
+        xfer->setColor(0.0f, Color("#c2b280"));
         xfer->setColor(50.0f, Color("#d6d63f"));
         xfer->setColor(250.0f, Color("#1f9e00"));
         xfer->setColor(1000.0f, Color("#30632c"));
@@ -166,21 +171,59 @@ ContourMapLayer::prepareForRendering(TerrainEngine* engine)
 {
     super::prepareForRendering(engine);
 
-    if (!engine->getResources()->reserveTextureImageUnitForLayer(_reservation, this, "ContourMap"))
+    if (!engine->getResources()->reserveTextureImageUnitForLayer(_reservationColorRamp, this, "ContourMap"))
     {
         setStatus(Status::ResourceUnavailable, "No texture image units available");
         return;
     }
 
     osg::StateSet* stateset = getOrCreateStateSet();
-    stateset->setTextureAttributeAndModes(_reservation.unit(), _xferTexture.get(), osg::StateAttribute::ON);
-    _xferSampler->set(_reservation.unit());
+    stateset->setTextureAttribute(_reservationColorRamp.unit(), _xferTexture.get());
+    _xferSampler->set(_reservationColorRamp.unit());
+
+    if (!options().vdatum()->empty())
+    {
+        auto* vdatum = VerticalDatum::get(options().vdatum().get());
+        if (vdatum)
+        {
+            if (!engine->getResources()->reserveTextureImageUnitForLayer(_reservationGeoid, this, "ContourMap"))
+            {
+                setStatus(Status::ResourceUnavailable, "No texture image units available");
+                return;
+            }
+
+            auto* geoid = vdatum->getGeoid();
+            OE_SOFT_ASSERT_AND_RETURN(geoid, void());
+
+            auto* hf = geoid->getHeightField();
+            OE_SOFT_ASSERT_AND_RETURN(hf, void());
+
+            auto image = new osg::Image();
+            image->allocateImage(hf->getNumColumns(), hf->getNumRows(), 1, GL_RED, GL_FLOAT);
+            memcpy(image->data(), hf->getFloatArray()->getDataPointer(), hf->getFloatArray()->getTotalDataSize());
+
+            auto tex = new osg::Texture2D(image);
+            tex->setInternalFormat(GL_R32F);
+            tex->setResizeNonPowerOfTwoHint(false);
+            tex->setWrap(osg::Texture::WRAP_S, osg::Texture::CLAMP_TO_EDGE);
+            tex->setWrap(osg::Texture::WRAP_T, osg::Texture::CLAMP_TO_EDGE);
+            tex->setFilter(osg::Texture::MIN_FILTER, osg::Texture::NEAREST);
+            tex->setFilter(osg::Texture::MAG_FILTER, osg::Texture::NEAREST);
+
+            stateset->setTextureAttribute(_reservationGeoid.unit(), tex);
+            stateset->addUniform(new osg::Uniform("oe_contour_geoid", _reservationGeoid.unit()));
+            stateset->setDefine("OE_USE_GEOID");
+
+            osgDB::writeImageFile(*image, "geoid.tif");
+        }
+    }
 }
 
 Status
 ContourMapLayer::closeImplementation()
 {
-    _reservation.release();
+    _reservationColorRamp.release();
+    _reservationGeoid.release();
 
     return super::closeImplementation();
 }

src/osgEarth/ContourMap.glsl

@@ -1,32 +1,59 @@
 #pragma vp_function oe_contour_vertex, vertex_view
 #pragma import_defines(OE_USE_GL4)
+#pragma import_defines(OE_USE_GEOID)
+
+out vec4 oe_layer_tilec;
+
+#ifdef OE_USE_GEOID
+  vec4 oe_tile_key;
+  uniform sampler2D oe_contour_geoid;
+  out float oe_contour_offset;
+#endif
 
 #ifdef OE_USE_GL4
 #pragma include RexEngine.GL4.glsl
-
 uint64_t oe_terrain_getElevationHandle();
 vec2 oe_terrain_getElevationCoord(in vec2);
-out vec4 oe_layer_tilec;
 out vec2 oe_elev_coord;
 flat out uint64_t oe_elev_tex;
 
 void oe_contour_vertex(inout vec4 not_used)
 {
     oe_elev_coord = oe_terrain_getElevationCoord(oe_layer_tilec.st);
     oe_elev_tex = oe_terrain_getElevationHandle();
+
+#ifdef OE_USE_GEOID
+    // calculate long and lat from [0..1] across the profile:
+    vec2 r = (oe_tile_key.xy + oe_layer_tilec.xy) / exp2(oe_tile_key.z);
+    vec2 uv = vec2(0.5 * r.x, r.y);
+    oe_contour_offset = texture(oe_contour_geoid, uv).r;
+#endif
 }
 #else
-void oe_contour_vertex(inout vec4 not_used) { }
+void oe_contour_vertex(inout vec4 not_used)
+{
+#ifdef OE_USE_GEOID
+    // calculate long and lat from [0..1] across the profile:
+    vec2 r = (oe_tile_key.xy + oe_layer_tilec.xy) / exp2(oe_tile_key.z);
+    vec2 uv = vec2(0.5 * r.x, r.y);
+    oe_contour_offset = texture(oe_contour_geoid, uv).r;
+#endif
+}
 #endif
 
 [break]
 #pragma vp_function oe_contour_fragment, fragment_coloring
 #pragma import_defines(OE_USE_GL4)
+#pragma import_defines(OE_USE_GEOID)
 
 uniform sampler1D oe_contour_xfer;
 uniform float oe_contour_min;
 uniform float oe_contour_range;
 
+#ifdef OE_USE_GEOID
+in float oe_contour_offset;
+#endif
+
 #ifdef OE_USE_GL4
 
 flat in uint64_t oe_elev_tex;
@@ -36,23 +63,32 @@ void oe_contour_fragment( inout vec4 color )
 {
     if (oe_elev_tex > 0)
     {
-        float height = texture(sampler2D(oe_elev_tex), oe_elev_coord).r;
+        float offset = 0;
+#ifdef OE_USE_GEOID
+        offset = oe_contour_offset;
+#endif
+
+        float height = texture(sampler2D(oe_elev_tex), oe_elev_coord).r - offset;
         float height_normalized = (height - oe_contour_min) / oe_contour_range;
         float lookup = clamp(height_normalized, 0.0, 1.0);
         vec4 texel = texture(oe_contour_xfer, lookup);
         color.rgb = mix(color.rgb, texel.rgb, texel.a);
     }
 }
 
-#else // OE_USE_GL4
+#else // not OE_USE_GL4
 
 // GL3 implementation:
 float oe_terrain_getElevation(in vec2 uv);
 in vec4 oe_layer_tilec;
 
 void oe_contour_fragment(inout vec4 color)
 {
-    float height = oe_terrain_getElevation(oe_layer_tilec.st);
+    float offset = 0;
+#ifdef OE_USE_GEOID
+    offset = oe_contour_offset;
+#endif
+    float height = oe_terrain_getElevation(oe_layer_tilec.st) - offset;
     float height_normalized = (height - oe_contour_min) / oe_contour_range;
     float lookup = clamp(height_normalized, 0.0, 1.0);
     vec4 texel = texture(oe_contour_xfer, lookup);

src/osgEarthDrivers/engine_rex/TileNode.cpp

@@ -69,14 +69,14 @@ TileNode::TileNode(const TileKey& key, TileNode* parent, EngineContext* context,
     unsigned tw, th;
     _key.getProfile()->getNumTiles(_key.getLOD(), tw, th);
 
-    const double m = 65536; //pow(2.0, 16.0);
+    //const double m = 65536; //pow(2.0, 16.0);
 
     double x = (double)_key.getTileX();
     double y = (double)(th - _key.getTileY() - 1);
 
     _tileKeyValue.set(
-        (float)(x-tw/2), //(int)fmod(x, m),
-        (float)(y-th/2), // (int)fmod(y, m),
+        (float)x, //(float)(x-tw/2), //(int)fmod(x, m),
+        (float)y, //(float)(y-th/2), // (int)fmod(y, m),
         (float)_key.getLOD(),
         -1.0f);
 

tests/mapzen.earth

@@ -1,12 +1,9 @@
-<Map name="MapZen Elevation">
-    <options screen_space_error="75">
-        <terrain skirt_ratio="0.05" lod_method="screen_space"/>
-    </options>
-
+<Map name="MapZen Elevation">    
     <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>
+        <vdatum>egm96</vdatum>
     </XYZElevation>
 
     <ContourMap name="Color Ramp"/>