search

nasa / HyperCP

Other
41 stars 22 forks source link

PIU Object too deep #251

Closed oceancolorcoder closed 4 weeks ago

oceancolorcoder commented 1 month ago

@ARamsay17 this is a new bug in PIU processing new pySAS data image

Uncertainty Update Elapsed Time: 11.0 s
Perform simple residual NIR subtraction.
Processing MODISA
Convolving MODIS Aqua (ir)radiances in the slice
28 spectra in slice (ensemble).
3 spectra remaining in slice to average after filtering to lowest 10.0%.
Calculating M99 glint correction with complete LUT
Calculating Zhang glint correction.
Zhang17 Elapsed Time: 4.1 s
Reading : Data/hybrid_reference_spectrum_p1nm_resolution_c2020-09-21_with_unc.nc
/opt/anaconda3/envs/hypercp/lib/python3.11/site-packages/comet_maths/linear_algebra/matrix_calculation.py:274: UserWarning: One of the correlation matrices is not positive definite. It has been slightly changed (maximum difference of 6.661338147750939e-16) to accomodate our method.
  warnings.warn(
/opt/anaconda3/envs/hypercp/lib/python3.11/site-packages/comet_maths/linear_algebra/matrix_calculation.py:274: UserWarning: One of the correlation matrices is not positive definite. It has been slightly changed (maximum difference of 4.440892098500626e-16) to accomodate our method.
  warnings.warn(
/opt/anaconda3/envs/hypercp/lib/python3.11/site-packages/comet_maths/linear_algebra/matrix_calculation.py:274: UserWarning: One of the correlation matrices is not positive definite. It has been slightly changed (maximum difference of 6.661338147750939e-16) to accomodate our method.
  warnings.warn(
Uncertainty Update Elapsed Time: 11.1 s
Perform simple residual NIR subtraction.
Processing MODISA
Traceback (most recent call last):
  File "/Users/daurin/GitRepos/HyperCP/Main.py", line 523, in singleL2Clicked
    self.processSingle("L2")
  File "/Users/daurin/GitRepos/HyperCP/Main.py", line 491, in processSingle
    Controller.processFilesSingleLevel(
  File "/Users/daurin/GitRepos/HyperCP/Source/Controller.py", line 897, in processFilesSingleLevel
    Controller.processSingleLevel(pathOut, fp, calibrationMap, level, flag_Trios)
  File "/Users/daurin/GitRepos/HyperCP/Source/Controller.py", line 757, in processSingleLevel
    root = Controller.processL2(root,outFilePath)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/daurin/GitRepos/HyperCP/Source/Controller.py", line 448, in processL2
    node = ProcessL2.processL2(root,station)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/daurin/GitRepos/HyperCP/Source/ProcessL2.py", line 2364, in processL2
    if not ProcessL2.stationsEnsemblesReflectance(node, root,station):
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/daurin/GitRepos/HyperCP/Source/ProcessL2.py", line 2253, in stationsEnsemblesReflectance
    if not ProcessL2.ensemblesReflectance(node, sasGroup, referenceGroup, ancGroup, 
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/daurin/GitRepos/HyperCP/Source/ProcessL2.py", line 1362, in ensemblesReflectance
    stats = instrument.generateSensorStats("SeaBird",
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/daurin/GitRepos/HyperCP/Source/ProcessInstrumentUncertainties.py", line 90, in generateSensorStats
    _, output[stype]['std_Signal_Interpolated'] = self.interp_common_wvls(
                                                  ^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/daurin/GitRepos/HyperCP/Source/ProcessInstrumentUncertainties.py", line 1461, in interp_common_wvls
    new_y = np.interp(newWavebands, x, y)  #InterpolatedUnivariateSpline(x, y, k=3)(newWavebands)
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "<__array_function__ internals>", line 180, in interp
  File "/opt/anaconda3/envs/hypercp/lib/python3.11/site-packages/numpy/lib/function_base.py", line 1594, in interp
    return interp_func(x, xp, fp, left, right)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
           ValueError: object too deep for desired array
oceancolorcoder commented 1 month ago

@ARamsay17 Data for testing provided in Data/PACE-PAX

oceancolorcoder commented 4 weeks ago

@ARamsay17 I still can’t run PIU on PACE-PAX data on the Box. During processing of the fifth ensemble while running ProcessL2.ensemblesReflectance, we get to L1363:

 stats = instrument.generateSensorStats("SeaBird",
                    dict(ES=esRawGroup, LI=liRawGroup, LT=ltRawGroup),
                    dict(ES=esRawSlice, LI=liRawSlice, LT=ltRawSlice),
                    instrument_wb)

This calls ProcessInstrumentUncertainties.lightDarkStats for HyperOCR:

           # rawData here is the group, passed along only for the purpose of
            # confirming "FrameTypes", i.e., ShutterLight or ShutterDark. Calculations
            # are performed on the Slice.
            # output contains:
            # ave_Light: (array 1 x number of wavebands)
            # ave_Dark: (array 1 x number of wavebands)
            # std_Light: (array 1 x number of wavebands)
            # std_Dark: (array 1 x number of wavebands)
            # std_Signal: OrdDict by wavebands: sqrt( (std(Light)^2 + std(Dark)^2)/ave(Light)^2 )
            output[sensortype] = self.lightDarkStats(
                [rawData[sensortype]['LIGHT'],
                rawData[sensortype]['DARK']],
                [rawSlice[sensortype]['LIGHT'],
                rawSlice[sensortype]['DARK']],
                sensortype
            )

As noted, these should all be 1x arrays plus the dictionary for std_Signal with waveband keys and a single spectrum. However, std_Signal for this ensemble yields waveband keys and an array of spectra. It may be important that this fifth ensemble has 25 spectra (the others are 32, 31, 35, 28) because HyperOCR.lightDarkStats (L1262) reads:

            if N > 25:  # normal case
                    std_Light.append(np.std(lightData[k])/np.sqrt(N))
                    std_Dark.append(np.std(darkData[k])/np.sqrt(Nd) )  # sigma here is essentially sigma**2 so N must sqrt
            else:  # few scans, use different statistics
                std_Light.append((N-1/N-3)*(lightData[k] / np.sqrt(N))**2)
                    std_Dark.append((Nd-1/Nd-3)*(darkData[k] / np.sqrt(Nd))**2)

So, the structural problem is in how the alternate statistics are being calculated for N<=25.

Indeed, the result of the upper clause of the if statement yields a single number each time, while the result of the second clause yields a vector 25 long. This is the result of vectorwise use of lightData[k] in the latter, but scalarwise use of np.std(lightData[k] in the former. Could you link to the intended equation for low N?

(Side note: how is it that the number of spectra within the darkSlice are equivalent to the lightSlice here? Darks are only collected intermittently after every 6 lights or so (see the L1AQC groups in this L1BQC. So, unless darks in the slice variable, slice[1] were interpolated in time to the lights, slice[0], this doesn’t make sense. If they were interpolated, then taking the averages and standard deviations will introduce error, or …?)

oceancolorcoder commented 4 weeks ago

Resolved with PR #263