improved read_data readability. added satellite/sensor in main input arguments

main.py

@@ -38,12 +38,13 @@ _eco_file = 'goge1_2_img.v1'
 # #################################################################### #
 
 
-def main(*,
+def main(satellite: str = 'snpp',
+         sensor: str = 'viirs',
          data_path: str = _datapath,
          mod02: str = _fname_mod02,
          mod03: str = _fname_mod03,
-         img02: str = _fname_img02,
-         img03: str = _fname_img03,
+         img02: str = None,
+         img03: str = None,
          threshold_file: str = _threshold_file,
          geos_atm_1: str = _geos_atm_1,
          geos_atm_2: str = _geos_atm_2,
@@ -56,15 +57,19 @@ def main(*,
 
     Parameters
     ----------
+    satellite: str
+        satellite name, not case-sensitive. Available options: [snpp, ]
+    sensor: str
+        sensor name, not case-sensitive. Available options: [viirs, ]
     data_path: str
         path where the data is stored
     mod02: str
         VIIRS MOD02 file name
     mod03: str
         VIIRS MOD03 file name
-    img02: str
+    img02: [optional] str
         VIIRS IMG02 file name
-    img03: str
+    img03: [optional] str
         VIIRS IMG03 file name
     threshold_file: str
         thresholds file name
@@ -108,7 +113,7 @@ def main(*,
 
     sunglint_angle = thresholds['Sun_Glint']['bounds'][3]
 
-    viirs_data = rd.get_data(file_names, sunglint_angle)
+    viirs_data = rd.get_data(satellite, sensor, file_names, sunglint_angle)
 
     # Initialize the confidence arrays for the various test groups
     cmin_G1 = np.ones(viirs_data.M01.shape)

read_data.py

@@ -5,6 +5,7 @@ import ancillary_data as anc
 import scene as scn
 
 import os
+import logging
 from datetime import datetime as dt
 from typing import Dict
 
@@ -12,30 +13,47 @@ _DTR = np.pi/180.
 _RTD = 180./np.pi
 _bad_data = -999.0
 
+logging.basicConfig(level=logging.INFO, format='%(name)s - %(levelname)s - %(message)s')
+# logging.basicConfig(level=logging.INFO, filename='logfile.log', 'filemode='w',
+#                       format='%(name)s %(levelname)s %(message)s')
 
-def read_data(sensor: str,
-              l1b_filename: str,
-              geo_filename: str) -> xr.Dataset:
+
+def read_viirs_data(l1b_filename: str,
+                    geo_filename: str) -> xr.Dataset:
+    """Read VIIRS MOD or IMG data
+
+    Parameters
+    ----------
+    l1b_filename: str
+        L1b VIIRS filename
+    geo_filename: str
+        geoocation VIIRS filename
+
+    Returns
+    -------
+    in_data: xarray.Dataset
+        dataset containing VIIRS bands, geolocation data and additional info (i.e. relative azimuth,
+        sunglint angle and scattering angle)
+    """
 
     in_data = xr.open_dataset(geo_filename, group='geolocation_data')
 
     data = xr.open_dataset(l1b_filename, group='observation_data', decode_cf=False)
 
-    if sensor.lower() == 'viirs':
-        for band in list(data.variables):
-            if 'reflectance' in data[band].long_name:
-                if hasattr(data[band], 'VCST_scale_factor'):
-                    scale_factor = data[band].VCST_scale_factor * data[band].bias_correction
-                else:
-                    scale_factor = data[band].scale_factor
-                in_data[band] = (('number_of_lines', 'number_of_pixels'),
-                                 data[band].values * scale_factor / np.cos(in_data.solar_zenith.values*_DTR))
-
-            elif 'radiance' in data[band].long_name:
-                in_data[band] = (('number_of_lines', 'number_of_pixels'),
-                                 data[f'{band}_brightness_temperature_lut'].values[data[band].values])
+    for band in list(data.variables):
+        if 'reflectance' in data[band].long_name:
+            if hasattr(data[band], 'VCST_scale_factor'):
+                scale_factor = data[band].VCST_scale_factor * data[band].bias_correction
             else:
-                pass
+                scale_factor = data[band].scale_factor
+            in_data[band] = (('number_of_lines', 'number_of_pixels'),
+                             data[band].values * scale_factor / np.cos(in_data.solar_zenith.values*_DTR))
+
+        elif 'radiance' in data[band].long_name:
+            in_data[band] = (('number_of_lines', 'number_of_pixels'),
+                             data[f'{band}_brightness_temperature_lut'].values[data[band].values])
+        else:
+            pass
 
     relazi = relative_azimuth_angle(in_data.sensor_azimuth.values, in_data.solar_azimuth.values)
     sunglint = sun_glint_angle(in_data.sensor_zenith.values, in_data.solar_zenith.values, relazi)
@@ -50,6 +68,19 @@ def read_data(sensor: str,
 
 def relative_azimuth_angle(sensor_azimuth: np.ndarray,
                            solar_azimuth: np.ndarray) -> np.ndarray:
+    """Computation of the relative azimuth angle
+
+    Parameters
+    ----------
+    sensor_azimuth: np.ndarray
+        sensor azimuth angle from the geolocation file
+    solar_azimuth: np.ndarray
+        solar azimuth angle from the geolocation file
+
+    Returns
+    -------
+    relative_azimuth: np.ndarray
+    """
     rel_azimuth = np.abs(180. - np.abs(sensor_azimuth - solar_azimuth))
     return rel_azimuth
 
@@ -57,6 +88,21 @@ def relative_azimuth_angle(sensor_azimuth: np.ndarray,
 def sun_glint_angle(sensor_zenith: np.ndarray,
                     solar_zenith: np.ndarray,
                     rel_azimuth: np.ndarray) -> np.ndarray:
+    """Computation of the sun glint angle
+
+    Parameters
+    ----------
+    sensor_zenith: np.ndarray
+        sensor zenith angle from the geolocation file
+    solar_zenith: np.ndarray
+        solar zenith angle from the geolocation file
+    relative_azimuth: np.ndarray
+        relative azimuth computed from function relative_azimuth_angle()
+
+    Returns
+    -------
+    sunglint_angle: 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
@@ -67,6 +113,21 @@ def sun_glint_angle(sensor_zenith: np.ndarray,
 def scattering_angle(solar_zenith: np.ndarray,
                      sensor_zenith: np.ndarray,
                      relative_azimuth: np.ndarray) -> np.ndarray:
+    """Computation of the scattering angle
+
+    Parameters
+    ----------
+    solar_zenith: np.ndarray
+        solar zenith angle from the geolocation file
+    sensor_zenith: np.ndarray
+        sensor zenith angle angle from the geolocation file
+    relative_azimuth: np.ndarray
+        relative azimuth computed from function relative_azimuth_angle()
+
+    Returns
+    -------
+    scattering_angle: 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))
@@ -79,9 +140,22 @@ def correct_reflectances():
 
 
 def read_ancillary_data(file_names: str,
-                        viirs_data: xr.Dataset,
+                        latitude: np.ndarray,
+                        longitude: np.ndarray,
                         resolution: int = 1) -> xr.Dataset:
-
+    """Read ancillary data using C functions from original MVCM
+
+    Parameters
+    ----------
+    file_names: str
+    latitude: np.ndarray
+    longitude: np.ndarray
+    resolution: int
+
+    Returns
+    -------
+    data: xarray.Dataset
+    """
     # Ancillary files temporarily defined here. Eventually we will find a better way to pass these
     anc_dir = file_names['ANC_DIR']
     sst_file = file_names['SST']
@@ -93,58 +167,51 @@ def read_ancillary_data(file_names: str,
     geos2 = file_names['GEOS_atm_2']
     vnptime = '.'.join(os.path.basename(file_names['MOD02']).split('.')[1:3])
     start_time = dt.strftime(dt.strptime(vnptime, 'A%Y%j.%H%M'), '%Y-%m-%d %H:%M:%S.000')
-    # start_time = '2022-06-22 14:54:00.000'
 
-    dim1 = viirs_data.latitude.shape[0]
-    dim2 = viirs_data.latitude.shape[1]
+    out_shape = latitude.shape
 
-    sst = np.empty((dim1*dim2, ), dtype=np.float32)
-    ndvi = np.empty((dim1*dim2, ), dtype=np.float32)
-    eco = np.empty((dim1*dim2, ), dtype=np.ubyte)
+    sst = np.empty((np.prod(out_shape), ), dtype=np.float32)
+    ndvi = np.empty((np.prod(out_shape), ), dtype=np.float32)
+    eco = np.empty((np.prod(out_shape), ), dtype=np.ubyte)
 
-    geos_data = {'tpw': np.empty((dim1*dim2, ), dtype=np.float32),
-                 'snowfr': np.empty((dim1*dim2, ), dtype=np.float32),
-                 'icefr': np.empty((dim1*dim2, ), dtype=np.float32),
-                 'ocnfr': np.empty((dim1*dim2, ), dtype=np.float32),
-                 'landicefr': np.empty((dim1*dim2, ), dtype=np.float32),
-                 'sfct': np.empty((dim1*dim2, ), dtype=np.float32),
+    geos_data = {'tpw': np.empty((np.prod(out_shape), ), dtype=np.float32),
+                 'snowfr': np.empty((np.prod(out_shape), ), dtype=np.float32),
+                 'icefr': np.empty((np.prod(out_shape), ), dtype=np.float32),
+                 'ocnfr': np.empty((np.prod(out_shape), ), dtype=np.float32),
+                 'landicefr': np.empty((np.prod(out_shape), ), dtype=np.float32),
+                 'sfct': np.empty((np.prod(out_shape), ), dtype=np.float32),
                  }
 
-    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 read successfully')
+    sst = anc.py_get_Reynolds_SST(latitude.ravel(), longitude.ravel(), resolution, anc_dir, sst_file, sst)
+    logging.info('SST read successfully')
 
-    ndvi = anc.py_get_NDVI_background(viirs_data.latitude.values.reshape((dim1*dim2, )),
-                                      viirs_data.longitude.values.reshape((dim1*dim2, )),
+    ndvi = anc.py_get_NDVI_background(latitude.ravel(), longitude.ravel(),
                                       resolution, anc_dir, ndvi_file, ndvi)
-    print('NDVI read successfully')
-
-    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')
-
-    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('GEOS5 data read successfully')
-
-    data = viirs_data
-    data['sst'] = (('number_of_lines', 'number_of_pixels'), sst.reshape((dim1, dim2)))
-    data['ndvi'] = (('number_of_lines', 'number_of_pixels'), ndvi.reshape((dim1, dim2)))
-    data['eco'] = (('number_of_lines', 'number_of_pixels'), eco.reshape((dim1, dim2)))
-    data['geos_tpw'] = (('number_of_lines', 'number_of_pixels'), geos_data['tpw'].reshape((dim1, dim2)))
-    data['geos_snowfr'] = (('number_of_lines', 'number_of_pixels'), geos_data['snowfr'].reshape((dim1, dim2)))
-    data['geos_icefr'] = (('number_of_lines', 'number_of_pixels'), geos_data['icefr'].reshape((dim1, dim2)))
-    data['geos_ocnfr'] = (('number_of_lines', 'number_of_pixels'), geos_data['ocnfr'].reshape((dim1, dim2)))
-    data['geos_landicefr'] = (('number_of_lines', 'number_of_pixels'), geos_data['landicefr'].reshape((dim1, dim2)))
-    data['geos_sfct'] = (('number_of_lines', 'number_of_pixels'), geos_data['sfct'].reshape((dim1, dim2)))
+    logging.info('NDVI read successfully')
+
+    eco = anc.py_get_Olson_eco(latitude.ravel(), longitude.ravel(), resolution, anc_dir, eco)
+    logging.info('Olson eco read successfully')
+
+    geos = anc.py_get_GEOS(latitude.ravel(), longitude.ravel(), resolution, start_time, anc_dir,
+                           geos1, geos2, geos_lnd, geos_ocn, geos_cnst, geos_data)
+    logging.info('GEOS5 data read successfully')
+
+    ancillary = {'sst': sst.reshape(out_shape),
+                 'ndvi': ndvi.reshape(out_shape),
+                 'eco': eco.reshape(out_shape)
+                 }
+    for var in list(geos):
+        ancillary[f'geos_{var}'] = geos[var].reshape(out_shape)
+
+    dims = ('number_of_lines', 'number_of_pixels')
+    data = xr.Dataset.from_dict({var: {'dims': dims, 'data': ancillary[var]} for var in list(ancillary)})
 
     return data
 
 
-def get_data(file_names: Dict[str, str],
+def get_data(satellite: str,
+             sensor: str,
+             file_names: Dict[str, str],
              sunglint_angle: float,
              hires: bool = False) -> xr.Dataset:
 
@@ -156,7 +223,7 @@ def get_data(file_names: Dict[str, str],
         img03 = file_names['IMG03']
 
     if hires is False:
-        viirs_data = read_data('viirs', f'{mod02}', f'{mod03}')
+        viirs_data = read_viirs_data(sensor, f'{mod02}', f'{mod03}')
         viirs_data = read_ancillary_data(file_names, viirs_data)
 
         if (('M05' in viirs_data) and ('M07' in viirs_data)):

setup.py

@@ -26,14 +26,14 @@ sourcefiles = ['src/get_Reynolds_SST.c',      # SST
                ]
 
 include_dirs = ['include',
-                '/opt/hdf4/4.2.9-gcc-4.9.2/include',
-                '/opt/hdfeos2/2.19-gcc-4.9.2/include',
-                '/opt/netcdf4/4.3.3-gcc-4.9.2/include',
+                '/opt/hdf4/4.2.14-gcc-8.3/include',
+                '/opt/hdfeos2/2.20-gcc-8.3/include',
+                '/opt/netcdf4/4.7.0-gcc-8.3/include',
                 numpy.get_include(), ]
 
-library_dirs = ['/opt/hdf4/4.2.9-gcc-4.9.2/lib',
-                '/opt/hdfeos2/2.19-gcc-4.9.2/lib',
-                '/opt/netcdf4/4.3.3-gcc-4.9.2/lib',
+library_dirs = ['/opt/hdf4/4.2.14-gcc-8.3/lib',
+                '/opt/hdfeos2/2.20-gcc-8.3/lib',
+                '/opt/netcdf4/4.7.0-gcc-8.3/lib',
                 ]
 
 extensions = [Extension('ancillary_data', sourcefiles,