Unverified Commit 538a6d37 authored by AustinSanders's avatar AustinSanders Committed by GitHub
Browse files

crism2isis conversion, notebooks, and gtests (#4169)

* crism2isis conversion,notebooks, and gtests

* Fixed error message typo

* Adds 4 new tests

* Added missing data

* More missing test data
parent 5622f85e
Loading
Loading
Loading
Loading
+326 −0
Original line number Diff line number Diff line
%% Cell type:code id: tags:

``` python
import pvl
import struct
import matplotlib.pyplot as plt
import numpy as np
import datetime
import os.path
import binascii
```

%% Cell type:code id: tags:

``` python
crism_file = '/home/arsanders/testData/crism/tsts/mrdr/input/t1865_mrrde_70n185_0256_1.lbl'
image_file = crism_file
```

%% Cell type:code id: tags:

``` python
header = pvl.load(crism_file)
```

%% Cell type:code id: tags:

``` python
header
```

%% Output

    PVLModule([
      ('PDS_VERSION_ID', 'PDS3')
      ('RECORD_TYPE', 'FIXED_LENGTH')
      ('RECORD_BYTES', 3920)
      ('FILE_RECORDS', 30720)
      ('LABEL_REVISION_NOTE',
       '2003-11-19, S. Slavney (GEO); 2005-09-25, S. Murchie (JHU/APL); 2007-03-09, '
       'E. Malaret (ACT Corp.); 2007-09-14, C. Hash (ACT Corp)')
      ('^IMAGE', 'T1865_MRRDE_70N185_0256_1.IMG')
      ('DATA_SET_ID', 'MRO-M-CRISM-5-RDR-MULTISPECTRAL-V1.0')
      ('PRODUCT_ID', 'T1865_MRRDE_70N185_0256_1')
      ('INSTRUMENT_HOST_NAME', 'MARS RECONNAISSANCE ORBITER')
      ('SPACECRAFT_ID', 'MRO')
      ('INSTRUMENT_NAME', 'COMPACT RECONNAISSANCE IMAGING SPECTROMETER FOR MARS')
      ('INSTRUMENT_ID', 'CRISM')
      ('TARGET_NAME', 'MARS')
      ('PRODUCT_TYPE', 'MAP_PROJECTED_MULTISPECTRAL_RDR')
      ('PRODUCT_CREATION_TIME',
       datetime.datetime(2007, 12, 22, 16, 50, 47, 432000, tzinfo=datetime.timezone.utc))
      ('START_TIME', 'N/A')
      ('STOP_TIME', 'N/A')
      ('SPACECRAFT_CLOCK_START_COUNT', 'N/A')
      ('SPACECRAFT_CLOCK_STOP_COUNT', 'N/A')
      ('PRODUCT_VERSION_ID', '1')
      ('PRODUCER_INSTITUTION_NAME', 'APPLIED PHYSICS LABORATORY')
      ('SOFTWARE_NAME', 'PIPE_create_crism_mrdr_product')
      ('SOFTWARE_VERSION_ID', '12.21.07')
      ('MRO:WAVELENGTH_FILE_NAME', 'T1865_MRRWV_70N185_0256_1.TAB')
      ('IMAGE',
       {'BANDS': 24,
        'BAND_NAME': ['Solar longitude, deg',
                      'Solar distance, AU',
                      'Observation ID, VNIR',
                      'Observation ID, IR',
                      'Ordinal counter, VNIR',
                      'Ordinal counter, IR',
                      'Column number, VNIR',
                      'Column number, IR',
                      'Line, or frame number, VNIR',
                      'Line, or frame number, IR',
                      'INA at areoid, deg',
                      'EMA at areoid, deg',
                      'Phase angle, deg',
                      'Latitude, areocentric, deg N',
                      'Longitude, areocentric, deg E',
                      'INA at surface from MOLA, deg',
                      'EMA at surface from MOLA, deg',
                      'Slope magnitude from MOLA, deg',
                      'MOLA slope azimuth, deg clkwise from N',
                      'Elevation, meters relative to MOLA',
                      'Thermal inertia, J m^-2 K^-1 s^-0.5',
                      'Bolometic albedo',
                      'Local solar time, hours',
                      'Spare'],
        'BAND_STORAGE_TYPE': 'BAND_SEQUENTIAL',
        'LINES': 1280,
        'LINE_SAMPLES': 980,
        'SAMPLE_BITS': 32,
        'SAMPLE_TYPE': 'PC_REAL'})
      ('IMAGE_MAP_PROJECTION',
       {'A_AXIS_RADIUS': Quantity(value=3396, units='KM'),
        'B_AXIS_RADIUS': Quantity(value=3396, units='KM'),
        'CENTER_LATITUDE': Quantity(value=67.5000001, units='DEGREE'),
        'CENTER_LONGITUDE': Quantity(value=185.0, units='DEGREE'),
        'COORDINATE_SYSTEM_NAME': 'PLANETOCENTRIC',
        'COORDINATE_SYSTEM_TYPE': 'BODY-FIXED ROTATING',
        'C_AXIS_RADIUS': Quantity(value=3396, units='KM'),
        'EASTERNMOST_LONGITUDE': Quantity(value=9.9999999, units='DEGREE'),
        'FIRST_STANDARD_PARALLEL': 'N/A',
        'LINE_FIRST_PIXEL': 1,
        'LINE_LAST_PIXEL': 1280,
        'LINE_PROJECTION_OFFSET': 18560.0,
        'MAP_PROJECTION_ROTATION': 0.0,
        'MAP_PROJECTION_TYPE': 'EQUIRECTANGULAR',
        'MAP_RESOLUTION': Quantity(value=256, units='PIXEL/DEGREE'),
        'MAP_SCALE': Quantity(value=0.231528833585, units='KM/PIXEL'),
        'MAXIMUM_LATITUDE': Quantity(value=72.5, units='DEGREE'),
        'MINIMUM_LATITUDE': Quantity(value=67.5000001, units='DEGREE'),
        'POSITIVE_LONGITUDE_DIRECTION': 'EAST',
        'REFERENCE_LATITUDE': 'N/A',
        'REFERENCE_LONGITUDE': 'N/A',
        'SAMPLE_FIRST_PIXEL': 1,
        'SAMPLE_LAST_PIXEL': 980,
        'SAMPLE_PROJECTION_OFFSET': 47360.0,
        'SECOND_STANDARD_PARALLEL': 'N/A',
        'WESTERNMOST_LONGITUDE': Quantity(value=360.0, units='DEGREE'),
        '^DATA_SET_MAP_PROJECTION': 'MRR_MAP.CAT'})
    ])

%% Cell type:code id: tags:

``` python
with open(crism_file, 'rb') as f:
    # If detached label, "^IMAGE" will be a list.
    image_file = os.path.dirname(crism_file) + "/" + header["^IMAGE"].lower()
    with open(image_file, 'rb') as im_f:
        b_image_data = im_f.read()
```

%% Cell type:code id: tags:

``` python
n_lines = 10
line_length = header['IMAGE']['LINE_SAMPLES'] * (header['IMAGE']['SAMPLE_BITS']//8)
```

%% Cell type:code id: tags:

``` python
image_data = []
for j in range(n_lines):
    image_sample = np.frombuffer(b_image_data[j*line_length:(j+1)*line_length], dtype=np.float32, count=int(line_length/4))
    image_data.append(image_sample)
image_data = np.array(image_data)
```

%% Cell type:code id: tags:

``` python
plt.imshow(image_data, cmap='magma')
```

%% Output

    <matplotlib.image.AxesImage at 0x2aff238d0850>


%% Cell type:code id: tags:

``` python
image_fn, image_ext = os.path.splitext(image_file)
mini_image_fn = image_fn + '_cropped' + image_ext
mini_image_bn = os.path.basename(mini_image_fn)

# Overwrite the number of lines in the label
header['IMAGE']['LINES'] = n_lines
header['IMAGE']['BANDS'] = 1
header['^IMAGE'] = mini_image_bn
header['FILE_RECORDS'] = n_lines -1
```

%% Cell type:code id: tags:

``` python
class RealIsisCubeLabelEncoder(pvl.encoder.ISISEncoder):
    def encode_time(self, value):
        if value.microsecond:
            second = u'%02d.%06d' % (value.second, value.microsecond)
        else:
            second = u'%02d' % value.second

        time = u'%02d:%02d:%s' % (value.hour, value.minute, second)
        return time
```

%% Cell type:code id: tags:

``` python
label_fn, label_ext = os.path.splitext(crism_file)
out_label = label_fn + '_cropped' + label_ext

grammar = pvl.grammar.ISISGrammar()
grammar.comments+=(("#", "\n"), )
encoder = RealIsisCubeLabelEncoder()
pvl.dump(header, out_label, encoder=encoder, grammar=grammar)
```

%% Output

    3539

%% Cell type:code id: tags:

``` python
with open(mini_image_fn, 'wb+') as f:
    b_reduced_image_data = image_data.tobytes()
    f.seek(0, 2)
    f.write(b_reduced_image_data)
```

%% Cell type:code id: tags:

``` python
```

%% Cell type:code id: tags:

``` python
header = pvl.load(out_label)
```

%% Cell type:code id: tags:

``` python
header
```

%% Output

    PVLModule([
      ('PDS_VERSION_ID', 'PDS3')
      ('RECORD_TYPE', 'FIXED_LENGTH')
      ('RECORD_BYTES', 3920)
      ('FILE_RECORDS', 9)
      ('LABEL_REVISION_NOTE',
       '2003-11-19, S. Slavney (GEO); 2005-09-25, S. Murchie (JHU/APL); 2007-03-09, '
       'E. Malaret (ACT Corp.); 2007-09-14, C. Hash (ACT Corp)')
      ('^IMAGE', 't1865_mrrde_70n185_0256_1_cropped.img')
      ('DATA_SET_ID', 'MRO-M-CRISM-5-RDR-MULTISPECTRAL-V1.0')
      ('PRODUCT_ID', 'T1865_MRRDE_70N185_0256_1')
      ('INSTRUMENT_HOST_NAME', 'MARS RECONNAISSANCE ORBITER')
      ('SPACECRAFT_ID', 'MRO')
      ('INSTRUMENT_NAME', 'COMPACT RECONNAISSANCE IMAGING SPECTROMETER FOR MARS')
      ('INSTRUMENT_ID', 'CRISM')
      ('TARGET_NAME', 'MARS')
      ('PRODUCT_TYPE', 'MAP_PROJECTED_MULTISPECTRAL_RDR')
      ('PRODUCT_CREATION_TIME',
       datetime.datetime(2007, 12, 22, 16, 50, 47, 432000, tzinfo=datetime.timezone.utc))
      ('START_TIME', 'N/A')
      ('STOP_TIME', 'N/A')
      ('SPACECRAFT_CLOCK_START_COUNT', 'N/A')
      ('SPACECRAFT_CLOCK_STOP_COUNT', 'N/A')
      ('PRODUCT_VERSION_ID', '1')
      ('PRODUCER_INSTITUTION_NAME', 'APPLIED PHYSICS LABORATORY')
      ('SOFTWARE_NAME', 'PIPE_create_crism_mrdr_product')
      ('SOFTWARE_VERSION_ID', '12.21.07')
      ('MRO:WAVELENGTH_FILE_NAME', 'T1865_MRRWV_70N185_0256_1.TAB')
      ('IMAGE',
       {'BANDS': 1,
        'BAND_NAME': ['Solar longitude, deg',
                      'Solar distance, AU',
                      'Observation ID, VNIR',
                      'Observation ID, IR',
                      'Ordinal counter, VNIR',
                      'Ordinal counter, IR',
                      'Column number, VNIR',
                      'Column number, IR',
                      'Line, or frame number, VNIR',
                      'Line, or frame number, IR',
                      'INA at areoid, deg',
                      'EMA at areoid, deg',
                      'Phase angle, deg',
                      'Latitude, areocentric, deg N',
                      'Longitude, areocentric, deg E',
                      'INA at surface from MOLA, deg',
                      'EMA at surface from MOLA, deg',
                      'Slope magnitude from MOLA, deg',
                      'MOLA slope azimuth, deg clkwise from N',
                      'Elevation, meters relative to MOLA',
                      'Thermal inertia, J m^-2 K^-1 s^-0.5',
                      'Bolometic albedo',
                      'Local solar time, hours',
                      'Spare'],
        'BAND_STORAGE_TYPE': 'BAND_SEQUENTIAL',
        'LINES': 10,
        'LINE_SAMPLES': 980,
        'SAMPLE_BITS': 32,
        'SAMPLE_TYPE': 'PC_REAL'})
      ('IMAGE_MAP_PROJECTION',
       {'A_AXIS_RADIUS': Quantity(value=3396, units='KM'),
        'B_AXIS_RADIUS': Quantity(value=3396, units='KM'),
        'CENTER_LATITUDE': Quantity(value=67.5000001, units='DEGREE'),
        'CENTER_LONGITUDE': Quantity(value=185.0, units='DEGREE'),
        'COORDINATE_SYSTEM_NAME': 'PLANETOCENTRIC',
        'COORDINATE_SYSTEM_TYPE': 'BODY-FIXED ROTATING',
        'C_AXIS_RADIUS': Quantity(value=3396, units='KM'),
        'EASTERNMOST_LONGITUDE': Quantity(value=9.9999999, units='DEGREE'),
        'FIRST_STANDARD_PARALLEL': 'N/A',
        'LINE_FIRST_PIXEL': 1,
        'LINE_LAST_PIXEL': 1280,
        'LINE_PROJECTION_OFFSET': 18560.0,
        'MAP_PROJECTION_ROTATION': 0.0,
        'MAP_PROJECTION_TYPE': 'EQUIRECTANGULAR',
        'MAP_RESOLUTION': Quantity(value=256, units='PIXEL/DEGREE'),
        'MAP_SCALE': Quantity(value=0.231528833585, units='KM/PIXEL'),
        'MAXIMUM_LATITUDE': Quantity(value=72.5, units='DEGREE'),
        'MINIMUM_LATITUDE': Quantity(value=67.5000001, units='DEGREE'),
        'POSITIVE_LONGITUDE_DIRECTION': 'EAST',
        'REFERENCE_LATITUDE': 'N/A',
        'REFERENCE_LONGITUDE': 'N/A',
        'SAMPLE_FIRST_PIXEL': 1,
        'SAMPLE_LAST_PIXEL': 980,
        'SAMPLE_PROJECTION_OFFSET': 47360.0,
        'SECOND_STANDARD_PARALLEL': 'N/A',
        'WESTERNMOST_LONGITUDE': Quantity(value=360.0, units='DEGREE'),
        '^DATA_SET_MAP_PROJECTION': 'MRR_MAP.CAT'})
    ])

%% Cell type:code id: tags:

``` python
```
+647 −0

File added.

Preview size limit exceeded, changes collapsed.

+173 −0
Original line number Diff line number Diff line
#include "crism2isis.h"
#include "ProcessImportPds.h"
#include "UserInterface.h"
#include "Pvl.h"
#include "FileName.h"
#include "TextFile.h"

using namespace std;

namespace Isis{
  void crism2isis(UserInterface &ui, Pvl *log) {
    ProcessImportPds p;
    Pvl pdsLabel;
    PvlGroup results;

    FileName inFile = ui.GetFileName("FROM");

    p.SetPdsFile(inFile.expanded(), "", pdsLabel);
    // 65535 is set to NULL
    p.SetNull(65535, 65535);

    CubeAttributeOutput &att = ui.GetOutputAttribute("TO");
    Cube *ocube = p.SetOutputCube(ui.GetFileName("TO"), att);

    Pvl outLabel;

    Pvl labelPvl(inFile.expanded());

    QString prodType;

    if (labelPvl.hasKeyword("PRODUCT_TYPE")) {
      prodType = (QString)labelPvl.findKeyword("PRODUCT_TYPE");
    }
    else {
      QString msg = "Unsupported CRISM file type, supported types are: DDR, MRDR, and TRDR";
      throw IException(IException::User, msg, _FILEINFO_);
    }

    if (prodType.toUpper() == "MAP_PROJECTED_MULTISPECTRAL_RDR") {
      QString prodId;
      if (labelPvl.hasKeyword("PRODUCT_ID")) {
        prodId = (QString)labelPvl.findKeyword("PRODUCT_ID");
        prodId = prodId.mid(prodId.indexOf("_") + 1, prodId.indexOf("_"));
      }
      else {
        QString msg = "Could not find label PRODUCT_ID, invalid MRDR";
        throw IException(IException::Unknown, msg, _FILEINFO_);
      }

      //If the product type is AL (Lambert albedo) or IF (I/F)
      //Read the wavelength table and put the corresponding
      //widths in the band bin group
      if (prodId.toUpper() == "MRRAL" || prodId.toUpper() == "MRRIF") {
        //If the wavelength file is specified in the label
        if (labelPvl.hasKeyword("MRO:WAVELENGTH_FILE_NAME")) {
          PvlGroup bandBin = PvlGroup("BandBin");
          PvlKeyword origBand = PvlKeyword("OriginalBand");
          PvlKeyword widths = PvlKeyword("Width");
          QString tablePath = (QString)labelPvl.findKeyword("MRO:WAVELENGTH_FILE_NAME");
          tablePath = tablePath.toLower();
          FileName tableFile(inFile.path() + "/" + tablePath);
          //Check if the wavelength file exists
          if (tableFile.fileExists()) {
            TextFile *fin = new TextFile(tableFile.expanded());
            // Open table file
            if (!fin->OpenChk()) {
              QString msg = "Cannot open wavelength table [" + tableFile.expanded() + "]";
              throw IException(IException::Io, msg, _FILEINFO_);
            }

            //For each line in the wavelength table, add the width to
            //The band bin group
            QString st;
            int band = 1;
            while (fin->GetLine(st)) {
              st = st.simplified().trimmed();
              QStringList cols = st.split(",");

              origBand += toString(band);
              widths += cols[2];
              band++;
            }
            delete fin;

            bandBin.addKeyword(origBand);
            bandBin.addKeyword(widths);
            ocube->putGroup(bandBin);
          }
          //Otherwise throw an error
          else {
            QString msg = "Cannot find wavelength table [" + tableFile.expanded() + "]";
            throw IException(IException::Io, msg, _FILEINFO_);
          }
        }
      }
      //If the product type is DE (Derived products for I/F) or DL
      //(Derived products for Lambert albedo) write the band names
      //to the band bin group
      else if (prodId.toUpper() == "MRRDE" || prodId.toUpper() == "MRRDL") {
        PvlGroup bandBin = PvlGroup("BandBin");
        PvlKeyword origBand = PvlKeyword("OriginalBand");
        PvlKeyword bandName = PvlKeyword("BandName");
        PvlKeyword bandNames = labelPvl.findObject("IMAGE").findKeyword("BAND_NAME");
        for (int i = 0; i < bandNames.size(); i++) {
          origBand += toString(i + 1);
          bandName += bandNames[i];
        }
        bandBin.addKeyword(origBand);
        bandBin.addKeyword(bandName);
        ocube->putGroup(bandBin);
      }
      //Translate the Mapping group
      p.TranslatePdsProjection(outLabel);
      ocube->putGroup(outLabel.findGroup("Mapping", Pvl::Traverse));

      //  Check for any change from the default projection offsets and multipliers to log

      if (p.GetProjectionOffsetChange()) {
        results = p.GetProjectionOffsetGroup();
        results[0].addComment("Projection offsets and multipliers have been changed from");
        results[0].addComment("defaults. New values are below.");
      }

    }
    else if (prodType.toUpper() == "TARGETED_RDR") {
    }
    else if (prodType.toUpper() == "DDR") {
      PvlGroup bandBin = PvlGroup("BandBin");
      PvlKeyword origBand = PvlKeyword("OriginalBand");
      PvlKeyword bandName = PvlKeyword("BandName");
      PvlKeyword bandNames = labelPvl.findObject("FILE").findObject("IMAGE").findKeyword("BAND_NAME");
      for (int i = 0; i < bandNames.size(); i++) {
        origBand += toString(i + 1);
        bandName += bandNames[i];
      }
      bandBin.addKeyword(origBand);
      bandBin.addKeyword(bandName);
      ocube->putGroup(bandBin);
    }
    else {
      QString msg = "Unsupported CRISM file type, supported types are: DDR, MRDR, and TRDR";
      throw IException(IException::User, msg, _FILEINFO_);
    }

    // Translate the Instrument group
    FileName transFile("$ISISROOT/appdata/translations/MroCrismInstrument.trn");
    PvlToPvlTranslationManager instrumentXlater(labelPvl, transFile.expanded());
    instrumentXlater.Auto(outLabel);

    // Translate the Archive group
    transFile  = "$ISISROOT/appdata/translations/MroCrismArchive.trn";
    PvlToPvlTranslationManager archiveXlater(labelPvl, transFile.expanded());
    archiveXlater.Auto(outLabel);

    ocube->putGroup(outLabel.findGroup("Instrument", Pvl::Traverse));
    ocube->putGroup(outLabel.findGroup("Archive", Pvl::Traverse));
    ocube->putGroup(outLabel.findGroup("Kernels", Pvl::Traverse));

    p.StartProcess();
    p.EndProcess();

    results.setName("Results");
    results += PvlKeyword("Warning",
                          "When using cam2map or cam2cam, images imported into "
                          "Isis3 using crism2isis should only be interpolated "
                          "using the nearest-neighbor algorithm due to gimble "
                          "jitter of the MRO CRISM instrument.");
    if (log){
      log->addGroup(results);
    }
    return;
  }
}
+11 −0
Original line number Diff line number Diff line
#ifndef crism2isis_h
#define crism2isis_h

#include "Cube.h"
#include "UserInterface.h"

namespace Isis{
  extern void crism2isis(UserInterface &ui, Pvl *log=nullptr);
}

#endif
+13 −158

File changed.

Preview size limit exceeded, changes collapsed.

Loading