Code to generate a low res index of the min/max elevation of the whole earth.

osgearth_computeminmaxelevation.cpp

/* -*-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.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
* IN THE SOFTWARE.
*
* 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/>
*/

#include <osg/ArgumentParser>
#include <osgEarth/Profile>
#include <osgEarth/Registry>
#include <osgEarth/TMS>

#include <iostream>
#include <fstream>
#include <map>

using namespace osgEarth;

int
usage(const std::string& msg)
{
    return -1;
}

bool computeMinMax(const TileKey& key, ElevationLayer* layer, short &min, short &max)
{
    min = SHRT_MAX;
    max = -SHRT_MAX;

    // Could  just use a gdal layer with max data level set to very high?
    GeoHeightField hf;
    TileKey keyToUse = key;

    while (keyToUse.valid())
    {
        hf = layer->createHeightField(keyToUse);
        if (hf.valid())
        {
            break;
        }
        keyToUse = keyToUse.createParentKey();
    }

    if (keyToUse != key)
    {     
        hf = hf.createSubSample(key.getExtent(), hf.getHeightField()->getNumColumns(), hf.getHeightField()->getNumRows(), osgEarth::RasterInterpolation::INTERP_BILINEAR);
    }
    

    for (unsigned int i = 0; i < hf.getHeightField()->getHeightList().size(); ++i)
    {
        float h = hf.getHeightField()->getHeightList()[i];
        if (h < min) min = h;
        if (h > max) max = h;
    }

    return true;
}

int main(int argc, char** argv)
{
    osg::ArgumentParser arguments(&argc, argv);

    const Profile* profile = osgEarth::Registry::instance()->getGlobalGeodeticProfile();

    // add a TMS elevation layer:
    TMSElevationLayer* elevation = new TMSElevationLayer();
    elevation->setURL("d:/geodata/116_tms/tms.xml");
    elevation->setVerticalDatum("egm96");
    elevation->open();

    std::ofstream fout("elevation_ranges.cpp");

    unsigned int maxLevel = 6;
    arguments.read("--maxLevel", maxLevel);


    fout << "namespace {" << std::endl;

    fout << "const unsigned int s_maxLevel=" << maxLevel << ";" << std::endl;

    for (unsigned int level = 0; level <= maxLevel; level++)
    {
        unsigned int wide, high;
        profile->getNumTiles(level, wide, high);

        unsigned int total = wide * high;

        OE_NOTICE << "Computing level " << level << " reserving " << total << std::endl;
        std::vector<short> minElevations;
        std::vector<short> maxElevations;
        minElevations.reserve(total);
        maxElevations.reserve(total);
        std::cout << "Reserved" << std::endl;

        
        unsigned int numComplete = 0;

        for (unsigned int y = 0; y < high; y++)
        {
            for (unsigned int x = 0; x < wide; x++)
            {
                TileKey key(level, x, y, profile);
                short min;
                short max;
                if (computeMinMax(key, elevation, min, max))
                {
                    minElevations.push_back(min);
                    maxElevations.push_back(max);
                    numComplete++;
                    if (numComplete % 100 == 0)
                    {
                        OE_NOTICE << "Computed " << numComplete << " of " << total << std::endl;
                    }
                }
                else
                {
                    minElevations.push_back(min);
                    maxElevations.push_back(max);
                    numComplete++;
                    //OE_NOTICE << "Failed to get elevation for " << key.str() << std::endl;
                }
            }
        }

        fout << "static short s_minElevations" << level << "[] = {" << std::endl;
        for (unsigned int i = 0; i < minElevations.size(); i++)
        {
            fout << (short)minElevations[i];
            if (i != minElevations.size() - 1) fout << ", ";
            if (i != 0 && (i + 1) % wide == 0) fout << std::endl;
        }
        fout << "};" << std::endl;

        fout << "static short s_maxElevations" << level << "[] = {" << std::endl;
        for (unsigned int i = 0; i < maxElevations.size(); i++)
        {
            fout << (short)maxElevations[i];
            if (i != maxElevations.size() - 1) fout << ", ";
            if (i != 0 && (i + 1) % wide == 0) fout << std::endl;
        }
        fout << "};" << std::endl;
    };

    fout << "static short* s_minElevationsLOD[] = {" << std::endl;
    for (unsigned int level = 0; level <= maxLevel; ++level)
    {
        fout << "s_minElevations" << level;
        if (level != maxLevel) fout << ",";
    }
    fout << "};" << std::endl;

    fout << "static short* s_maxElevationsLOD[] = {" << std::endl;
    for (unsigned int level = 0; level <= maxLevel; ++level)
    {
        fout << "s_maxElevations" << level;
        if (level != maxLevel) fout << ",";
    }
    fout << "};" << std::endl;

    fout << "}" << std::endl;
    return 0;
}