Merge pull request #160 from makaylas/cassis

ifeq ($(ISISROOT), $(BLANK))

.SILENT:

error:

	echo "Please set ISISROOT";

else

	include $(ISISROOT)/make/isismake.apps

endif

 No newline at end of file

#include "Isis.h"

#include "Application.h"

#include "Camera.h"

#include "Cube.h"

#include "FileList.h"

#include "IException.h"

#include "Pvl.h"

#include "IString.h"

#include "Longitude.h"

#include "OriginalLabel.h"

#include "Process.h"

#include "ProgramLauncher.h"

#include "TProjection.h"

#include "UserInterface.h"

using namespace std;

using namespace Isis;

#include <vector>

//functions in the code

void CompareLabels(Pvl &match, Pvl &comp);

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) {

      QString msg = "the list file [" + ui.GetFileName("FROMLIST") +

                    "does not contain any data";

      throw IException(IException::User, msg, _FILEINFO_);

    }

    // open all the cube and place in vector clist

    for(int i = 0; i < flist.size(); i++) {

      Cube *c = new Cube();

      clist.push_back(c);

      c->open(flist[i].toString());

    }

    // run the compare function here.  This will conpare the

    // labels of the first cube to the labels of each following cube.

    PvlKeyword sourceProductId("SourceProductId");

    QString 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 = (QString)pmatch->findGroup("Archive", Pvl::Traverse)["ObservationId"];

      QString bandname = (QString)pmatch->findGroup("BandBin", Pvl::Traverse)["Name"];

      bandname = bandname.toUpper();

      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++) {

      TProjection *proj = (TProjection *) 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;

    TProjection *proj = (TProjection *) 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().degrees();

        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;

      double endY =  DBL_MIN;

      for(int i = 0; i < (int)clist.size(); i++) {

        TProjection *proj = (TProjection *) clist[i]->projection();

        proj->SetWorld(0.5, 0.5);

        if(i == 0) {

          startX = proj->XCoord();

          endY = proj->YCoord();

        }

        else {

          if(proj->XCoord() < startX) startX =  proj->XCoord();

          if(proj->YCoord() > endY) endY = proj->YCoord();

        }

        Pvl *p = clist[i]->label();

        double nlines = p->findGroup("Dimensions", Pvl::Traverse)["Lines"];

        double nsamps = p->findGroup("Dimensions", Pvl::Traverse)["Samples"];

        proj->SetWorld((nsamps + 0.5), (nlines + 0.5));

        if(i == 0) {

          endX = proj->XCoord();

          startY = proj->YCoord();

        }

        else {

          if(proj->XCoord() > endX) endX =  proj->XCoord();

          if(proj->YCoord() < startY) startY = proj->YCoord();

        }

      }

      double avgX = (startX + endX) / 2;

      double avgY = (startY + endY) / 2;

      double sample = proj->ToWorldX(avgX);

      double line = proj->ToWorldY(avgY);

      for(int i = 0; i < (int)clist.size(); i++) {

        Camera *cam = clist[i]->camera();

        if(cam->SetImage(sample, line)) {

          Cemiss = cam->EmissionAngle();

          Cphase = cam->PhaseAngle();

          Cincid = cam->IncidenceAngle();

          ClocalSolTime = cam->LocalSolarTime();

          CsolarLong = cam->solarLongitude().degrees();

          CsunAzimuth = cam->SunAzimuth();

          CnorthAzimuth = cam->NorthAzimuth();

          runXY = false;

          break;

        }

      }

    }

    if(runXY) {

      QString msg = "Camera did not intersect images to gather stats";

      throw IException(IException::User, msg, _FILEINFO_);

    }

    // get the min and max SCLK values ( do this with QString comp.)

    // get the value from the original label blob

    QString startClock;

    QString stopClock;

    QString startTime;

    QString stopTime;

    for(int i = 0; i < (int)clist.size(); i++) {

      OriginalLabel origLab;

      clist[i]->read(origLab);

      PvlGroup timegrp = origLab.ReturnLabels().findGroup("TIME_PARAMETERS", Pvl::Traverse);

      if(i == 0) {

        startClock = (QString)timegrp["SpacecraftClockStartCount"];

        stopClock = (QString)timegrp["SpacecraftClockStopCount"];

        startTime = (QString)timegrp["StartTime"];

        stopTime = (QString)timegrp["StopTime"];

      }

      else {

        QString testStartTime = (QString)timegrp["StartTime"];

        QString testStopTime = (QString)timegrp["StopTime"];

        if(testStartTime < startTime) {

          startTime = testStartTime;

          startClock = (QString)timegrp["SpacecraftClockStartCount"];

        }

        if(testStopTime > stopTime) {

          stopTime = testStopTime;

          stopClock = (QString)timegrp["spacecraftClockStopCount"];

        }

      }

    }

    //  Concatenate all TDI's and summing and specialProcessingFlat into one keyword

    PvlKeyword cpmmTdiFlag("cpmmTdiFlag");

    PvlKeyword cpmmSummingFlag("cpmmSummingFlag");

    PvlKeyword specialProcessingFlag("SpecialProcessingFlag");

    for(int i = 0; i < 14; i++) {

      cpmmTdiFlag += (QString)"";

      cpmmSummingFlag += (QString)"";

      specialProcessingFlag += (QString)"";

    }

    for(int i = 0; i < (int)clist.size(); i++) {

      Pvl *clab = clist[i]->label();

      PvlGroup cInst = clab->findGroup("Instrument", Pvl::Traverse);

      OriginalLabel cOrgLab;

      clist[i]->read(cOrgLab);

      PvlGroup cGrp = cOrgLab.ReturnLabels().findGroup("INSTRUMENT_SETTING_PARAMETERS", Pvl::Traverse);

      cpmmTdiFlag[(int)cInst["CpmmNumber"]] = (QString) cGrp["MRO:TDI"];

      cpmmSummingFlag[(int)cInst["CpmmNumber"]] = (QString) cGrp["MRO:BINNING"];

      if(cInst.hasKeyword("Special_Processing_Flag")) {

        specialProcessingFlag[cInst["CpmmNumber"]] = (QString) cInst["Special_Processing_Flag"];

      }

      else {

        // there may not be the keyword Special_Processing_Flag if no

        //keyword then set the output to NOMINAL

        specialProcessingFlag[cInst["CpmmNumber"]] = "NOMINAL";

      }

    }

    // Get the blob of original labels from first image in list

    OriginalLabel org;

    clist[0]->read(org);

    //close all cubes

    for(int i = 0; i < (int)clist.size(); i++) {

      clist[i]->close();

      delete clist[i];

    }

    clist.clear();

    // automos step

    QString list = ui.GetFileName("FROMLIST");

    QString toMosaic = ui.GetFileName("TO");

    QString MosaicPriority = ui.GetString("PRIORITY");

    QString parameters = "FROMLIST=" + list + " MOSAIC=" + toMosaic + " PRIORITY=" + MosaicPriority;

    ProgramLauncher::RunIsisProgram("automos", parameters);

    // write out new information to new group mosaic

    PvlGroup mos("Mosaic");

    mos += PvlKeyword("ProductId ", ProdId);

    mos += PvlKeyword(sourceProductId);

    mos += PvlKeyword("StartTime ", startTime);

    mos += PvlKeyword("SpacecraftClockStartCount ", startClock);

    mos += PvlKeyword("StopTime ", stopTime);

    mos += PvlKeyword("SpacecraftClockStopCount ", stopClock);

    mos += PvlKeyword("IncidenceAngle ", toString(Cincid), "DEG");

    mos += PvlKeyword("EmissionAngle ", toString(Cemiss), "DEG");

    mos += PvlKeyword("PhaseAngle ", toString(Cphase), "DEG");

    mos += PvlKeyword("LocalTime ", toString(ClocalSolTime), "LOCALDAY/24");

    mos += PvlKeyword("SolarLongitude ", toString(CsolarLong), "DEG");

    mos += PvlKeyword("SubSolarAzimuth ", toString(CsunAzimuth), "DEG");

    mos += PvlKeyword("NorthAzimuth ", toString(CnorthAzimuth), "DEG");

    mos += cpmmTdiFlag;

    mos += cpmmSummingFlag;

    mos += specialProcessingFlag;

    Cube mosCube;

    mosCube.open(ui.GetFileName("TO"), "rw");

    PvlObject &lab = mosCube.label()->findObject("IsisCube");

    lab.addGroup(mos);

    //add orginal label blob to the output cube

    mosCube.write(org);

    mosCube.close();

  }

  catch(IException &e) {

    for(int i = 0; i < (int)clist.size(); i++) {

      clist[i]->close();

      delete clist[i];

    }

    QString msg = "The mosaic [" + ui.GetFileName("TO") + "] was NOT created";

    throw IException(IException::User, msg, _FILEINFO_);

  }

} // end of isis main

//Function to compare label - CompareLabels

void CompareLabels(Pvl &pmatch, Pvl &pcomp) {

  // test of the ObservationId

  PvlGroup matchgrp = pmatch.findGroup("Archive", Pvl::Traverse);

  PvlGroup compgrp = pcomp.findGroup("Archive", Pvl::Traverse);

  QString obsMatch = matchgrp["ObservationId"];

  QString obsComp = compgrp["ObservationId"];

  if(obsMatch != obsComp) {

    QString msg = "Images not from the same observation";

    throw IException(IException::User, msg, _FILEINFO_);

  }

  // Test of the BandBin filter name

  PvlGroup bmatchgrp = pmatch.findGroup("BandBin", Pvl::Traverse);

  PvlGroup bcompgrp = pcomp.findGroup("BandBin", Pvl::Traverse);

  QString bandMatch = bmatchgrp["Name"];

  QString bandComp = bcompgrp["Name"];

  if(bandMatch != bandComp) {

    QString msg = "Images not the same filter";

    throw IException(IException::User, msg, _FILEINFO_);

  }

}

<?xml version="1.0" encoding="UTF-8"?>

<application name="tgocassismos" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation="http://isis.astrogeology.usgs.gov/Schemas/Application/application.xsd">

  <brief>

    Produce an output mosaic of TGO CaSSIS images from the same observation and filter.

  </brief>

  <description>

    <p>

  If the final product wil not be an RDR, mapmos or automos can be used to produce a mosaic.  

  This program produces a mosaic and adds keywords to the label group "Mosaic" that are required 

  by hirdrgen. <br/>

  The  input cubes  to this program must be from the same observation, same filter, and map projected.<br/>  

   A new label group, required by tgocassisrdrgen, called "Mosaic"  is added to the image labels of 

   the output cube.  The original label blob is also propagated from the first file in the input list.<br/>

   </p>

<p>

  <b> processing sequence for single filter products that will become RDRs </b><br/>

         <b>cam2map</b> - project the cube.<br/>

         <b>tgocassismos</b> - mosaic files and add keywords to the mosaic group.<br/>

         <b>tgocassisrdrgen</b> - produce a RDR product.<br/>

</p>

 </description>

  <category>

    <missionItem>Trace Gas Orbiter</missionItem>

  </category>

  <history>

    <change name="Makayla Shepherd" date="2018-05-09">

      Original version.

    </change>

  </history>

  <groups>

    <group name="Files">

      <parameter name="FROMLIST">

        <type>filename</type>

        <fileMode>input</fileMode>

        <brief>

          Input list of files

        </brief>

        <description>

          A single colume list of TGO CaSSIS files that will be mosaiced together.  The input

          files must be from the some observation and be from the same filter.  The files must also

          be map projected.

        </description>

        <filter>

          *.lis *.lst *.txt

        </filter>

      </parameter>

      <parameter name="TO">

        <type>cube</type>

        <fileMode>output</fileMode>

        <brief>

          Output mosaic

        </brief>

        <description>

          The output mosaic cube with the new mosaic label group added.

        </description>

        <filter>

          *.cub

        </filter>

      </parameter>

      <parameter name="PRIORITY">

        <type>string</type>

        <default>

          <item>ONTOP</item>

        </default>

        <brief>The priority of pixel placement</brief>

        <description>

          This parameter is used to select one of two ways to masaic the

          pixels in areas of overlap.

        </description>

        <list>

          <option value="ONTOP">

            <brief>Input cube will be placed on top of the mosaic </brief>

            <description>

              When the input pixel is not Null, the input pixel

              value will be written to the output cube.

            </description>

          </option>

          <option value="BENEATH">

            <brief>Input cube will be placed beneath the mosaic</brief>

            <description>

              When the mosaic pixel is not Null it will be left unchanged.

            </description>

          </option>

        </list>

      </parameter>

    </group>

  </groups>

</application>

BLANKS = "%-6s"    

LENGTH = "%-40s"

include $(ISISROOT)/make/isismake.tststree

APPNAME = tgocassismos

include $(ISISROOT)/make/isismake.tsts

commands:

	$(LS) -1 $(INPUT)/*.cub > $(OUTPUT)/cube.lis;

	$(APPNAME) fromlist=$(OUTPUT)/cube.lis \

		   to=$(OUTPUT)/test_tgocassismos_070418.cub > /dev/null;

	catlab from=$(OUTPUT)/test_tgocassismos_070418.cub to=$(OUTPUT)/tgocassismos.pvl > /dev/null;

	$(RM) $(OUTPUT)/cube.lis;