Make sure MISR outputs are masked for night columns

[brhillman]

Previously, the fq_MISR_TAU_v_CTH output variable from the MISR_COLUMN routine was left unmasked for night columns, and due to a typo most of misr_meanztop was left unmasked as well. The result of this was that the output joint histogram variable clMISR and the output cloud top height misr_meanztop would be left unmasked and instead misleadingly set to 0 for night columns. Since MISR retrievals should be set to fillvalue for night columns, we need to make sure all MISR outputs are masked appropriately when sunlit == 0. This PR applies the fix for these two fields in MISR_COLUMN, and also makes sure that output variables from the MISR_SUBCOLUMN routine are similarly masked (previously only dist_model_layertops was masked, while box_MISR_ztop and tauOUT were left unmasked. Fixes #30

[brhillman]

Note that because the example input data does not include any night columns (sunlit == 1 for all columns in the input file), the results of the regression test suggest that this change is BFB with master:

python test_cosp2Imp.py ../../fix-misr-sunlit-mask/driver/data/outputs/UKMO/cosp2_output_um.nc data/outputs/UKMO/cosp2_output_um.nc
############################################################################################
Treating relative differences less than 0.0010000000% as insignificant
All fields match

However, to test that the fix really does affect the output variable for night columns, I created a new input file with the sunlit variable set to 0 for the first 20 columns via

ncap2 -s "sunlit(0:19) = 0;" cosp_input_um.nc cosp_input_um_sunlit.nc

Re-running both examples (master and this branch) and comparing results, we see that we have indeed affected the clMISR and misr_meanztop outputs:

(python2) [bhillma@ghost-login2 build (fix-misr-sunlit-mask)]$ ( cd ../driver; python2 test_cosp2Imp.py data/outputs/UKMO/cosp2_output_um.nc ${reference} )
############################################################################################
Treating relative differences less than 0.0010000000% as insignificant
  clMISR:              13.07 % of values differ, relative range:  -1.00e+00 to -1.00e+00
  misr_meanztop:       13.07 % of values differ, relative range:  -1.00e+00 to -1.00e+00
All other fields match.

We can also plot the differences to observe them visually. We should expect to be able to recover misr_cldarea from the clMISR output by summing all optical depth and cloud top height bins. Since misr_cldarea was set to R_UNDEF in the code previously, we should see a difference when performing this exercise on master and the new branch (ran with the non-trivial inputdata where the first 20 values have been set to night time). The following python code can be used to confirm this:

import xarray as xr
from matplotlib import pyplot

# Open output data
cases = ('master', 'fix-misr-sunlit-mask')
outputs = {
    'master': '~bhillma/codes/cosp/branches/master/driver/data/outputs/UKMO/cosp2_output_um.nc',
    'fix-misr-sunlit-mask': '~bhillma/codes/cosp/branches/fix-misr-sunlit-mask/driver/data/outputs/UKMO/cosp2_output_um.nc'
}
datasets = {c: xr.open_dataset(outputs[c]) for c in cases}

# Plot MISR cloud area for example case
figure, axes = pyplot.subplots(len(cases), 1, sharex=True, sharey=True)
for icase, case in enumerate(cases):
    ax = figure.add_axes(axes[icase])

    # Select data
    clmisr = datasets[case]['clMISR']
    misr_cldarea = datasets[case]['misr_cldarea']

    # We need to mask the data ourselves, since the netcdf files do not define a fillvalue
    clmisr = clmisr.where(clmisr >= 0)
    misr_cldarea = misr_cldarea.where(misr_cldarea >= 0)

    # Calculate cloud area; note that we need to convert misr_cldarea to percent for comparison
    misr_cldarea = 100 * misr_cldarea
    cldarea_from_clmisr = clmisr.sum(dim=('hgt16', 'tau7'), skipna=False)

    # Plot both versions of cloud area
    pl = ax.plot(misr_cldarea, label='misr_cldarea')
    pl = ax.plot(cldarea_from_clmisr, label='cldarea from clmisr', linestyle='dashed')
    ax.set_title('branch: %s'%(case))

    # Add legend to plot
    pyplot.legend(loc='best')

The resulting figure shows that on master, clMISR is not being properly masked, while on the new branch clMISR is being masked and the resulting cloud fraction is properly masked rather than zero where sunlit is 0:

We can also plot the misr_meanztop output variable from both master and the PR branch here:

figure, ax = pyplot.subplots(1, 1, sharex=True, sharey=True)
for icase, case in enumerate(cases):

    # Select data
    data = datasets[case]['misr_meanztop']

    # We need to mask the data ourselves, since the netcdf files do not define a fillvalue
    data = data.where(data >= 0)

    # Plot
    pl = ax.plot(data, label=case)

ax.set_xlabel('Column')
ax.set_ylabel('%s (%s)'%(data.long_name, data.units))
pyplot.legend(loc='best')

and we see that the output is now correctly masked for the first 20 columns, while it was not in master:

[brhillman]

@dustinswales I think this is ready for your eyes. Would you like me to add the input file with the non-trivial sunlit flag to designate some night columns to this PR to properly test masking, per discussion in #32 ?

[dustinswales]

@brhillman Thanks for this. I'll take a gander now. It would be fantastic if you could modify the sunlit field in driver/data/inputs/UKMO/cosp_input_um.nc to be =0 for a few points. The need for some nighttime points in this file is long overdue. The sample file has 1236 points, looks like a single overpass by a polar orbiting satellite. Maybe it makes sense to set half of the sunlit flags to 0? The driver, driver/src/cosp2_test.f90, as it stands only uses the first 153 points, but we can extend the default to cover all of the points 618 daytime, 618 nighttime?

[brhillman]

@dustinswales yeah that'd all be easy to add here. Question is, do you want me to make the change to process all 1236 points here, or in a separate PR? It would make the diffs look pretty obnoxious for this PR, but maybe that's okay?

[dustinswales]

@brhillman Maybe let's do it in a separate PR. Cleaner. And easier to dig up in the future.

Also, I just checked out your code and ran the tests. All is good. Let's wait for Roj's input on this before merging.

[brhillman]

@dustinswales okay cool. Want me to open another then to add the sunlit points?

[dustinswales]

@brhillman That would be fantastic.

[dustinswales]

@roj222