C++ GeoPoint::y方法代码示例

float
ElevationQuery::getElevation(const GeoPoint& point,
                             double          desiredResolution,
                             double*         out_actualResolution)
{
    float result = NO_DATA_VALUE;

    sync();
    if ( point.altitudeMode() == ALTMODE_ABSOLUTE )
    {
        getElevationImpl( point, result, desiredResolution, out_actualResolution );
    }
    else
    {
        GeoPoint point_abs( point.getSRS(), point.x(), point.y(), 0.0, ALTMODE_ABSOLUTE );
        getElevationImpl( point_abs, result, desiredResolution, out_actualResolution );
    }

    return result;
}

void
MouseCoordsLabelCallback::set( const GeoPoint& mapCoords, osg::View* view, MapNode* mapNode )
{
    if ( _label.valid() )
    {
        if ( _formatter )
        {
            _label->setText( Stringify()
                <<  _formatter->format( mapCoords )
                << ", " << mapCoords.z() );
        }
        else
        {
            _label->setText( Stringify()
                << std::fixed
                << mapCoords.x()
                << ", " << mapCoords.y()
                << ", " << mapCoords.z() );
        }
    }
}

GeoPoint
RectangleNode::getCorner( Corner corner ) const
{
    GeoPoint center = getPosition();

    double earthRadius = center.getSRS()->getEllipsoid()->getRadiusEquator();
    double lat = osg::DegreesToRadians(center.y());
    double lon = osg::DegreesToRadians(center.x());
    double halfWidthMeters  =  _width.as(Units::METERS) / 2.0;
    double halfHeightMeters  = _height.as(Units::METERS) / 2.0;   

    double eastLon, eastLat;
    double westLon, westLat;
    double northLon, northLat;
    double southLon, southLat;
        
    GeoMath::destination( lat, lon, osg::DegreesToRadians( 90.0 ), halfWidthMeters, eastLat, eastLon, earthRadius );
    GeoMath::destination( lat, lon, osg::DegreesToRadians( -90.0 ), halfWidthMeters, westLat, westLon, earthRadius );
    GeoMath::destination( lat, lon, osg::DegreesToRadians( 0.0 ),  halfHeightMeters, northLat, northLon, earthRadius );
    GeoMath::destination( lat, lon, osg::DegreesToRadians( 180.0 ), halfHeightMeters, southLat, southLon, earthRadius );

    if (corner == CORNER_LOWER_LEFT)
    {
        return GeoPoint(center.getSRS(), osg::RadiansToDegrees(westLon), osg::RadiansToDegrees(southLat), 0, ALTMODE_RELATIVE);
    }
    else if (corner == CORNER_LOWER_RIGHT)
    {
        return GeoPoint(center.getSRS(), osg::RadiansToDegrees(eastLon), osg::RadiansToDegrees(southLat), 0, ALTMODE_RELATIVE);
    }
    else if (corner == CORNER_UPPER_LEFT)
    {
        return GeoPoint(center.getSRS(), osg::RadiansToDegrees(westLon), osg::RadiansToDegrees(northLat), 0, ALTMODE_RELATIVE);
    }
    else if (corner == CORNER_UPPER_RIGHT)
    {
        return GeoPoint(center.getSRS(), osg::RadiansToDegrees(eastLon), osg::RadiansToDegrees(northLat), 0, ALTMODE_RELATIVE);
    }
    return GeoPoint();
}

void LOSCreationDialog::centerMapOnNode(osg::Node* node)
{
  if (node && _map.valid() && _manager.valid() && _views)
  {
    AnnotationNode* annoNode = dynamic_cast<AnnotationNode*>(node);
    if (annoNode && annoNode->getAnnotationData() && annoNode->getAnnotationData()->getViewpoint())
    {
      _manager->doAction(this, new SetViewpointAction(osgEarth::Viewpoint(*annoNode->getAnnotationData()->getViewpoint()), *_views));
    }
    else
    {
      osg::Vec3d center = node->getBound().center();

      GeoPoint output;
      output.fromWorld( _map->getSRS(), center );
      //_map->worldPointToMapPoint(center, output);

      _manager->doAction(this, new SetViewpointAction(osgEarth::Viewpoint(
          "center", output.x(), output.y(), output.z(), 0.0, -90.0, 1e5), *_views));
    }
  }
}

bool
ElevationQuery::getElevationImpl(const GeoPoint& point,
                                 double&         out_elevation,
                                 double          desiredResolution,
                                 double*         out_actualResolution)
{
    osg::Timer_t start = osg::Timer::instance()->tick();

    if ( _mapf.elevationLayers().empty() )
    {
        // this means there are no heightfields.
        out_elevation = 0.0;
        return true;        
    }

    // tile size (resolution of elevation tiles)
    unsigned tileSize = std::max(_mapf.getMapOptions().elevationTileSize().get(), 2u);

    //This is the max resolution that we actually have data at this point
    unsigned int bestAvailLevel = getMaxLevel( point.x(), point.y(), point.getSRS(), _mapf.getProfile());

    if (desiredResolution > 0.0)
    {
        unsigned int desiredLevel = _mapf.getProfile()->getLevelOfDetailForHorizResolution( desiredResolution, tileSize );
        if (desiredLevel < bestAvailLevel) bestAvailLevel = desiredLevel;
    }

    OE_DEBUG << LC << "Best available data level " << point.x() << ", " << point.y() << " = "  << bestAvailLevel << std::endl;

    // transform the input coords to map coords:
    GeoPoint mapPoint = point;
    if ( point.isValid() && !point.getSRS()->isHorizEquivalentTo( _mapf.getProfile()->getSRS() ) )
    {
        mapPoint = point.transform(_mapf.getProfile()->getSRS());
        if ( !mapPoint.isValid() )
        {
            OE_WARN << LC << "Fail: coord transform failed" << std::endl;
            return false;
        }
    }    

    // get the tilekey corresponding to the tile we need:
    TileKey key = _mapf.getProfile()->createTileKey( mapPoint.x(), mapPoint.y(), bestAvailLevel );
    if ( !key.valid() )
    {
        OE_WARN << LC << "Fail: coords fall outside map" << std::endl;
        return false;
    }
        
    bool result = false;      
    while (!result)
    {      
        GeoHeightField geoHF;
        TileCache::Record record;
        // Try to get the hf from the cache
        if ( _cache.get( key, record ) )
        {                        
            geoHF = record.value();
        }
        else
        {
            // Create it            
            osg::ref_ptr<osg::HeightField> hf = new osg::HeightField();
            hf->allocate( tileSize, tileSize );

            // Initialize the heightfield to nodata
            for (unsigned int i = 0; i < hf->getFloatArray()->size(); i++)
            {
                hf->getFloatArray()->at( i ) = NO_DATA_VALUE;
            }   

            if (_mapf.populateHeightField( hf, key ) )
            {                
                geoHF = GeoHeightField( hf.get(), key.getExtent() );
                _cache.insert( key, geoHF );
            }
        }

        if (geoHF.valid())
        {            
            float elevation = 0.0f;                 
            result = geoHF.getElevation( mapPoint.getSRS(), mapPoint.x(), mapPoint.y(), _mapf.getMapInfo().getElevationInterpolation(), mapPoint.getSRS(), elevation);                              
            if (result && elevation != NO_DATA_VALUE)
            {                        
                // see what the actual resolution of the heightfield is.
                if ( out_actualResolution )
                    *out_actualResolution = geoHF.getXInterval(); 
                out_elevation = (double)elevation;                
                break;
            }
            else
            {                               
                result = false;
            }
        }

        if (!result)
        {
            key = key.createParentKey();                        
            if (!key.valid())
//.........这里部分代码省略.........

bool
OGR_SpatialReference::transformInPlace( GeoPoint& input ) const
{
    if ( !handle || !input.isValid() ) {
        osgGIS::notify( osg::WARN ) << "Spatial reference or input point is invalid" << std::endl;
        return false;
    }

	OGR_SpatialReference* input_sr = (OGR_SpatialReference*)input.getSRS();
    if ( !input_sr ) {
        osgGIS::notify( osg::WARN ) << "SpatialReference: input point has no SRS" << std::endl;
        return false;
    }

    // first check whether the input point is geocentric - and if so, pre-convert it to geographic:
    if ( input_sr->isGeocentric() )
    {
        input.set( input * input_sr->getInverseReferenceFrame() );
        osg::Vec3d temp = input_sr->getEllipsoid().geocentricToLatLong( input );
        input = GeoPoint( temp, input_sr->getGeographicSRS() ); 
        input_sr = static_cast<OGR_SpatialReference*>( input.getSRS() );
    }

    osg::Vec3d input_vec = input;
    bool crs_equiv = false;
    bool mat_equiv = false;
    testEquivalence( input_sr, /*out*/crs_equiv, /*out*/mat_equiv );

    // pull it out of its source frame:
    if ( !mat_equiv )
    {
        input.set( input * input_sr->inv_ref_frame );
    }    

    bool result = false;

    if ( !crs_equiv )
    {
        OGR_SCOPE_LOCK();

        //TODO: some kind of per-thread cache

        void* xform_handle = OCTNewCoordinateTransformation( input_sr->handle, this->handle );
        if ( !xform_handle ) {
            osgGIS::notify( osg::WARN ) << "Spatial Reference: SRS xform not possible" << std::endl
                << "    From => " << input_sr->getWKT() << std::endl
                << "    To   => " << this->getWKT() << std::endl;
            return false;
        }

        //TODO: figure out why xforming GEOCS x=-180 to another GEOCS doesn't work
        if ( OCTTransform( xform_handle, 1, &input.x(), &input.y(), &input.z() ) )
        {
            result = true;
        }
        else
        {
            osgGIS::notify( osg::WARN ) << "Spatial Reference: Failed to xform a point from "
                << input_sr->getName() << " to " << this->getName()
                << std::endl;
        }

        OCTDestroyCoordinateTransformation( xform_handle );
    }
    else
    {
        result = true;
    }

    // put it into the new ref frame:
    if ( !mat_equiv )
    {
        input.set( input * ref_frame );
    }

    if ( result == true )
    {
        applyTo( input );
    }

    return result;
}

void GraticuleNode::traverse(osg::NodeVisitor& nv)
{
    if (nv.getVisitorType() == osg::NodeVisitor::UPDATE_VISITOR)
    {
        updateLabels();
    }
    else if (nv.getVisitorType() == osg::NodeVisitor::CULL_VISITOR)
    {
        osgUtil::CullVisitor* cv = static_cast<osgUtil::CullVisitor*>(&nv);

        osg::Vec3d vp = cv->getViewPoint();

        osg::Matrixd viewMatrix = *cv->getModelViewMatrix();

        // Only update if the view matrix has changed.
        if (viewMatrix != _viewMatrix)
        {            
            GeoPoint eyeGeo;
            eyeGeo.fromWorld( _mapNode->getMapSRS(), vp );
            _lon = eyeGeo.x();
            _lat = eyeGeo.y();            

            osg::Viewport* viewport = cv->getViewport();

            float centerX = viewport->x() + viewport->width() / 2.0;
            float centerY = viewport->y() + viewport->height() / 2.0;

            float offsetCenterX = centerX;
            float offsetCenterY = centerY;

            bool hitValid = false;

            // Try the center of the screen.
            if(_mapNode->getTerrain()->getWorldCoordsUnderMouse(cv->getCurrentCamera()->getView(), centerX, centerY, _focalPoint))
            {
                hitValid = true;
            }

            if (hitValid)
            {
                GeoPoint focalGeo;
                focalGeo.fromWorld( _mapNode->getMapSRS(), _focalPoint );
                _lon = focalGeo.x();
                _lat = focalGeo.y();
                // We only store the previous view matrix if we actually got a hit.  Otherwise we still need to update.
                _viewMatrix = viewMatrix;
            }


            double targetResolution = (_viewExtent.height() / 180.0) / _options.gridLines().get();

            double resolution = _resolutions[0];
            for (unsigned int i = 0; i < _resolutions.size(); i++)
            {
                resolution = _resolutions[i];
                if (resolution <= targetResolution)
                {
                    break;
                }
            }

            // Trippy
            //resolution = targetResolution;

            _viewExtent = getViewExtent( cv );

            // Try to compute an approximate meters to pixel value at this view.
            double fovy, aspectRatio, zNear, zFar;
            cv->getProjectionMatrix()->getPerspective(fovy, aspectRatio, zNear, zFar);
            double dist = osg::clampAbove(eyeGeo.z(), 1.0);
            double halfWidth = osg::absolute( tan(osg::DegreesToRadians(fovy/2.0)) * dist );
            _metersPerPixel = (2.0 * halfWidth) / (double)viewport->height();

            if (_resolution != resolution)
            {
                setResolution(resolution);
            }
        }
    }
    osg::Group::traverse(nv);
}

bool
ElevationQuery::getElevationImpl(const GeoPoint& point,
                                 float&          out_elevation,
                                 double          desiredResolution,
                                 double*         out_actualResolution)
{
    // assertion.
    if ( !point.isAbsolute() )
    {
        OE_WARN << LC << "Assertion failure; input must be absolute" << std::endl;
        return false;
    }

    osg::Timer_t begin = osg::Timer::instance()->tick();

    // first try the terrain patches.
    if ( _terrainModelLayers.size() > 0 )
    {
        osgUtil::IntersectionVisitor iv;

        if ( _ivrc.valid() )
            iv.setReadCallback(_ivrc.get());

        for(LayerVector::iterator i = _terrainModelLayers.begin(); i != _terrainModelLayers.end(); ++i)
        {
            // find the scene graph for this layer:
            Layer* layer = i->get();
            osg::Node* node = layer->getNode();
            if ( node )
            {
                // configure for intersection:
                osg::Vec3d surface;
                point.toWorld( surface );

                // trivial bounds check:
                if ( node->getBound().contains(surface) )
                {
                    osg::Vec3d nvector;
                    point.createWorldUpVector(nvector);

                    osg::Vec3d start( surface + nvector*5e5 );
                    osg::Vec3d end  ( surface - nvector*5e5 );

                    // first time through, set up the intersector on demand
                    if ( !_lsi.valid() )
                    {
                        _lsi = new osgUtil::LineSegmentIntersector(start, end);
                        _lsi->setIntersectionLimit( _lsi->LIMIT_NEAREST );
                    }
                    else
                    {
                        _lsi->reset();
                        _lsi->setStart( start );
                        _lsi->setEnd  ( end );
                    }

                    // try it.
                    iv.setIntersector( _lsi.get() );
                    node->accept( iv );

                    // check for a result!!
                    if ( _lsi->containsIntersections() )
                    {
                        osg::Vec3d isect = _lsi->getIntersections().begin()->getWorldIntersectPoint();

                        // transform back to input SRS:
                        GeoPoint output;
                        output.fromWorld( point.getSRS(), isect );
                        out_elevation = (float)output.z();
                        if ( out_actualResolution )
                            *out_actualResolution = 0.0;

                        return true;
                    }
                }
                else
                {
                    //OE_INFO << LC << "Trivial rejection (bounds check)" << std::endl;
                }
            }
        }
    } 

    if (_elevationLayers.empty())
    {
        // this means there are no heightfields.
        out_elevation = NO_DATA_VALUE;
        return true;
    }

    // secure map pointer:
    osg::ref_ptr<const Map> map;
    if (!_map.lock(map))
    {
        return false;
    }    

    // tile size (resolution of elevation tiles)
    unsigned tileSize = 257; // yes?

//.........这里部分代码省略.........

    /**
     * Create and return an image for the given TileKey.
     */
    osg::Image* createImage( const TileKey& key, ProgressCallback* progress )
    {
        if (_debugDirect)
        {
            //osg::Image* image = new osg::Image;
            //image->allocateImage(256,256,1, GL_RGB, GL_UNSIGNED_BYTE);
            //return image;

            return osgDB::readImageFile( getDirectURI(key) );
            //return URI(getDirectURI(key)).getImage(_dbOptions.get(), progress);
        }

        // center point of the tile (will be in spherical mercator)
        double x, y;
        key.getExtent().getCentroid(x, y);

        // transform it to lat/long:
        GeoPoint geo;

        GeoPoint( getProfile()->getSRS(), x, y ).transform(
            getProfile()->getSRS()->getGeographicSRS(),
            geo );

        // contact the REST API. Docs are here:
        // http://msdn.microsoft.com/en-us/library/ff701716.aspx

        // construct the request URI:
        std::string request = Stringify()
            << std::setprecision(12)
            << _options.imageryMetadataAPI().get()     // base REST API
            << "/"    << _options.imagerySet().get()   // imagery set to use
            << "/"    << geo.y() << "," << geo.x()     // center point in lat/long
            << "?zl=" << key.getLOD() + 1              // zoom level
            << "&o=json"                               // response format
            << "&key=" << _options.key().get();        // API key

        // check the URI cache.
        URI                  location;
        TileURICache::Record rec;

        if ( _tileURICache.get(request, rec) )
        {
            location = URI(rec.value());
            //CacheStats stats = _tileURICache.getStats();
            //OE_INFO << "Ratio = " << (stats._hitRatio*100) << "%" << std::endl;
        }
        else
        {
            unsigned c = ++_apiCount;
            if ( c % 25 == 0 )
                OE_INFO << LC << "API calls = " << c << std::endl;
            
            // fetch it:
            ReadResult metadataResult = URI(request).readString(_dbOptions, progress);

            if ( metadataResult.failed() )
            {
                // check for a REST error:
                if ( metadataResult.code() == ReadResult::RESULT_SERVER_ERROR )
                {
                    OE_WARN << LC << "REST API request error!" << std::endl;

                    Config metadata;
                    std::string content = metadataResult.getString();
                    metadata.fromJSON( content );
                    ConfigSet errors = metadata.child("errorDetails").children();
                    for(ConfigSet::const_iterator i = errors.begin(); i != errors.end(); ++i )
                    {
                        OE_WARN << LC << "REST API: " << i->value() << std::endl;
                    }
                    return 0L;
                }
                else
                {
                    OE_WARN << LC << "Request error: " << metadataResult.getResultCodeString() << std::endl;
                }
                return 0L;
            }

            // decode it:
            Config metadata;
            if ( !metadata.fromJSON(metadataResult.getString()) )
            {
                OE_WARN << LC << "Error decoding REST API response" << std::endl;
                return 0L;
            }

            // check the vintage field. If it's empty, that means we got a "no data" tile.
            Config* vintageEnd = metadata.find("vintageEnd");
            if ( !vintageEnd || vintageEnd->value().empty() )
            {
                OE_DEBUG << LC << "NO data image encountered." << std::endl;
                return 0L;
            }

            // find the tile URI:
            Config* locationConf= metadata.find("imageUrl");
//.........这里部分代码省略.........

bool
ElevationQuery::getElevationImpl(const GeoPoint& point,
                                 double&         out_elevation,
                                 double          desiredResolution,
                                 double*         out_actualResolution)
{
    osg::Timer_t start = osg::Timer::instance()->tick();

    if ( _maxDataLevel == 0 || _tileSize == 0 )
    {
        // this means there are no heightfields.
        out_elevation = 0.0;
        return true;
    }

    
    //This is the max resolution that we actually have data at this point
    unsigned int bestAvailLevel = getMaxLevel( point.x(), point.y(), point.getSRS(), _mapf.getProfile());

    if (desiredResolution > 0.0)
    {
        unsigned int desiredLevel = _mapf.getProfile()->getLevelOfDetailForHorizResolution( desiredResolution, _tileSize );
        if (desiredLevel < bestAvailLevel) bestAvailLevel = desiredLevel;
    }

    OE_DEBUG << "Best available data level " << point.x() << ", " << point.y() << " = "  << bestAvailLevel << std::endl;

    // transform the input coords to map coords:
    GeoPoint mapPoint = point;
    if ( point.isValid() && !point.getSRS()->isEquivalentTo( _mapf.getProfile()->getSRS() ) )
    {
        mapPoint = point.transform(_mapf.getProfile()->getSRS());
        if ( !mapPoint.isValid() )
        {
            OE_WARN << LC << "Fail: coord transform failed" << std::endl;
            return false;
        }
    }

    osg::ref_ptr<osg::HeightField> tile;

    // get the tilekey corresponding to the tile we need:
    TileKey key = _mapf.getProfile()->createTileKey( mapPoint.x(), mapPoint.y(), bestAvailLevel );
    if ( !key.valid() )
    {
        OE_WARN << LC << "Fail: coords fall outside map" << std::endl;
        return false;
    }

    // Check the tile cache. Note that the TileSource already likely has a MemCache
    // attached to it. We employ a secondary cache here for a couple reasons. One, this
    // cache will store not only the heightfield, but also the tesselated tile in the event
    // that we're using GEOMETRIC mode. Second, since the call the getHeightField can 
    // fallback on a lower resolution, this cache will hold the final resolution heightfield
    // instead of trying to fetch the higher resolution one each item.

    TileCache::Record record;
    if ( _tileCache.get(key, record) )
    {
        tile = record.value().get();
    }

    // if we didn't find it, build it.
    if ( !tile.valid() )
    {
        // generate the heightfield corresponding to the tile key, automatically falling back
        // on lower resolution if necessary:
        _mapf.getHeightField( key, true, tile, 0L );

        // bail out if we could not make a heightfield a all.
        if ( !tile.valid() )
        {
            OE_WARN << LC << "Unable to create heightfield for key " << key.str() << std::endl;
            return false;
        }

        _tileCache.insert(key, tile.get());
    }

    OE_DEBUG << LC << "LRU Cache, hit ratio = " << _tileCache.getStats()._hitRatio << std::endl;

    // see what the actual resolution of the heightfield is.
    if ( out_actualResolution )
        *out_actualResolution = (double)tile->getXInterval();

    bool result = true;

    const GeoExtent& extent = key.getExtent();
    double xInterval = extent.width()  / (double)(tile->getNumColumns()-1);
    double yInterval = extent.height() / (double)(tile->getNumRows()-1);
    
    out_elevation = (double) HeightFieldUtils::getHeightAtLocation( 
        tile.get(), 
        mapPoint.x(), mapPoint.y(), 
        extent.xMin(), extent.yMin(), 
        xInterval, yInterval, _mapf.getMapInfo().getElevationInterpolation() );

    osg::Timer_t end = osg::Timer::instance()->tick();
    _queries++;
    _totalTime += osg::Timer::instance()->delta_s( start, end );
//.........这里部分代码省略.........

bool
MGRSFormatter::transform( const GeoPoint& input, MGRSCoord& out ) const
{
    if ( !input.isValid() )
        return false;

    // convert to lat/long if necessary:
    GeoPoint inputGeo = input;
    if ( !inputGeo.makeGeographic() )
        return false;

    unsigned    zone;
    char        gzd;
    unsigned    x=0, y=0;
    char        sqid[3];
    std::string space;

    if ( _options & USE_SPACES )
        space = " ";

    sqid[0] = '?';
    sqid[1] = '?';
    sqid[2] = 0;

    double latDeg = inputGeo.y();
    double lonDeg = inputGeo.x();

    if ( latDeg >= 84.0 || latDeg <= -80.0 ) // polar projection
    {
        bool isNorth = latDeg > 0.0;
        zone = 0;
        gzd = isNorth ? (lonDeg < 0.0 ? 'Y' : 'Z') : (lonDeg < 0.0? 'A' : 'B');

        osg::ref_ptr<const SpatialReference> ups =
            const_cast<MGRSFormatter*>(this)->_srsCache[ s_polarZoneSpecs[isNorth?0:1] ];
        if (!ups.valid())
            ups = SpatialReference::create( s_polarZoneSpecs[isNorth?0:1] );

        if ( !ups.valid() )
        {
            OE_WARN << LC << "Failed to create UPS SRS" << std::endl;
            return false;
        }

        osg::Vec3d upsCoord;
        if ( _refSRS->transform(osg::Vec3d(lonDeg,latDeg,0), ups.get(), upsCoord) == false )
        {
            OE_WARN << LC << "Failed to transform lat/long to UPS" << std::endl;
            return false;
        }

        int sqXOffset = upsCoord.x() >= 0.0 ? (int)floor(upsCoord.x()/100000.0) : -(int)floor(1.0-(upsCoord.x()/100000.0));
        int sqYOffset = upsCoord.y() >= 0.0 ? (int)floor(upsCoord.y()/100000.0) : -(int)floor(1.0-(upsCoord.y()/100000.0));

        int alphaOffset = isNorth ? 7 : 12;

        sqid[0] = UPS_COL_ALPHABET[ (UPS_COL_ALPHABET_SIZE+sqXOffset) % UPS_COL_ALPHABET_SIZE ];
        sqid[1] = UPS_ROW_ALPHABET[alphaOffset + sqYOffset];

        x = (unsigned)(upsCoord.x() - (100000.0*(double)sqXOffset));
        y = (unsigned)(upsCoord.y() - (100000.0*(double)sqYOffset));
    }

    else // UTM
    {
        // figure out the grid zone designator
        unsigned gzdIndex = ((unsigned)(latDeg+80.0))/8;
        gzd = GZD_ALPHABET[gzdIndex];

        // figure out the UTM zone:
        zone = (unsigned)floor((lonDeg+180.0)/6.0);   // [0..59]
        bool north = latDeg >= 0.0;

        // convert the input coordinates to UTM:
        // yes, always use +north so we get Y relative to equator

        // using an SRS cache speed things up a lot..
        osg::ref_ptr<const SpatialReference>& utm =
            const_cast<MGRSFormatter*>(this)->_srsCache[s_lateralZoneSpecs[zone]];
        if ( !utm.valid() )
            utm = SpatialReference::create( s_lateralZoneSpecs[zone] );

        osg::Vec3d utmCoord;
        if ( _refSRS->transform( osg::Vec3d(lonDeg,latDeg,0), utm.get(), utmCoord) == false )
        {
            OE_WARN << LC << "Error transforming lat/long into UTM" << std::endl;
            return false;
        }

        // the alphabet set:
        unsigned set = zone % 6; // [0..5]

        // find the horizontal SQID offset (100KM increments) from the central meridian:
        unsigned xSetOffset = 8 * (set % 3);
        double xMeridianOffset = utmCoord.x() - 500000.0;
        int sqMeridianOffset = xMeridianOffset >= 0.0 ? (int)floor(xMeridianOffset/100000.0) : -(int)floor(1.0-(xMeridianOffset/100000.0));
        unsigned indexOffset = (4 + sqMeridianOffset);
        sqid[0] = UTM_COL_ALPHABET[xSetOffset + indexOffset];
        double xWest = 500000.0 + (100000.0*(double)sqMeridianOffset);
        x = (unsigned)(utmCoord.x() - xWest);
//.........这里部分代码省略.........

void
ModelSplatter::operator()(const TileKey& key, osg::Node* node)
{
    TerrainTileNode* tile = osgEarth::findTopMostNodeOfType<TerrainTileNode>(node);
    if ( !tile )
        return;

    if ( key.getLOD() >= _minLOD && _model.valid() )
    {
        // make sure the correct model is loaded
        establish();

        // elevation texture and matrix are required
        osg::Texture* elevationTex = tile->getElevationTexture();
        if ( !elevationTex )
        {
            //OE_WARN << LC << "No elevation texture for key " << key.str() << "\n";
            return;
        }

        osg::RefMatrix* elevationTexMat = tile->getElevationTextureMatrix();
        if ( !elevationTexMat )
        {
            //OE_WARN << LC << "No elevation texture matrix for key " << key.str() << "\n";
            return;
        }
        
        tile->addChild( _model.get() );

        osg::StateSet* ss = tile->getOrCreateStateSet();

        // first, a rotation vector to make trees point up.
        GeoPoint p;
        key.getExtent().getCentroid(p);
        osg::Vec3d up;
        p.createWorldUpVector(up);
        osg::Quat q;
        q.makeRotate(osg::Vec3d(0,0,1), up);
        osg::Matrixd zup = osg::Matrixd::rotate(q);

        // matrices to resolve the weird terrain localization into a usable LTP.
        osg::Matrix tile2world = tile->getMatrix();
        osg::Matrix world2ltp;
        p.createWorldToLocal(world2ltp);
        osg::Matrix local2ltp = tile2world * world2ltp;
        osg::Matrix ltp2local;
        ltp2local.invert(local2ltp);

        // after inverting the matrix, combine the ZUP (optimization)
        local2ltp.preMult( zup );

        ss->addUniform( new osg::Uniform("oe_trees_local2ltp", osg::Matrixf(local2ltp)) );
        ss->addUniform( new osg::Uniform("oe_trees_ltp2local", osg::Matrixf(ltp2local)) );

        // calculate the scatter area:
        float h = key.getExtent().height() * 111320.0f;
        float w = key.getExtent().width() * 111320.0f * cos(fabs(osg::DegreesToRadians(p.y())));
        ss->addUniform( new osg::Uniform("oe_trees_span", osg::Vec2f(w,h)) );
        
        ss->setTextureAttributeAndModes(2, tile->getElevationTexture(), 1);
        ss->addUniform( new osg::Uniform("oe_terrain_tex_matrix", osg::Matrixf(*elevationTexMat)) );        
    }
}

    void update( float x, float y, osgViewer::View* view )
    {
        bool yes = false;

        // look under the mouse:
        osg::Vec3d world;
        osgUtil::LineSegmentIntersector::Intersections hits;
        if ( view->computeIntersections(x, y, hits) )
        {
            world = hits.begin()->getWorldIntersectPoint();

            // convert to map coords:
            GeoPoint mapPoint;
            mapPoint.fromWorld( s_mapNode->getMapSRS(), world );

            // Depending on the level of detail key you request, you will get a mesh that should line up exactly with the highest resolution mesh that the terrain engine will draw.
            // At level 15 that is a 257x257 heightfield.  If you select a higher lod, the mesh will be less dense.
            TileKey key = s_mapNode->getMap()->getProfile()->createTileKey(mapPoint.x(), mapPoint.y(), 15);
            OE_NOTICE << "Creating tile " << key.str() << std::endl;
            osg::ref_ptr<osg::Node> node = s_mapNode->getTerrainEngine()->createTile(key);
            if (node.valid())
            {
                // Extract the triangles from the node that was created and do our own rendering.  Simulates what you would do when passing in the triangles to a physics engine.
                OE_NOTICE << "Created tile for " << key.str() << std::endl;
                CollectTrianglesVisitor v;
                node->accept(v);
                node = v.buildNode();

                if (_node.valid())
                {
                    s_root->removeChild( _node.get() );
                }

                osg::Group* group = new osg::Group;

                // Show the actual mesh.
                group->addChild( node.get() );

                _node = group;

                // Clamp the marker to the intersection of the triangles created by osgEarth.  This should line up with the mesh that is actually rendered.
                double z = 0.0;
                s_mapNode->getTerrain()->getHeight( node, s_mapNode->getMapSRS(), mapPoint.x(), mapPoint.y(), &z);

                GeoTransform* xform = new GeoTransform();
                xform->setPosition( osgEarth::GeoPoint(s_mapNode->getMapSRS(),mapPoint.x(),  mapPoint.y(), z, ALTMODE_ABSOLUTE) );
                xform->addChild( marker.get() );
                group->addChild( xform );

                s_root->addChild( _node.get() );
            }
            else
            {
                OE_NOTICE << "Failed to create tile for " << key.str() << std::endl;
            }
        }
    }

// Calculates a sub-extent of a larger extent, given the number of children and
// the child number. This currently assumes the subdivision ordering used by
// VirtualPlanetBuilder.
GeoExtent
TerrainUtils::getSubExtent(const GeoExtent& extent,
                           int              num_children,
                           int              child_no)
{
    GeoPoint centroid = extent.getCentroid();
    GeoExtent sub_extent;

    switch( num_children )
    {
    case 0:
    case 1:
        sub_extent = extent;
        break;

    case 2:
        if ( child_no == 0 )
        {
            sub_extent = GeoExtent(
                             extent.getXMin(),
                             extent.getYMin(),
                             centroid.x(),
                             extent.getYMax(),
                             extent.getSRS() );
        }
        else
        {
            sub_extent = GeoExtent(
                             centroid.x(),
                             extent.getYMin(),
                             extent.getXMax(),
                             extent.getYMax(),
                             extent.getSRS() );
        }
        break;

    case 4:
        if ( child_no == 2 )
        {
            sub_extent = GeoExtent(
                             extent.getXMin(),
                             centroid.y(),
                             centroid.x(),
                             extent.getYMax(),
                             extent.getSRS() );
        }
        else if ( child_no == 3 )
        {
            sub_extent = GeoExtent(
                             centroid.x(),
                             centroid.y(),
                             extent.getXMax(),
                             extent.getYMax(),
                             extent.getSRS() );
        }
        else if ( child_no == 0 )
        {
            sub_extent = GeoExtent(
                             extent.getXMin(),
                             extent.getYMin(),
                             centroid.x(),
                             centroid.y(),
                             extent.getSRS() );
        }
        else if ( child_no == 1 )
        {
            sub_extent = GeoExtent(
                             centroid.x(),
                             extent.getYMin(),
                             extent.getXMax(),
                             centroid.y(),
                             extent.getSRS() );
        }
    }

    return sub_extent;
}

bool
getIsectPoint( const GeoPoint& p1, const GeoPoint& p2, const GeoPoint& p3, const GeoPoint& p4, GeoPoint& out )
{
    double denom = (p4.y()-p3.y())*(p2.x()-p1.x()) - (p4.x()-p3.x())*(p2.y()-p1.y());
    if ( denom == 0.0 )
    {
        out = GeoPoint::invalid(); // parallel lines
        return false;
    }

    double ua_num = (p4.x()-p3.x())*(p1.y()-p3.y()) - (p4.y()-p3.y())*(p1.x()-p3.x());
    double ub_num = (p2.x()-p1.x())*(p1.y()-p3.y()) - (p2.y()-p1.y())*(p1.x()-p3.x());

    double ua = ua_num/denom;
    double ub = ub_num/denom;

    if ( ua < 0.0 || ua > 1.0 ) // || ub < 0.0 || ub > 1.0 )
    {
        out = GeoPoint::invalid(); // isect point is on line, but not on line segment
        return false;
    }

    double x = p1.x() + ua*(p2.x()-p1.x());
    double y = p1.y() + ua*(p2.y()-p1.y());
    double z = p1.getDim() > 2? p1.z() + ua*(p2.z()-p1.z()) : 0.0; //right?
    out = GeoPoint( x, y, z, p1.getSRS() );
    return true;
}