diff --git a/mvcm/main.py b/mvcm/main.py index 46fb12a4a4675c85a94da902f18ab1ee46cb1702..31871d10f4fc7ca9629c36ccc5fb0495b3b7fe26 100644 --- a/mvcm/main.py +++ b/mvcm/main.py @@ -6,15 +6,15 @@ except ImportError: import ruamel.yaml as yml import argparse - -# import cProfile import functools import logging +import os import time import numpy as np import xarray as xr from pkg_resources import get_distribution # type: ignore +from rich.logging import RichHandler import mvcm.read_data as rd import mvcm.scene as scn @@ -22,10 +22,12 @@ import mvcm.spectral_tests as tst import mvcm.utility_functions as utils from mvcm.restoral import Restoral -# from glob import glob - +_LOG_FORMAT = "%(message)s" +logging.basicConfig(level="NOTSET", datefmt="[%X]", format=_LOG_FORMAT, handlers=[RichHandler()]) +logger = logging.getLogger(__name__) LOG_LEVELS = ["CRITICAL", "ERROR", "WARNING", "INFO", "DEBUG"] + # #################################################################### # # TEST CASE # data: @@ -104,7 +106,7 @@ def timer(func): value = func(*args, **kwargs) end_time = time.time() run_time = end_time - start_time - print(f"Finished {func.__name__!r} in {run_time:.2f} secs\n") + logger.info(f"Finished {func.__name__!r} in {run_time:.2f} secs\n") return value return wrapper_timer @@ -115,10 +117,10 @@ def main( satellite: str = "snpp", sensor: str = "viirs", data_path: str = _datapath, - mod02: str = "", # _fname_mod02, - mod03: str = "", # _fname_mod03, - img02: str = None, # _fname_img02, - img03: str = None, # _fname_img03, + mod02: str = "", + mod03: str = "", + img02: str = "", + img03: str = "", threshold_file: str = _threshold_file, geos_atm_1: str = _geos_atm_1, geos_atm_2: str = _geos_atm_2, @@ -130,7 +132,7 @@ def main( eco_file: str = _eco_file, out_file: str = "temp_out.nc", verbose: int = 0, - debug: bool = False, + debug: bool = False, # noqa ) -> None: """MVCM main function. @@ -170,6 +172,8 @@ def main( Ecosystem File out_file: str Output file name + verbose: int + Verbosity level Returns ------- @@ -178,6 +182,48 @@ def main( verbose_level = np.minimum(verbose, 4) logging.basicConfig(level=LOG_LEVELS[verbose_level]) + if os.path.exists(mod02) is False: + mod02_err = f"File {mod02} not found" + raise FileNotFoundError(mod02_err) + if os.path.exists(mod03) is False: + mod03_err = f"File {mod03} not found" + raise FileNotFoundError(mod03_err) + if img02 != "" and os.path.exists(img02) is False: + img02_err = f"File {img02} not found" + raise FileNotFoundError(img02_err) + if img03 != "" and os.path.exists(img03) is False: + img03_err = f"File {img03} not found" + raise FileNotFoundError(img03_err) + if os.path.exists(geos_atm_1) is False: + geos_atm_1_err = f"File {geos_atm_1} not found" + raise FileNotFoundError(geos_atm_1_err) + if os.path.exists(geos_atm_2) is False: + geos_atm_2_err = f"File {geos_atm_2} not found" + raise FileNotFoundError(geos_atm_2_err) + if os.path.exists(geos_land) is False: + geos_land_err = f"File {geos_land} not found" + raise FileNotFoundError(geos_land_err) + if os.path.exists(geos_ocean) is False: + geos_ocean_err = f"File {geos_ocean} not found" + raise FileNotFoundError(geos_ocean_err) + if os.path.exists(geos_constants) is False: + geos_constants_err = f"File {geos_constants} not found" + raise FileNotFoundError(geos_constants_err) + if os.path.exists(ndvi_file) is False: + ndvi_file_err = f"File {ndvi_file} not found" + raise FileNotFoundError(ndvi_file_err) + if os.path.exists(sst_file) is False: + sst_file_err = f"File {sst_file} not found" + raise FileNotFoundError(sst_file_err) + if os.path.exists(eco_file) is False: + eco_file_err = f"File {eco_file} not found" + raise FileNotFoundError(eco_file_err) + + if img02 != "" or img03 != "": + use_hires = False + else: + use_hires = True + file_names = { "MOD02": f"{mod02}", "MOD03": f"{mod03}", @@ -198,32 +244,29 @@ def main( text = f.read() thresholds = yml.safe_load(text) - ################# - ################# - if img02 is None or img03 is None: - use_hires = False - else: - use_hires = True - ################# - ################# + # We are not processing night granules + if use_hires is True: + with xr.open_dataset(file_names["IMG02"], group="observation_data") as vnp02: + for b in _mod_bands: + try: + vnp02[b] + except KeyError: + logger.info(f"Band {b} not found in file. No output will be written.") + return - # for the time being we're not processing night granules - # vnp02 = xr.open_dataset(file_names["MOD02"], group="observation_data") - # for b in _mod_bands: - # try: - # vnp02[b] - # except KeyError: - # logging.info(f"Band {b} not found in file. No output will be written.") - # return - - logging.info(f"Reading VNP02: {file_names['MOD02']}") + logger.info(f"Reading VNP02: {file_names['MOD02']}") viirs_data = rd.get_data(satellite, sensor, file_names, hires=use_hires) sunglint_angle = thresholds["Sun_Glint"]["bounds"][3] scene_flags = scn.find_scene(viirs_data, sunglint_angle) - Restore = Restoral(data=viirs_data, thresholds=thresholds, scene_flag=scene_flags) + restore = Restoral( + data=viirs_data, + thresholds=thresholds, + scene_flag=scene_flags, + log_level=LOG_LEVELS[verbose_level], + ) cmin_3d = np.ones((18, viirs_data.M11.shape[0], viirs_data.M11.shape[1])) @@ -274,13 +317,14 @@ def main( scene_types = np.zeros(viirs_data.M11.shape) for _scene_i, scene_name in enumerate(_scene_list): scene_types[viirs_data[scene_name].values == 1] = i - print(f"Processing {scene_name}") + + logger.info(f"Processing {scene_name}") if np.all(viirs_data[scene_name].values == 0): - logging.info("Skipping, no pixels in scene.") + logger.info("Skipping, no pixels in scene.") continue - MyScene = tst.CloudTests(data=viirs_data, scene_name=scene_name, thresholds=thresholds) + my_scene = tst.CloudTests(data=viirs_data, scene_name=scene_name, thresholds=thresholds) # Initialize the confidence arrays for the various test groups cmin_g1 = np.ones(viirs_data.M11.shape) @@ -296,41 +340,41 @@ def main( m15_name = "M15" # Group 1 - cmin_g1, bits["01"] = MyScene.test_11um(m15_name, cmin_g1, bits["01"]) + cmin_g1, bits["01"] = my_scene.test_11um(m15_name, cmin_g1, bits["01"]) # this test needs to be implemented properly, # for now is commented since it's only used for night granules - # cmin_g1, bits["02"] = MyScene.surface_temperature_test( + # cmin_g1, bits["02"] = my_scene.surface_temperature_test( # "", viirs_data, cmin_g1, bits["02"] # ) - cmin_g1, bits["03"] = MyScene.sst_test("M15", "M16", cmin_g1, bits["03"]) + cmin_g1, bits["03"] = my_scene.sst_test("M15", "M16", cmin_g1, bits["03"]) # Group 2 - cmin_g2, bits["04"] = MyScene.bt_diff_86_11um("M14-M15", cmin_g2, bits["04"]) - cmin_g2, bits["05"] = MyScene.test_11_12um_diff("M15-M16", cmin_g2, bits["05"]) - cmin_g2, bits["06"] = MyScene.variability_11um_test(m15_name, cmin_g2, bits["06"]) - cmin_g2, bits["07"] = MyScene.bt_difference_11_4um_test_ocean( + cmin_g2, bits["04"] = my_scene.bt_diff_86_11um("M14-M15", cmin_g2, bits["04"]) + cmin_g2, bits["05"] = my_scene.test_11_12um_diff("M15-M16", cmin_g2, bits["05"]) + cmin_g2, bits["06"] = my_scene.variability_11um_test(m15_name, cmin_g2, bits["06"]) + cmin_g2, bits["07"] = my_scene.bt_difference_11_4um_test_ocean( "M15-M12", cmin_g2, bits["07"] ) - cmin_g2, bits["08"] = MyScene.bt_difference_11_4um_test_land( + cmin_g2, bits["08"] = my_scene.bt_difference_11_4um_test_land( "M15-M16", "M15-M12", cmin_g2, bits["08"] ) - cmin_g2, bits["09"] = MyScene.oceanic_stratus_11_4um_test("M15-M12", cmin_g2, bits["09"]) + cmin_g2, bits["09"] = my_scene.oceanic_stratus_11_4um_test("M15-M12", cmin_g2, bits["09"]) # Group 3 - cmin_g3, bits["10"] = MyScene.nir_reflectance_test("M07", cmin_g3, bits["10"]) + cmin_g3, bits["10"] = my_scene.nir_reflectance_test("M07", cmin_g3, bits["10"]) # I need to make sure that this test is performing as intended. # Right now it doesn't seem to have an impact in any way whatsoever - cmin_g3, bits["11"] = MyScene.vis_nir_ratio_test("M07-M05ratio", cmin_g3, bits["11"]) - cmin_g3, bits["12"] = MyScene.test_16_21um_reflectance("M10", cmin_g3, bits["12"]) + cmin_g3, bits["11"] = my_scene.vis_nir_ratio_test("M07-M05ratio", cmin_g3, bits["11"]) + cmin_g3, bits["12"] = my_scene.test_16_21um_reflectance("M10", cmin_g3, bits["12"]) # this test seems to be almost working, but not 100% there yet - cmin_g3, bits["13"] = MyScene.visible_reflectance_test("M05", cmin_g3, bits["13"]) - cmin_g3, bits["11"] = MyScene.gemi_test("GEMI", cmin_g3, bits["11"]) + cmin_g3, bits["13"] = my_scene.visible_reflectance_test("M05", cmin_g3, bits["13"]) + cmin_g3, bits["11"] = my_scene.gemi_test("GEMI", cmin_g3, bits["11"]) # Group 4 - cmin_g4, bits["15"] = MyScene.test_1_38um_high_clouds("M09", cmin_g4, bits["15"]) + cmin_g4, bits["15"] = my_scene.test_1_38um_high_clouds("M09", cmin_g4, bits["15"]) # Group 5 - cmin_g5, bits["16"] = MyScene.thin_cirrus_4_12um_BTD_test("M13-M16", cmin_g5, bits["16"]) + cmin_g5, bits["16"] = my_scene.thin_cirrus_4_12um_BTD_test("M13-M16", cmin_g5, bits["16"]) bit = {} qabit = {} @@ -371,7 +415,7 @@ def main( & (cmin >= 0.05) & (cmin <= 0.99) ) - cmin[idx], bits["r01"]["test"][idx] = Restore.spatial_variability(cmin, idx, bits["r01"]) + cmin[idx], bits["r01"]["test"][idx] = restore.spatial_variability(cmin, idx, bits["r01"]) # Sun Glint sunglint_bits = ( @@ -386,7 +430,7 @@ def main( & (cmin <= 0.95) & (sunglint_bits == 1) ) - cmin[idx], bits["r02"]["test"][idx] = Restore.sunglint(cmin, idx, bits["r02"]) + cmin[idx], bits["r02"]["test"][idx] = restore.sunglint(cmin, idx, bits["r02"]) # Shallow Water # ... @@ -399,7 +443,7 @@ def main( idx = np.nonzero( (scene_flags["day"] == 1) & (scene_flags["water"] == 1) & (scene_flags["ice"] == 0) ) - cmin[idx], bits["r03"]["test"][idx] = Restore.shallow_water( + cmin[idx], bits["r03"]["test"][idx] = restore.shallow_water( sh_water_bits, cmin, idx, bits["r03"] ) @@ -444,7 +488,7 @@ def main( & (cmin <= 0.95) & (land_bits == 1) ) - cmin[idx], bits["r04"]["test"][idx] = Restore.land(land_bits, cmin, idx, bits["r04"]) + cmin[idx], bits["r04"]["test"][idx] = restore.land(land_bits, cmin, idx, bits["r04"]) # idx = np.nonzero( # (scene_flags["day"] == 1) @@ -455,7 +499,7 @@ def main( # & (cmin <= 0.95) # & (desert_bits == 1) # ) - # cmin[idx], bits["r04"]["test"][idx] = Restore.land(desert_bits, cmin, idx, bits["r04"]) + # cmin[idx], bits["r04"]["test"][idx] = restore.land(desert_bits, cmin, idx, bits["r04"]) coast_bits = restoral_bits["05"] # land_night_bits = restoral_bit["16"] @@ -467,12 +511,12 @@ def main( & (scene_flags["ice"] == 0) & (scene_flags["coast"] == 1) ) - cmin[idx], bits["r05"]["test"][idx] = Restore.coast(coast_bits, cmin, idx, bits["r05"]) + cmin[idx], bits["r05"]["test"][idx] = restore.coast(coast_bits, cmin, idx, bits["r05"]) # idx = np.nonzero((scene_flags['night'] == 1) & (scene_flags['land'] == 1) & # (scene_flags['snow'] == 0) & (scene_flags['ice'] == 0) & # (cmin <= 0.95)) - # cmin[idx], bits['r06']['test'][idx] = Restore.land_night(land_night_bits, + # cmin[idx], bits['r06']['test'][idx] = restore.land_night(land_night_bits, # cmin, # bits['r06'])[idx] @@ -486,22 +530,22 @@ def main( """ idx = np.nonzero((scene_flags['day'] == 1) & (scene_flags['water'] == 1) & (scene_flags['ice'] == 0)) - cmin[idx] = Restore.shallow_water(sh_water_bits, cmin)[idx] + cmin[idx] = restore.shallow_water(sh_water_bits, cmin)[idx] idx = np.nonzero((scene_flags['day'] == 1) & (scene_flags['land'] == 1) & (scene_flags['snow'] == 0) & (scene_flags['ice'] == 0) & (cmin <= 0.95)) - cmin[idx] = Restore.land(land_bits, cmin)[idx] + cmin[idx] = restore.land(land_bits, cmin)[idx] idx = np.nonzero((scene_flags['day'] == 1) & (scene_flags['land'] == 1) & (scene_flags['snow'] == 0) & (scene_flags['ice'] == 0) & (scene_flags['coast'] == 1)) - cmin[idx] = Restore.coast(coast_bits, cmin)[idx] + cmin[idx] = restore.coast(coast_bits, cmin)[idx] idx = np.nonzero((scene_flags['night'] == 1) & (scene_flags['land'] == 1) & (scene_flags['snow'] == 0) & (scene_flags['ice'] == 0) & (cmin <= 0.95)) - cmin[idx] = Restore.land_night(land_night_bits, cmin)[idx] + cmin[idx] = restore.land_night(land_night_bits, cmin)[idx] """ # bits translation MVCM-python -> MVCM-C # 01 13 test_11um @@ -729,10 +773,10 @@ if __name__ == "__main__": satellite = args.satellite or "snpp" sensor = args.sensor or "viirs" data_path = args.path or _datapath - mod02 = args.l1b or "" # _fname_mod02 - mod03 = args.geolocation or "" # _fname_mod03 - img02 = args.hires_l1b or None - img03 = args.hires_geo or None + mod02 = args.l1b + mod03 = args.geolocation + img02 = args.hires_l1b or "" + img03 = args.hires_geo or "" threshold_file = args.threshold or _threshold_file geos_atm_1 = args.atmos_1 or _geos_atm_1 geos_atm_2 = args.atmos_2 or _geos_atm_2 diff --git a/mvcm/read_data.py b/mvcm/read_data.py index 3885b8573c203561af17bd13fd3f7c54c0ab3259..004f975534e765f7faeeb6dba875ca003b738ace 100644 --- a/mvcm/read_data.py +++ b/mvcm/read_data.py @@ -3,12 +3,12 @@ import logging import os from datetime import datetime as dt -from typing import Dict import numpy as np import numpy.typing as npt import xarray as xr from attrs import Factory, define, field, validators +from rich.logging import RichHandler import ancillary_data as anc import mvcm.scene as scn @@ -86,15 +86,13 @@ _emissive_bands = [ "I05", ] +_LOG_FORMAT = "%(message)s" +logging.basicConfig(level="NOTSET", datefmt="[%X]", format=_LOG_FORMAT, handlers=[RichHandler()]) logger = logging.getLogger(__name__) -# 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') @define(kw_only=True, slots=True) -class CollectInputs(object): +class CollectInputs: """Class that collects all input files. Parameters @@ -190,6 +188,8 @@ class CollectInputs(object): validators.instance_of(str), ], ) + log_level: str = field(default="INFO") + logging.basicConfig(level=log_level) dims: tuple = field( default=("number_of_lines", "number_of_pixels"), @@ -654,7 +654,7 @@ class ReadAncillary(CollectInputs): logger.debug("Olson ecosystem file read successfully") return eco.reshape(self.out_shape) - def get_geos(self) -> Dict: + def get_geos(self) -> dict: """Read GEOS-5 data and interpolate the fields to the sensor resolution. Parameters @@ -663,7 +663,7 @@ class ReadAncillary(CollectInputs): Returns ------- - geos_data: Dict + geos_data: dict dictionary containing all quantities required by MVCM (see geos_variables here below) """ @@ -736,7 +736,7 @@ class ReadAncillary(CollectInputs): return ancillary_data -def get_data(satellite: str, sensor: str, file_names: Dict, hires: bool = False) -> xr.Dataset: +def get_data(satellite: str, sensor: str, file_names: dict, hires: bool = False) -> xr.Dataset: """Store every variable in one xarray dataset. Parameters @@ -745,7 +745,7 @@ def get_data(satellite: str, sensor: str, file_names: Dict, hires: bool = False) satellite name sensor: str sensor name - file_names: Dict + file_names: dict dictionary containing all input file names hires: bool use high resolution data @@ -832,7 +832,7 @@ def get_data(satellite: str, sensor: str, file_names: Dict, hires: bool = False) def get_data(satellite: str, sensor: str, - file_names: Dict[str, str], + file_names: dict[str, str], sunglint_angle: float, hires: bool = False) -> xr.Dataset: diff --git a/mvcm/restoral.py b/mvcm/restoral.py index ce42c414ac0f806686fa7ff72682945f408bb7c6..463d1543fd6eef331c7f0c17f846629d940fce90 100644 --- a/mvcm/restoral.py +++ b/mvcm/restoral.py @@ -1,7 +1,6 @@ """Restoral tests for MVCM.""" import logging -from typing import Dict, Tuple import numpy as np import xarray as xr @@ -72,7 +71,7 @@ logger = logging.getLogger(__name__) @define(kw_only=True, slots=True) -class Restoral(object): +class Restoral: """Initialize class.""" data: xr.Dataset = field( @@ -80,16 +79,17 @@ class Restoral(object): validators.instance_of(xr.Dataset), ] ) - thresholds: Dict = field( + thresholds: dict = field( validator=[ - validators.instance_of(Dict), + validators.instance_of(dict), ] ) - scene_flag: Dict = field( + scene_flag: dict = field( validator=[ - validators.instance_of(Dict), + validators.instance_of(dict), ] ) + log_level: str = field(default="INFO") # confidence: np.ndarray = field(validator=[validators.instance_of(np.ndarray), ]) # VIIRS does not have water vapor channels (6.7um and 7.2um) for this test @@ -99,8 +99,8 @@ class Restoral(object): # chk_spatial_var def spatial_variability( - self, confidence: np.ndarray, scn_idx: tuple, bits: Dict - ) -> Tuple[np.ndarray, np.ndarray]: + self, confidence: np.ndarray, scn_idx: tuple, bits: dict + ) -> tuple[np.ndarray, np.ndarray]: """Check spatial variability.""" threshold = self.thresholds["Daytime_Ocean_Spatial_Variability"] bt = self.data.M15.values @@ -144,8 +144,8 @@ class Restoral(object): # chk_sunglint def sunglint( - self, confidence: np.ndarray, scn_idx: tuple, bits: Dict - ) -> Tuple[np.ndarray, np.ndarray]: + self, confidence: np.ndarray, scn_idx: tuple, bits: dict + ) -> tuple[np.ndarray, np.ndarray]: """Restore Sun glint pixels.""" threshold = self.thresholds["Sun_Glint"] variability_threshold = self.thresholds["Daytime_Ocean_Spatial_Variability"]["var_11um"] @@ -194,8 +194,8 @@ class Restoral(object): # chk_shallow_water def shallow_water( - self, bits: Dict, confidence: np.ndarray, scn_idx: tuple, bits_test: Dict - ) -> Tuple[np.ndarray, np.ndarray]: + self, bits: dict, confidence: np.ndarray, scn_idx: tuple, bits_test: dict + ) -> tuple[np.ndarray, np.ndarray]: """Restore shallow water pixels.""" threshold = self.thresholds["Ocean_NDVI_and_Shallow_Water"] @@ -236,8 +236,8 @@ class Restoral(object): # chk_land def land( - self, bit: np.ndarray, confidence: np.ndarray, scn_idx: tuple, bits: Dict - ) -> Tuple[np.ndarray, np.ndarray]: + self, bit: np.ndarray, confidence: np.ndarray, scn_idx: tuple, bits: dict + ) -> tuple[np.ndarray, np.ndarray]: """Restore land pixels.""" rad = self.data.M15.values m05m04ratio = self.data.M08.values / self.data.M04.values @@ -294,8 +294,8 @@ class Restoral(object): # chk_coast def coast( - self, bit: np.ndarray, confidence: np.ndarray, scn_idx: tuple, bits: Dict - ) -> Tuple[np.ndarray, np.ndarray]: + self, bit: np.ndarray, confidence: np.ndarray, scn_idx: tuple, bits: dict + ) -> tuple[np.ndarray, np.ndarray]: """Restore coastal pixels.""" ndvi = utils.compute_ndvi(self.data.M05.values, self.data.M07.values) threshold = self.thresholds["Coastal_NDVI_Thresholds"] @@ -319,7 +319,7 @@ class Restoral(object): return confidence[scn_idx], bits["test"][scn_idx] # chk_land_night - def land_night(self, bit: np.ndarray, confidence: np.ndarray, bits: Dict) -> np.ndarray: + def land_night(self, bit: np.ndarray, confidence: np.ndarray, bits: dict) -> np.ndarray: """Restore land night pixels.""" rad = self.data.M15.values threshold = self.thresholds["Land_Restoral"] @@ -392,7 +392,7 @@ def local_nxn_standard_deviation(arr: np.ndarray, # eventually move this function to utility_functions def spatial_var(rad: np.ndarray, - threshold: Dict) -> np.ndarray: + threshold: dict) -> np.ndarray: """Compute 3x3 spatial variability and returns number of pixels within certain threshold""" test = sliding_window_view(np.pad(rad, [1, 1], mode='constant'), (3, 3)) -