read_data tested. working as intended

read_data.py

@@ -55,8 +55,8 @@ class CollectInputs(object):
     ancillary_dir: str = field(default='', validator=[validators.instance_of(str), ])
     sst_file: str = field(default='', validator=[validators.instance_of(str), ])
     ndvi_file: str = field(default='', validator=[validators.instance_of(str), ])
-    geos_file_1: str = field(default='', validator=[validators.instance_of(str), ])
-    geos_file_2: str = field(default='', validator=[validators.instance_of(str), ])
+    geos_file_before: str = field(default='', validator=[validators.instance_of(str), ])
+    geos_file_after: str = field(default='', validator=[validators.instance_of(str), ])
     geos_land: str = field(default='', validator=[validators.instance_of(str), ])
     geos_ocean: str = field(default='', validator=[validators.instance_of(str), ])
     geos_constants: str = field(default='', validator=[validators.instance_of(str), ])
@@ -147,11 +147,14 @@ class ReadData(CollectInputs):
         return rad_data
 
     def preprocess_viirs(self,
+                         geo_data: xr.Dataset,
                          viirs: xr.Dataset) -> xr.Dataset:
         """Create band combinations (e.g. ratios, differences) that are used by some spectral tests.
 
         Parameters
         ----------
+        geo_data: xarray.Dataset
+            dataset containing all geolocation data
         viirs: xr.Dataset
             VIIRS L1b data
 
@@ -209,6 +212,9 @@ class ReadData(CollectInputs):
         # temp value to force the code to work
         viirs['M128'] = (self.dims, np.zeros(viirs.M15.shape))
 
+        viirs.update(geo_data)
+        viirs = viirs.set_coords(['latitude', 'longitude'])
+
         logging.debug('Viirs preprocessing completed successfully.')
         return viirs
 
@@ -368,7 +374,7 @@ class ReadAncillary(CollectInputs):
         """
         geos_variables = ['tpw', 'snow_fraction', 'ice_fraction', 'ocean_fraction',
                           'land_ice_fraction', 'surface_temperature']
-        geos_data = {var: np.empty(self.out_shape, dtype=np.float32) for var in geos_variables}
+        geos_data = {var: np.empty(self.out_shape, dtype=np.float32).ravel() for var in geos_variables}
 
         geos_data = anc.py_get_GEOS(self.latitude.ravel(), self.longitude.ravel(), self.resolution,
                                     self.get_granule_time(), self.ancillary_dir, self.geos_file_before,
@@ -381,79 +387,40 @@ class ReadAncillary(CollectInputs):
         logging.debug('GEOS data read successfully')
         return geos_data
 
+    def pack_data(self) -> xr.Dataset:
+        """Pack all the ancillary data into a single dataset
 
-def correct_reflectances():
-    pass
+        Parameters
+        ----------
 
+        Returns
+        -------
+        ancillary_data: xr.Dataset
+            dataset containing all the ancillary data
+        """
+        ancillary_data = xr.Dataset()
+        ancillary_data['sst'] = (self.dims, self.get_sst())
+        ancillary_data['ndvi'] = (self.dims, self.get_ndvi())
 
-def read_ancillary_data(file_names: str,
-                        latitude: np.ndarray,
-                        longitude: np.ndarray,
-                        resolution: int = 1) -> xr.Dataset:
-    """Read ancillary data using C functions from original MVCM
+        geos_tmp = self.get_geos()
+        for var in list(geos_tmp.keys()):
+            ancillary_data[f'geos_{var}'] = (self.dims, geos_tmp[var])
 
-    Parameters
-    ----------
-    file_names: str
-    latitude: np.ndarray
-    longitude: np.ndarray
-    resolution: int
+        return ancillary_data
 
-    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']
-    ndvi_file = file_names['NDVI']
-    geos_cnst = file_names['GEOS_constants']
-    geos_lnd = file_names['GEOS_land']
-    geos_ocn = file_names['GEOS_ocean']
-    geos1 = file_names['GEOS_atm_1']
-    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')
-
-    out_shape = latitude.shape
-
-    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((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(latitude.ravel(), longitude.ravel(), resolution, anc_dir, sst_file, sst)
-    logging.info('SST read successfully')
-
-    ndvi = anc.py_get_NDVI_background(latitude.ravel(), longitude.ravel(),
-                                      resolution, anc_dir, ndvi_file, ndvi)
-    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)})
+    scene_flags = scn.find_scene(viirs, sunglint_angle)
+    scene = scn.scene_id(scene_flags)
+
+    scene_xr = xr.Dataset()
+    for s in scn._scene_list:
+        scene_xr[s] = (('number_of_lines', 'number_of_pixels'), scene[s])
+
+    scene['lat'] = viirs.latitude
+    scene['lon'] = viirs.longitude
+
+    data = xr.Dataset(viirs, coords=scene_xr)
+    data.drop_vars(['latitude', 'longitude'])
 
-    return data
 
 
 def get_data(satellite: str,
@@ -567,22 +534,3 @@ def get_data(satellite: str,
     data.drop_vars(['latitude', 'longitude'])
 
     return data
-
-
-def get_data_new(satellite: str,
-                 sensor: str,
-                 file_names: Dict[str, str],
-                 sunglint_angle: float,
-                 hires: bool = False) -> xr.Dataset:
-    """Reads satellite and ancillary data and packs everything in an xarray dataset for the main function
-    """
-    mod02 = file_names['MOD02']
-    mod03 = file_names['MOD03']
-
-    viirs = read_viirs(sensor, f'{mod02}', f'{mod03}')
-    ancillary_data = read_ancillary_data(file_names, viirs.latitude.values, viirs.longitude.values)
-
-
-
-
-