C++ Projection::UniversalLatitude方法代码示例

int main() {
  Isis::Preference::Preferences(true);
  try {
    cout << "Test 1, create a seeder" << endl;

    PvlGroup alg("PolygonSeederAlgorithm");

    if(!alg.hasKeyword("Name")) {
      cout << "Test without subgrid" << endl;
      alg += PvlKeyword("Name", "Strip");
      alg += PvlKeyword("MinimumThickness", toString(0.3));
      alg += PvlKeyword("MinimumArea", toString(10));
      alg += PvlKeyword("XSpacing", toString(1500));
      alg += PvlKeyword("YSpacing", toString(1500));
    }

    PvlObject o("AutoSeed");
    o.addGroup(alg);

    Pvl pvl;
    pvl.addObject(o);
    cout << pvl << endl << endl;

    PolygonSeeder *ps = PolygonSeederFactory::Create(pvl);

    std::cout << "Test to make sure Parse did it's job" << std::endl;
    std::cout << "MinimumThickness = " << ps->MinimumThickness() << std::endl;
    std::cout << "MinimumArea = " << ps->MinimumArea() << std::endl;

    cout << "Test 2, test a square polygon" << endl;
    try {
      // Call the seed member with a polygon
      geos::geom::CoordinateSequence *pts;
      vector<geos::geom::Geometry *> polys;

      // Create the A polygon
      pts = new geos::geom::CoordinateArraySequence();
      pts->add(geos::geom::Coordinate(0, 0));
      pts->add(geos::geom::Coordinate(0, 1.5));
      pts->add(geos::geom::Coordinate(0.5, 1.5));
      pts->add(geos::geom::Coordinate(0.5, 0));
      pts->add(geos::geom::Coordinate(0, 0));

      polys.push_back(Isis::globalFactory.createPolygon(
                        Isis::globalFactory.createLinearRing(pts), NULL));

      geos::geom::MultiPolygon *mp = Isis::globalFactory.createMultiPolygon(polys);

      cout << "Lon/Lat polygon = " << mp->toString() << endl;
      // Create the projection necessary for seeding
      PvlGroup radii = Projection::TargetRadii("MARS");
      Isis::Pvl maplab;
      maplab.addGroup(Isis::PvlGroup("Mapping"));
      Isis::PvlGroup &mapGroup = maplab.findGroup("Mapping");
      mapGroup += Isis::PvlKeyword("EquatorialRadius", (QString)radii["EquatorialRadius"]);
      mapGroup += Isis::PvlKeyword("PolarRadius", (QString)radii["PolarRadius"]);
      mapGroup += Isis::PvlKeyword("LatitudeType", "Planetocentric");
      mapGroup += Isis::PvlKeyword("LongitudeDirection", "PositiveEast");
      mapGroup += Isis::PvlKeyword("LongitudeDomain", toString(360));
      mapGroup += Isis::PvlKeyword("CenterLatitude", "0");
      mapGroup += Isis::PvlKeyword("CenterLongitude", "0");
      mapGroup += Isis::PvlKeyword("ProjectionName", "Sinusoidal");

      Projection *proj = Isis::ProjectionFactory::Create(maplab);

      /*
      This test doesn't make sense because there is no ground range on this
      projection.


      double x1,x2,y1,y2;
      proj->XYRange(x1,x2,y1,y2);
      if(fabs(x1) < 0.00000001) x1 = 0.0;
      if(fabs(x2) < 0.00000001) x2 = 0.0;
      if(fabs(y1) < 0.00000001) y1 = 0.0;
      if(fabs(y2) < 0.00000001) y2 = 0.0;
      std::cout << "X: " << x1 << "-" << x2 << " Y: " << y1 << "-" << y2 << std::endl;
      */

      geos::geom::MultiPolygon *xymp = PolygonTools::LatLonToXY(*mp, proj);
      vector<geos::geom::Point *> seedValues = ps->Seed(xymp);

      vector<geos::geom::Point *> points;
      for(unsigned int pt = 0; pt < seedValues.size(); pt ++) {
        if(proj->SetCoordinate(seedValues[pt]->getX(), seedValues[pt]->getY())) {
          points.push_back(Isis::globalFactory.createPoint(
                             geos::geom::Coordinate(proj->UniversalLongitude(),
                                 proj->UniversalLatitude())));
        }
        else {
          IString msg = "Unable to convert to a (lon,lat)";
          throw IException(IException::Programmer, msg, _FILEINFO_);
        }
      }

      cout << setprecision(13);
      for(unsigned int i = 0; i < points.size(); i++) {
        cout << "  POINT (";
        cout << points[i]->getX() << " " << points[i]->getY() << ")" << endl;
      }
//.........这里部分代码省略.........

void IsisMain() {

  // Get the list of cubes to mosaic

  UserInterface &ui = Application::GetUserInterface();
  FileList flist(ui.GetFilename("FROMLIST"));


  vector<Cube *> clist;
  try {
    if (flist.size() < 1) {
      string msg = "the list file [" +ui.GetFilename("FROMLIST") +
                   "does not contain any data";
      throw iException::Message(iException::User,msg,_FILEINFO_);
    }

    // open all the cube and place in vector clist  

    for (int i=0; i<(int)flist.size(); i++) {
      Cube *c = new Cube();
      clist.push_back(c);
      c->Open(flist[i]);
    }



    // run the compair function here.  This will conpair the 
    // labels of the first cube to the labels of each following cube. 
    PvlKeyword sourceProductId("SourceProductId");
    string ProdId;
    for (int i=0; i<(int)clist.size(); i++) {
      Pvl *pmatch = clist[0]->Label();
      Pvl *pcomp = clist[i]->Label();
      CompareLabels(*pmatch, *pcomp);
      PvlGroup g = pcomp->FindGroup("Instrument",Pvl::Traverse);
      if (g.HasKeyword("StitchedProductIds")) {
        PvlKeyword k = g["StitchedProductIds"];
        for (int j=0; j<(int)k.Size(); j++) {
          sourceProductId += g["stitchedProductIds"][j];
        }     
      }
      ProdId = (string)pmatch->FindGroup("Archive",Pvl::Traverse)["ObservationId"];
      iString bandname = (string)pmatch->FindGroup("BandBin",Pvl::Traverse)["Name"];
      bandname = bandname.UpCase();
      ProdId = ProdId + "_" + bandname;
    }
    bool runXY=true;

    //calculate the min and max lon
    double minLat = DBL_MAX;
    double maxLat = -DBL_MAX;
    double minLon = DBL_MAX;
    double maxLon = -DBL_MAX;
    double avgLat;
    double avgLon;
    for (int i=0; i<(int)clist.size(); i++) {
      Projection *proj = clist[i]->Projection();
      if (proj->MinimumLatitude() < minLat) minLat = proj->MinimumLatitude();
      if (proj->MaximumLatitude() > maxLat) maxLat = proj->MaximumLatitude();
      if (proj->MinimumLongitude() < minLon) minLon = proj->MinimumLongitude();
      if (proj->MaximumLongitude() > maxLon) maxLon = proj->MaximumLongitude();
    }
    avgLat = (minLat + maxLat) / 2;
    avgLon = (minLon + maxLon) / 2;
    Projection *proj = clist[0]->Projection();
    proj->SetGround(avgLat,avgLon);
    avgLat = proj->UniversalLatitude();
    avgLon = proj->UniversalLongitude();

    // Use camera class to get Inc., emi., phase, and other values
    double Cemiss;
    double Cphase;
    double Cincid;
    double ClocalSolTime;
    double CsolarLong;
    double CsunAzimuth;
    double CnorthAzimuth;
    for (int i=0; i<(int)clist.size(); i++) {
      Camera *cam = clist[i]->Camera();
      if (cam->SetUniversalGround(avgLat,avgLon)) {
        Cemiss = cam->EmissionAngle();
        Cphase = cam->PhaseAngle();
        Cincid = cam->IncidenceAngle();
        ClocalSolTime = cam->LocalSolarTime();
        CsolarLong = cam->SolarLongitude();
        CsunAzimuth = cam->SunAzimuth();
        CnorthAzimuth = cam->NorthAzimuth();
        runXY = false;
        break;
      }
    }

    //The code within the if runXY was added in 10/07 to find an intersect with
    //pole images that would fail when using projection set universal ground.  
    // This is run if no intersect is found when using lat and lon in 
    // projection space.
    if (runXY) {
      double startX = DBL_MAX;
      double endX = DBL_MIN;
      double startY = DBL_MAX;
//.........这里部分代码省略.........