started cleaning up the code a little

main.py

@@ -37,15 +37,55 @@ _eco_file = 'goge1_2_img.v1'
 # #################################################################### #
 
 
-def main(*, data_path=_datapath, mod02=_fname_mod02, mod03=_fname_mod03,
-         img02=_fname_img02, img03=_fname_img03, threshold_file=_threshold_file,
-         geos_atm_1=_geos_atm_1, geos_atm_2=_geos_atm_2, geos_land=_geos_land,
-         geos_ocean=_geos_ocean, geos_constants=_geos_constants, ndvi_file=_ndvi_file, sst_file=_sst_file):
-
-    # datapath = '/ships19/hercules/pveglio/neige_data/snpp_test_input'
-    # fname_l1b = 'VNP02MOD.A2014213.1548.001.2017301015346.uwssec.bowtie_restored_scaled.nc'
-    # fname_geo = 'VNP03MOD.A2014213.1548.001.2017301015705.uwssec.nc'
-    # thresh_file = '/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.yaml'
+def main(*,
+         data_path: str = _datapath,
+         mod02: str = _fname_mod02,
+         mod03: str = _fname_mod03,
+         img02: str = _fname_img02,
+         img03: str = _fname_img03,
+         threshold_file: str = _threshold_file,
+         geos_atm_1: str = _geos_atm_1,
+         geos_atm_2: str = _geos_atm_2,
+         geos_land: str = _geos_land,
+         geos_ocean: str = _geos_ocean,
+         geos_constants: str = _geos_constants,
+         ndvi_file: str = _ndvi_file,
+         sst_file: str = _sst_file) -> None:
+    """Main function for the MVCM.
+
+    Parameters
+    ----------
+    data_path: str
+        path where the data is stored
+    mod02: str
+        VIIRS MOD02 file name
+    mod03: str
+        VIIRS MOD03 file name
+    img02: str
+        VIIRS IMG02 file name
+    img03: str
+        VIIRS IMG03 file name
+    threshold_file: str
+        thresholds file name
+    geos_atm_1: str
+        file name of first GEOS5 file for atmospheric parameters
+    goes_atm_2: str
+        file name of second GEOS5 file for atmospheric parameters
+    geos_land: str
+        file name of GEOS5 land parameters
+    geos_ocean: str
+        file name of GEOS5 ocean parameters
+    geos_constants: str
+        file name of GEOS5 constants
+    ndvi_file: str
+        NDVI file name
+    sst_file:
+        Sea surface temperature file name
+
+    Returns
+    -------
+    None
+    """
 
     file_names = {'MOD02': f'{mod02}',
                   'MOD03': f'{mod03}',
@@ -69,36 +109,21 @@ def main(*, data_path=_datapath, mod02=_fname_mod02, mod03=_fname_mod03,
 
     viirs_data = rd.get_data(file_names, sunglint_angle)
 
-    # scene_xr = xr.Dataset()
-    # for s in scn._scene_list:
-    #    scene_xr[s] = (('number_of_lines', 'number_of_pixels'), scn.scene_id[s])
-    # scene_xr['latitude'] = viirs_xr.latitude
-    # scene_xr['longitude'] = viirs_xr.longitude
-    #
-    # viirs_data = xr.Dataset(viirs_xr, coords=scene_xr)
-    # viirs_data.drop_vars(['latitude', 'longitude'])
-
+    # Initialize the confidence arrays for the various test groups
     cmin_G1 = np.ones(viirs_data.M01.shape)
     cmin_G2 = np.ones(viirs_data.M01.shape)
     cmin_G3 = np.ones(viirs_data.M01.shape)
     cmin_G4 = np.ones(viirs_data.M01.shape)
     cmin_G5 = np.ones(viirs_data.M01.shape)
 
-#    cmin_test = {'Ocean_Day': np.ones(viirs_data.M01.shape),
-#                 'Polar_Ocean_Day': np.ones(viirs_data.M01.shape),
-#                 'Polar_Ocean_Night': np.ones(viirs_data.M01.shape)
-#                 }
-#    cmin2 = np.ones(viirs_data.M01.shape)
-#    cmin3 = np.ones(viirs_data.M01.shape)
-#    cmin4 = np.ones(viirs_data.M01.shape)
-
     # ------------------- #
     # ### Testing Setup ###
     # ------------------- #
-    perform = {'11um BT Test': False,
+
+    perform = {'11um BT Test': True,
                'CO2 High Clouds Test': False,
                'Water Vapor High Clouds Test': False,
-               'Surface Temperature Test': True,
+               'Surface Temperature Test': False,
                'SST Test': False,
                '8.6-11um BT Difference Test': False,
                '11-12um BTD Thin Cirrus Test': False,

read_data.py

@@ -4,11 +4,15 @@ import numpy as np
 import ancillary_data as anc
 import scene as scn
 
+from typing import Dict
+
 _DTR = np.pi/180.
 _RTD = 180./np.pi
 
 
-def read_data(sensor: str, l1b_filename: str, geo_filename: str):
+def read_data(sensor: str,
+              l1b_filename: str,
+              geo_filename: str) -> xr.Dataset:
 
     in_data = xr.open_dataset(geo_filename, group='geolocation_data')
 
@@ -41,12 +45,15 @@ def read_data(sensor: str, l1b_filename: str, geo_filename: str):
     return in_data
 
 
-def relative_azimuth_angle(sensor_azimuth: float, solar_azimuth: float) -> float:
+def relative_azimuth_angle(sensor_azimuth: np.ndarray,
+                           solar_azimuth: np.ndarray) -> np.ndarray:
     rel_azimuth = np.abs(180. - np.abs(sensor_azimuth - solar_azimuth))
     return rel_azimuth
 
 
-def sun_glint_angle(sensor_zenith: float, solar_zenith: float, rel_azimuth: float) -> float:
+def sun_glint_angle(sensor_zenith: np.ndarray,
+                    solar_zenith: np.ndarray,
+                    rel_azimuth: np.ndarray) -> np.ndarray:
     cossna = (np.sin(sensor_zenith*_DTR) * np.sin(solar_zenith*_DTR) * np.cos(rel_azimuth*_DTR) +
               np.cos(sensor_zenith*_DTR) * np.cos(solar_zenith*_DTR))
     cossna[cossna > 1] = 1
@@ -54,7 +61,9 @@ def sun_glint_angle(sensor_zenith: float, solar_zenith: float, rel_azimuth: floa
     return sunglint_angle
 
 
-def scattering_angle(solar_zenith, sensor_zenith, relative_azimuth):
+def scattering_angle(solar_zenith: np.ndarray,
+                     sensor_zenith: np.ndarray,
+                     relative_azimuth: np.ndarray) -> np.ndarray:
     cos_scatt_angle = -1. * (np.cos(solar_zenith*_DTR) * np.cos(sensor_zenith*_DTR) -
                              np.sin(solar_zenith*_DTR) * np.sin(sensor_zenith*_DTR) *
                              np.cos(relative_azimuth*_DTR))
@@ -66,7 +75,9 @@ def correct_reflectances():
     pass
 
 
-def read_ancillary_data(file_names: str, viirs_data: float, resolution=1):
+def read_ancillary_data(file_names: str,
+                        viirs_data: xr.Dataset,
+                        resolution: int = 1) -> xr.Dataset:
 
     # Ancillary files temporarily defined here. Eventually we will find a better way to pass these
     start_time = '2022-06-22 14:54:00.000'
@@ -81,7 +92,7 @@ def read_ancillary_data(file_names: str, viirs_data: float, resolution=1):
 
     dim1 = viirs_data.latitude.shape[0]
     dim2 = viirs_data.latitude.shape[1]
-    print(dim1, dim2)
+
     sst = np.empty((dim1*dim2, ), dtype=np.float32)
     ndvi = np.empty((dim1*dim2, ), dtype=np.float32)
     eco = np.empty((dim1*dim2, ), dtype=np.ubyte)
@@ -97,22 +108,22 @@ def read_ancillary_data(file_names: str, viirs_data: float, resolution=1):
     sst = anc.py_get_Reynolds_SST(viirs_data.latitude.values.reshape((dim1*dim2, )),
                                   viirs_data.longitude.values.reshape((dim1*dim2, )),
                                   resolution, anc_dir, sst_file, sst)
-    print('sst')
+    print('SST read successfully')
 
     ndvi = anc.py_get_NDVI_background(viirs_data.latitude.values.reshape((dim1*dim2, )),
                                       viirs_data.longitude.values.reshape((dim1*dim2, )),
                                       resolution, anc_dir, ndvi_file, ndvi)
+    print('NDVI read successfully')
 
-    print('ndvi')
     eco = anc.py_get_Olson_eco(viirs_data.latitude.values.reshape((dim1*dim2, )),
                                viirs_data.longitude.values.reshape((dim1*dim2, )),
                                resolution, anc_dir, eco)
+    print('Olson eco read successfully')
 
-    print('eco')
     geos_data = anc.py_get_GEOS(viirs_data.latitude.values.reshape((dim1*dim2, )),
                                 viirs_data.longitude.values.reshape((dim1*dim2, )),
                                 resolution, start_time, anc_dir, geos1, geos2, geos_lnd, geos_ocn, geos_cnst, geos_data)
-    print('geos')
+    print('GEOS5 data read successfully')
 
     data = viirs_data
     data['sst'] = (('number_of_lines', 'number_of_pixels'), sst.reshape((dim1, dim2)))
@@ -128,7 +139,8 @@ def read_ancillary_data(file_names: str, viirs_data: float, resolution=1):
     return data
 
 
-def get_data(file_names, sunglint_angle):
+def get_data(file_names: Dict[str, str],
+             sunglint_angle: float) -> xr.Dataset:
 
     mod02 = file_names['MOD02']
     mod03 = file_names['MOD03']
@@ -136,7 +148,7 @@ def get_data(file_names, sunglint_angle):
     viirs_data = read_data('viirs', f'{mod02}', f'{mod03}')
     viirs_data = read_ancillary_data(file_names, viirs_data)
 
-    # Include channel differences and ratios that are used in the tests
+    # Compute channel differences and ratios that are used in the tests
     viirs_data['M15-M14'] = (('number_of_lines', 'number_of_pixels'),
                              viirs_data.M15.values - viirs_data.M14.values)
     viirs_data['M15-M16'] = (('number_of_lines', 'number_of_pixels'),