diff --git a/mvcm/conf.py b/mvcm/conf.py index 797f278275465073d26462acc387eefcd1da0e79..d2b8676b4412833e83267ac0d582b65c9ea2203f 100644 --- a/mvcm/conf.py +++ b/mvcm/conf.py @@ -1,6 +1,11 @@ """Module defining confidence computations.""" + +import logging + import numpy as np +logger = logging.getLogger(__name__) + def conf_test_new(rad: np.ndarray, thr: np.ndarray) -> np.ndarray: """Compute confidence based on thresholds.""" @@ -220,7 +225,7 @@ def conf_test_dble(rad, coeffs): ### # Inner region passes test - print("I NEED TO REVIEW THIS TO WRITE IT MORE CLEARLY") + logger.warning("I NEED TO REVIEW THIS TO WRITE IT MORE CLEARLY") # FOR NOW ALPHA AND GAMMA ARE SWITCHED BECAUSE OF HOW THE ARRAYS ARE DEFINED. # THINK ON HOW THIS COULD BE WRITTEN SO THAT IT'S EASIER TO UNDERSTAND (AND DEBUG) # Value is within range of lower set of limits diff --git a/mvcm/main.py b/mvcm/main.py index 10dea73520feb57a03db09fbf5c588f04c0a3254..51c284a11d06dad48554131ea8877f5963c6e405 100644 --- a/mvcm/main.py +++ b/mvcm/main.py @@ -6,13 +6,11 @@ except ImportError: import ruamel.yaml as yml import argparse -import functools import logging -import time import numpy as np import xarray as xr -from netCDF4 import Dataset +from netCDF4 import Dataset # type: ignore from pkg_resources import get_distribution # type: ignore from rich.logging import RichHandler @@ -100,22 +98,6 @@ _mod_bands = [ _img_bands = ["I01", "I02", "I03", "I04", "I05"] -def timer(func): - """Compute elapsed time.""" - - @functools.wraps(func) - def wrapper_timer(*args, **kwargs): - start_time = time.time() - value = func(*args, **kwargs) - end_time = time.time() - run_time = end_time - start_time - logger.info(f"Finished {func.__name__!r} in {run_time:.2f} secs\n") - return value - - return wrapper_timer - - -@timer def main( satellite: str = "snpp", sensor: str = "viirs", @@ -182,8 +164,8 @@ def main( ------- None """ - verbose_level = np.minimum(verbose, 4) - logging.basicConfig(level=LOG_LEVELS[verbose_level]) + verbose_level = np.minimum(verbose + 1, 4) + logger.setLevel(LOG_LEVELS[verbose_level]) if img02 is None or img03 is None: use_hires = False @@ -231,9 +213,6 @@ def main( logger.info(f"All bands found in file {file_names['MOD02']}. The code will run.") with Dataset(file_names["MOD03"]) as f: - # time_coverage_start = f.getncattr("time_coverage_start") - # time_coverage_end = f.getncattr("time_coverage_end") - # orbit_number = f.getncattr("OrbitNumber") attrs = {attr: f.getncattr(attr) for attr in f.ncattrs()} viirs_data = rd.get_data(satellite, sensor, file_names, thresholds, hires=use_hires) @@ -241,12 +220,7 @@ def main( sunglint_angle = thresholds["Sun_Glint"]["bounds"][3] scene_flags = scn.find_scene(viirs_data, sunglint_angle, thresholds["Snow_Mask"]) - restore = Restoral( - data=viirs_data, - thresholds=thresholds, - scene_flag=scene_flags, - log_level=LOG_LEVELS[verbose_level], - ) + restore = Restoral(data=viirs_data, thresholds=thresholds, scene_flag=scene_flags) cmin_3d = np.ones((18, viirs_data.M11.shape[0], viirs_data.M11.shape[1]), dtype=np.float32) @@ -668,7 +642,9 @@ if __name__ == "__main__": version=_VERSION, help="print version and exit", ) - parser.add_argument("-v", "--verbose", action="store_true", help="print verbose information") + parser.add_argument( + "-v", "--verbose", action="count", default=1, help="print verbose information" + ) parser.add_argument("-d", "--debug", action="store_true", help="activate debug mode") args = parser.parse_args() diff --git a/mvcm/read_data.py b/mvcm/read_data.py index 80b634b9b9c86eb607d776f934797c346b9e6b8f..59e8ed6e415c240b11de499f9ec15724c9a84d5b 100644 --- a/mvcm/read_data.py +++ b/mvcm/read_data.py @@ -1,6 +1,6 @@ """Read l1b, geo and ancillary data.""" -import logging +import logging # noqa - pyright complains about imports not being sorted but pysort seem happy with the current sorting import os from datetime import datetime as dt @@ -8,7 +8,6 @@ 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 @@ -16,7 +15,7 @@ import mvcm.scene as scn _DTR = np.pi / 180.0 _RTD = 180.0 / np.pi _bad_data = -999.0 -_datapath = "/ships19/hercules/pveglio/mvcm_test_data" + _band_names = [ "um0_41", "um0_44", @@ -86,8 +85,6 @@ _emissive_bands = [ "I05", ] -_LOG_FORMAT = "%(message)s" -logging.basicConfig(level="NOTSET", datefmt="[%X]", format=_LOG_FORMAT, handlers=[RichHandler()]) logger = logging.getLogger(__name__) @@ -133,12 +130,9 @@ class CollectInputs: geos_land: str = "" geos_ocean: str = "" geos_constants: str = "GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4" - log_level: str = "INFO" dims: tuple = ("number_of_lines", "number_of_pixels") - logging.basicConfig(level=log_level) - @define(slots=True, kw_only=True) class ReadData(CollectInputs): @@ -178,8 +172,8 @@ class ReadData(CollectInputs): logger.debug(f"Reading {self.file_name_geo}") if os.path.exists(self.file_name_geo) is False: - logger.error("Geolocation file not found") err_msg = f"Could not find the file {self.file_name_geo}" + logger.error(err_msg) raise FileNotFoundError(err_msg) geo_data = xr.open_dataset(self.file_name_geo, group="geolocation_data", engine="netcdf4") @@ -215,8 +209,8 @@ class ReadData(CollectInputs): logger.debug(f"Reading {self.file_name_l1b}") if os.path.exists(self.file_name_l1b) is False: - logger.error("L1b file not found") err_msg = f"Could not find the file {self.file_name_l1b}" + logger.error(err_msg) raise FileNotFoundError(err_msg) l1b_data = xr.open_dataset( @@ -311,48 +305,6 @@ class ReadData(CollectInputs): else: gemi = np.full((viirs_out.M16.shape), _bad_data) - for band in _mod_bands[:11]: - if band in viirs: - if band == "M02": - pass - # idx = np.nonzero((viirs.M02.values > 1.3) & (viirs.M05 < 2.0)) - # viirs["M02"].values[idx] = 1.3 - # idx = np.nonzero((viirs.M02.values < -99) | (viirs.M02.values > 1.3)) - # viirs["M02"].values[idx] = _bad_data - else: - pass - # idx = np.nonzero((viirs[band].values < -99) | (viirs[band].values > 2.0)) - # viirs[band].values[idx] = _bad_data - for band in _mod_bands[11:]: - if band in viirs: - pass - # idx = np.nonzero((viirs[band].values < 0) | (viirs[band].values > 1000)) - # viirs[band].values[idx] = _bad_data - - """ - if 'M05' in viirs: - idx = np.nonzero((viirs.M05.values < -99) | (viirs.M05.values > 2)) - viirs['M05'].values[idx] = _bad_data - else: - viirs['M05'] = (self.dims, np.full(viirs.M15.shape, _bad_data)) - if 'M07' in viirs: - idx = np.nonzero((viirs.M07.values < -99) | (viirs.M07.values > 2)) - viirs['M07'].values[idx] = _bad_data - else: - viirs['M07'] = (self.dims, np.full(viirs.M15.shape, _bad_data)) - - idx = np.nonzero((viirs.M12.values < 0) | (viirs.M12.values > 1000)) - viirs['M12'].values[idx] = _bad_data - idx = np.nonzero((viirs.M13.values < 0) | (viirs.M13.values > 1000)) - viirs['M13'].values[idx] = _bad_data - idx = np.nonzero((viirs.M14.values < 0) | (viirs.M14.values > 1000)) - viirs['M14'].values[idx] = _bad_data - idx = np.nonzero((viirs.M15.values < 0) | (viirs.M15.values > 1000)) - viirs['M15'].values[idx] = _bad_data - idx = np.nonzero((viirs.M16.values < 0) | (viirs.M16.values > 1000)) - viirs['M16'].values[idx] = _bad_data - """ - # Compute channel differences and ratios that are used in the tests if ("M05" in viirs) and ("M07" in viirs): viirs_out["M07-M05ratio"] = ( @@ -384,10 +336,6 @@ class ReadData(CollectInputs): viirs_out["M10"] = (self.dims, hires_data.I03.values) viirs_out["M12"] = (self.dims, hires_data.I04.values) viirs_out["M15hi"] = (self.dims, hires_data.I05.values) - # b12 = viirs_out["M12"].values - # b12[b12 < 0] = _bad_data - # b12[b12 > 330] = _bad_data - # viirs_out["M12"].values = b12 # temp value to force the code to work viirs_out["M128"] = (self.dims, np.zeros(viirs_out.M16.shape)) @@ -524,7 +472,7 @@ class ReadAncillary(CollectInputs): def get_granule_time(self): """Get granule timestamp and format it for MVCM.""" vnp_time = ".".join(os.path.basename(self.file_name_l1b).split(".")[1:3]) - return dt.strftime(dt.strptime(vnp_time, "A%Y%j.%H%M"), "%Y-%m-%d %H:%M:%S.000") + return dt.strftime(dt.strptime(vnp_time, "A%Y%j.%H%M"), "%Y-%m-%d %H:%M:%S.000") # noqa def get_sst(self) -> npt.NDArray: """Read Reynolds SST file. @@ -540,7 +488,9 @@ class ReadAncillary(CollectInputs): """ logger.debug(f"Reading SST file {self.sst_file}") if not os.path.isfile(self.sst_file): - logger.error("SST file not found") + err_msg = f"SST file {self.sst_file} not found" + logger.error(err_msg) + raise FileNotFoundError(err_msg) sst = np.empty(self.out_shape, dtype=np.float32).ravel() sst = anc.py_get_Reynolds_SST( self.latitude.ravel(), @@ -568,7 +518,9 @@ class ReadAncillary(CollectInputs): """ logger.debug(f"Reading NDVI file {self.ndvi_file}") if not os.path.isfile(self.ndvi_file): - logger.error("NDVI file not found") + err_msg = f"NDVI file {self.ndvi_file} not found" + logger.error(err_msg) + raise FileNotFoundError(err_msg) ndvi = np.empty(self.out_shape, dtype=np.float32).ravel() ndvi = anc.py_get_NDVI_background( self.latitude.ravel(), @@ -596,7 +548,9 @@ class ReadAncillary(CollectInputs): """ logger.debug(f"Reading Ecosystem file {self.eco_file}") if not os.path.isfile(self.eco_file): - logger.error("Ecosystem file not found") + err_msg = f"Ecosystem file {self.eco_file} not found" + logger.error(err_msg) + raise FileNotFoundError(err_msg) eco = np.empty(self.out_shape, dtype=np.ubyte).ravel() eco = anc.py_get_Olson_eco( @@ -630,19 +584,27 @@ class ReadAncillary(CollectInputs): logger.debug(f"Reading GEOS-5 land file {self.geos_land}") logger.debug(f"Reading GEOS-5 ocean file {self.geos_ocean}") logger.debug(f"Reading GEOS-5 constants file {self.geos_constants}") + if not os.path.isfile(self.geos_file_1): + err_msg = f"GEOS-5 file 1 {self.geos_file_1} not found" logger.error("GEOS-5 file 1 not found") - - # if not os.path.isfile(os.path.join(self.ancillary_dir, self.geos_file_1)): - # logger.error("GEOS-5 file 1 not found") - # if not os.path.isfile(os.path.join(self.ancillary_dir, self.geos_file_2)): - # logger.error("GEOS-5 file 2 not found") - # if not os.path.isfile(os.path.join(self.ancillary_dir, self.geos_land)): - # logger.error("GEOS-5 land file not found") - # if not os.path.isfile(os.path.join(self.ancillary_dir, self.geos_ocean)): - # logger.error("GEOS-5 ocean file not found") - # if not os.path.isfile(os.path.join(self.ancillary_dir, self.geos_constants)): - # logger.error("GEOS-5 constants file not found") + raise FileNotFoundError(err_msg) + if not os.path.isfile(self.geos_file_2): + err_msg = f"GEOS-5 file 2 {self.geos_file_2} not found" + logger.error("GEOS-5 file 2 not found") + raise FileNotFoundError(err_msg) + if not os.path.isfile(self.geos_land): + err_msg = f"GEOS-5 land file {self.geos_land} not found" + logger.error("GEOS-5 land file not found") + raise FileNotFoundError(err_msg) + if not os.path.isfile(self.geos_ocean): + err_msg = f"GEOS-5 ocean file {self.geos_ocean} not found" + logger.error("GEOS-5 ocean file not found") + raise FileNotFoundError(err_msg) + if not os.path.isfile(self.geos_constants): + err_msg = f"GEOS-5 constants file {self.geos_constants} not found" + logger.error("GEOS-5 constants file not found") + raise FileNotFoundError(err_msg) geos_variables = [ "tpw", @@ -702,7 +664,11 @@ class ReadAncillary(CollectInputs): def get_data( - satellite: str, sensor: str, file_names: dict, thresholds: dict, hires: bool = False + satellite: str, + sensor: str, + file_names: dict, + thresholds: dict, + hires: bool = False, ) -> xr.Dataset: """Store every variable in one xarray dataset. @@ -722,14 +688,15 @@ def get_data( input_data: xr.Dataset dataset containing all the satellite and ancillary data """ - Viirs = ReadData( + viirs = ReadData( file_name_geo=file_names["MOD03"], file_name_l1b=file_names["MOD02"], satellite=satellite, sensor=sensor, ) + if hires is True: - Viirs_hires = ReadData( + viirs_hires = ReadData( file_name_geo=file_names["IMG03"], file_name_l1b=file_names["IMG02"], satellite=satellite, @@ -737,19 +704,19 @@ def get_data( hires=True, ) - geo_data_mod = Viirs.read_viirs_geo() - viirs_data_mod = Viirs.read_viirs_l1b(geo_data_mod.solar_zenith.values) + geo_data_mod = viirs.read_viirs_geo() + viirs_data_mod = viirs.read_viirs_l1b(geo_data_mod.solar_zenith.values) if hires is True: - geo_data = Viirs_hires.read_viirs_geo() - viirs_data_img = Viirs_hires.read_viirs_l1b(geo_data.solar_zenith.values) - viirs_data = Viirs_hires.preprocess_viirs(geo_data, viirs_data_mod, viirs_data_img) + geo_data = viirs_hires.read_viirs_geo() + viirs_data_img = viirs_hires.read_viirs_l1b(geo_data.solar_zenith.values) + viirs_data = viirs_hires.preprocess_viirs(geo_data, viirs_data_mod, viirs_data_img) res = 1 else: - viirs_data = Viirs.preprocess_viirs(geo_data_mod, viirs_data_mod) + viirs_data = viirs.preprocess_viirs(geo_data_mod, viirs_data_mod) geo_data = geo_data_mod res = 1 - Ancillary = ReadAncillary( + ancillary = ReadAncillary( file_name_l1b=file_names["MOD02"], latitude=geo_data.latitude.values, longitude=geo_data.longitude.values, @@ -765,7 +732,7 @@ def get_data( eco_file=file_names["ECO"], ) - viirs_data.update(Ancillary.pack_data()) + viirs_data.update(ancillary.pack_data()) scene_flags = scn.find_scene( viirs_data, geo_data.sunglint_angle.values, thresholds["Snow_Mask"] @@ -780,136 +747,3 @@ def get_data( input_data.drop_vars(["latitude", "longitude"]) return input_data - - -""" - - 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']) - - - -def get_data(satellite: str, - sensor: str, - file_names: dict[str, str], - sunglint_angle: float, - hires: bool = False) -> xr.Dataset: - - mod02 = file_names['MOD02'] - mod03 = file_names['MOD03'] - - if hires is True: - img02 = file_names['IMG02'] - img03 = file_names['IMG03'] - - if hires is False: - viirs = read_viirs(sensor, f'{mod02}', f'{mod03}') - viirs = read_ancillary_data(file_names, viirs) - - if (('M05' in viirs) and ('M07' in viirs)): - m01 = viirs.M05.values - m02 = viirs.M07.values - r1 = 2.0 * (np.power(m02, 2.0) - - np.power(m01, 2.0)) + (1.5 * m02) + (0.5 * m01) - r2 = m02 + m01 + 0.5 - r3 = r1 / r2 - gemi = r3 * (1.0 - 0.25*r3) - ((m01 - 0.125) / (1.0 - m01)) - else: - gemi = np.full((viirs.M15.shape), _bad_data) - - if 'M05' in viirs: - idx = np.nonzero((viirs.M05.values < -99) | (viirs.M05.values > 2)) - viirs['M05'].values[idx] = _bad_data - else: - viirs['M05'] = (('number_of_lines', 'number_of_pixels'), - np.full(viirs.M15.shape, _bad_data)) - if 'M07' in viirs: - idx = np.nonzero((viirs.M07.values < -99) | (viirs.M07.values > 2)) - viirs['M07'].values[idx] = _bad_data - else: - viirs['M07'] = (('number_of_lines', 'number_of_pixels'), - np.full(viirs.M15.shape, _bad_data)) - - idx = np.nonzero((viirs.M12.values < 0) | (viirs.M12.values > 1000)) - viirs['M12'].values[idx] = _bad_data - idx = np.nonzero((viirs.M13.values < 0) | (viirs.M13.values > 1000)) - viirs['M13'].values[idx] = _bad_data - idx = np.nonzero((viirs.M14.values < 0) | (viirs.M14.values > 1000)) - viirs['M14'].values[idx] = _bad_data - idx = np.nonzero((viirs.M15.values < 0) | (viirs.M15.values > 1000)) - viirs['M15'].values[idx] = _bad_data - idx = np.nonzero((viirs.M16.values < 0) | (viirs.M16.values > 1000)) - viirs['M16'].values[idx] = _bad_data - - # Compute channel differences and ratios that are used in the tests - viirs['M15-M13'] = (('number_of_lines', 'number_of_pixels'), - viirs.M15.values - viirs.M13.values) - viirs['M14-M15'] = (('number_of_lines', 'number_of_pixels'), - viirs.M14.values - viirs.M15.values) - viirs['M15-M16'] = (('number_of_lines', 'number_of_pixels'), - viirs.M15.values - viirs.M16.values) - viirs['M15-M12'] = (('number_of_lines', 'number_of_pixels'), - viirs.M15.values - viirs.M12.values) - viirs['M13-M16'] = (('number_of_lines', 'number_of_pixels'), - viirs.M13.values - viirs.M16.values) - viirs['M07-M05ratio'] = (('number_of_lines', 'number_of_pixels'), - viirs.M07.values / viirs.M05.values) - viirs['GEMI'] = (('number_of_lines', 'number_of_pixels'), gemi) - - # temp value to force the code to work - viirs['M128'] = (('number_of_lines', 'number_of_pixels'), - np.zeros(viirs.M15.shape)) - - else: - viirs = read_viirs('viirs', f'{img02}', f'{img03}') - viirs['M05'] = viirs.I01 - viirs['M07'] = viirs.I02 - - idx = np.nonzero((viirs.M05.values < -99) | (viirs.M05.values > 2)) - viirs['M05'].values[idx] = _bad_data - idx = np.nonzero((viirs.M07.values < -99) | (viirs.M07.values > 2)) - viirs['M07'].values[idx] = _bad_data - - idx = np.nonzero((viirs.I01.values < 0) | (viirs.I01.values > 1000)) - viirs['I01'].values[idx] = _bad_data - idx = np.nonzero((viirs.I02.values < 0) | (viirs.I02.values > 1000)) - viirs['I02'].values[idx] = _bad_data - idx = np.nonzero((viirs.I03.values < 0) | (viirs.I03.values > 1000)) - viirs['I03'].values[idx] = _bad_data - idx = np.nonzero((viirs.I04.values < 0) | (viirs.I04.values > 1000)) - viirs['I04'].values[idx] = _bad_data - idx = np.nonzero((viirs.I05.values < 0) | (viirs.I05.values > 1000)) - viirs['I05'].values[idx] = _bad_data - - viirs = read_ancillary_data(file_names, viirs, resolution=2) - - viirs['I05-I04'] = (('number_of_lines_2', 'number_of_pixels_2'), - viirs.I05.values - viirs.I04.values) - viirs['I02-I01ratio'] = (('number_of_lines_2', 'number_of_pixels_2'), - viirs.I02.values / viirs.I01.values) - - 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 -""" diff --git a/mvcm/restoral.py b/mvcm/restoral.py index 463d1543fd6eef331c7f0c17f846629d940fce90..95288b642c17cf6f2ac4dbbf837fe65ca682dcc0 100644 --- a/mvcm/restoral.py +++ b/mvcm/restoral.py @@ -12,63 +12,6 @@ _bad_data = -999.0 logger = logging.getLogger(__name__) -# tracking down the list of restoral tests that need to be implemented -# and the function that gets called in the C code -# Also list the set of conditions to run the restoral test -# -# Night, snow/ice temperature inversion tests -# - NightSnow_Inv_check() -# - pxin.night -# - pxin.snow || pxin.ice -# -# Water -# - chk_spatial_var() -# - pxin.water -# - pxin.uniform -# - !pxin.ice -# - pxout.intermediate_conf <= 0.99 -# - pxout.intermediate_conf >= 0.05 -# -# Water Day -# - chck_sunglint() -# - pxin.day -# - pxin.water -# - pxin.uniform -# - pxin.sunglint -# - pxout.intermediate_conf <= 0.95 -# -# Daytime water NDVI test and turbid water test for shallow water -# - chk_shallow_water() -# - pxin.day -# - pxin.water -# - !pxin.ice -# -# Daytime land -# - chk_land() -# - pxin.day -# - pxin.land -# - !pxin.snow -# - !pxin.ice -# - pxout.intermediate_conf <= 0.95 - -# -# Daytime coast(land) -# - chk_land() -# - pxin.day -# - pxin.land -# - !pxin.snow -# - !pxin.ice -# - pxin.coast -# -# Nighttime land test -# - chck_land_night() -# - pxin.night -# - pxin.land -# - !pxin.snow -# - !pxin.ice -# - pxout.intermediate_conf <= 0.95 -# - @define(kw_only=True, slots=True) class Restoral: @@ -89,7 +32,6 @@ class Restoral: 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 diff --git a/mvcm/write_output.py b/mvcm/write_output.py index 7934d1af11d0d0144066881928f21eba93ca9a68..98aebe93c9277d0354e71d8aa09ffbf84d3796fe 100644 --- a/mvcm/write_output.py +++ b/mvcm/write_output.py @@ -223,7 +223,8 @@ def save_output(data: xr.Dataset, attrs: dict, fname: str, compression: int) -> out_ds = xr.Dataset() for attr in attrs: - out_ds.attrs[attr] = attrs[attr] + if attr in _attributes_list: + out_ds.attrs[attr] = attrs[attr] for attr in attributes: out_ds.attrs[attr] = attributes[attr] logger.debug("attributes added to datasets")