diff --git a/mvcm/main.py b/mvcm/main.py
index c0832d67be72e81b711b542bc851d5c0ac1e4ee3..cf59473ba41521b0e2c3ab88e195a1aafec47ccf 100644
--- a/mvcm/main.py
+++ b/mvcm/main.py
@@ -168,19 +168,7 @@ def main(
try:
vnp02[b]
except KeyError:
- logger.info(f"Band {b} not found in file. No output will be written.")
- with Dataset(file_names["MOD03"]) as f:
- attrs = {attr: f.getncattr(attr) for attr in f.ncattrs()}
- shape = f["geolocation_data/latitude"][:].shape
- save_output(
- None,
- attrs,
- out_file,
- compression=5,
- debug=debug,
- shape=shape,
- )
- return
+ logger.info(f"Band {b} not found in file. Fill values will be used instead.")
logger.info(f"All bands found in file {file_names['MOD02']}. The code will run.")
with Dataset(file_names["MOD03"]) as f:
diff --git a/mvcm/main_prep_only.py b/mvcm/main_prep_only.py
new file mode 100644
index 0000000000000000000000000000000000000000..eb24c3cf174b5ae64a236b0393089666c77149b4
--- /dev/null
+++ b/mvcm/main_prep_only.py
@@ -0,0 +1,295 @@
+"""Main function for MVCM."""
+
+import argparse
+import logging
+import os
+from datetime import datetime as dt
+
+import numpy as np
+import psutil
+import xarray as xr
+from netCDF4 import Dataset # type: ignore
+from pkg_resources import get_distribution # type: ignore
+from rich.logging import RichHandler
+from ruamel.yaml import YAML
+
+import mvcm.read_data as rd
+from mvcm.constants import SensorConstants
+from mvcm.write_output import run_mvcm, save_output
+
+_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"]
+
+proc = psutil.Process(os.getpid())
+
+
+def main(
+ satellite: str = "snpp",
+ sensor: str = "viirs",
+ data_path: str = "",
+ mod02: str = "",
+ mod03: str = "",
+ img02: str | None = "",
+ img03: str | None = "",
+ threshold_file: str = "",
+ geos_atm_1: str = "",
+ geos_atm_2: str = "",
+ geos_land: str = "",
+ geos_ocean: str = "",
+ geos_constants: str = "",
+ ndvi_file: str = "",
+ sst_file: str = "",
+ eco_file: str = "",
+ out_file: str = "",
+ exe_path: str = "/home/pveglio/mvcm/mvcm",
+ verbose: int = 0,
+ *,
+ debug: bool = False,
+ hires_only: bool = False,
+ mvcm_threshold_file: str | None,
+) -> None:
+ """MVCM main function.
+
+ 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: [optional] str
+ VIIRS IMG02 file name
+ img03: [optional] 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: str
+ Sea surface temperature file name
+ eco_file: str
+ Ecosystem File
+ out_file: str
+ Output file name
+ verbose: int
+ Verbosity level
+ exe_path: str
+ Path to the MVCM executable
+ debug: bool
+ Debug mode
+ hires_only: bool
+ Only run MVCM with high resolution
+
+ Returns
+ -------
+ None
+ """
+ verbose_level = np.minimum(verbose + 1, 4)
+ logger.setLevel(LOG_LEVELS[verbose_level])
+
+ if img02 is None or img03 is None:
+ use_hires = False
+ else:
+ use_hires = True
+
+ file_names = {
+ "MOD02": f"{mod02}",
+ "MOD03": f"{mod03}",
+ "IMG02": f"{img02}",
+ "IMG03": f"{img03}",
+ "GEOS_atm_1": f"{geos_atm_1}",
+ "GEOS_atm_2": f"{geos_atm_2}",
+ "GEOS_land": f"{geos_land}",
+ "GEOS_ocean": f"{geos_ocean}",
+ "GEOS_constants": f"{geos_constants}",
+ "NDVI": f"{ndvi_file}",
+ "SST": f"{sst_file}",
+ "ECO": f"{eco_file}",
+ "ANC_DIR": f"{data_path}",
+ }
+ logger.info(f"Memory usage #1: {proc.memory_info().rss / 1e6} MB")
+
+ if hires_only is False:
+ logger.info("Running regular MVCM before high resolution")
+ run_mvcm(file_names, mvcm_threshold_file, out_file, exe_path)
+
+ with open(threshold_file) as f:
+ text = f.read()
+ thresholds = YAML(typ="safe").load(text)
+
+ viirs_m_bands = SensorConstants.VIIRS_VIS_BANDS + SensorConstants.VIIRS_IR_BANDS
+
+ with Dataset(file_names["MOD03"]) as f:
+ attrs = {attr: f.getncattr(attr) for attr in f.ncattrs()}
+ lat_shape = f["geolocation_data/latitude"][:].shape
+
+ fmt_in = "%Y-%m-%dT%H:%M:%S.%fZ"
+ fmt_out = "%Y-%m-%dT%H%M"
+
+ # 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 SensorConstants.VIIRS_IMG_BANDS:
+ try:
+ vnp02[b]
+ except KeyError:
+ logger.info(f"Band {b} not found in file. Hires data will be fill values.")
+ attrs = {"version": _VERSION, "satellite": "VIIRS", "sensor": "VIIRS"}
+ save_output(None, attrs, out_file, compression=5, debug=debug, shape=lat_shape)
+ return
+ logger.info(f"All bands found in file {file_names['IMG02']}. The code will run.")
+ else:
+ with xr.open_dataset(file_names["MOD02"], group="observation_data") as vnp02:
+ for b in viirs_m_bands:
+ try:
+ vnp02[b]
+ except KeyError:
+ logger.info(f"Band {b} not found in file. Hires data will be fill values.")
+ attrs = {"version": _VERSION, "satellite": "VIIRS", "sensor": "VIIRS"}
+ save_output(None, attrs, out_file, compression=5, debug=debug, shape=lat_shape)
+ return
+ logger.info(f"All bands found in file {file_names['MOD02']}. The code will run.")
+
+ logger.info(f"Memory usage #2: {proc.memory_info().rss / 1e6} MB")
+ # rd.get_data(satellite, sensor, file_names, thresholds, hires=use_hires)
+ viirs_data, pixel_type = rd.get_data(satellite, sensor, file_names, thresholds, hires=use_hires)
+ # viirs_data = rd.get_data(satellite, sensor, file_names, thresholds, hires=use_hires)
+
+ # viirs_data = xr.open_dataset("input_data.nc")
+ logger.info(f"Memory usage #3: {proc.memory_info().rss / 1e6} MB")
+
+ with xr.open_dataset(file_names["MOD02"]) as vnp02:
+ time_str = dt.strftime(dt.strptime(vnp02.time_coverage_start, fmt_in), fmt_out)
+
+ viirs_data.to_netcdf(
+ f"{os.path.dirname(out_file)}/input_data_{time_str}.nc", mode="w", format="NETCDF4"
+ )
+ pixel_type.to_netcdf(
+ f"{os.path.dirname(out_file)}/pixel_type_{time_str}.nc", mode="w", format="NETCDF4"
+ )
+
+ logger.info(f"Memory usage #4: {proc.memory_info().rss / 1e6} MB")
+ ###
+
+
+if __name__ == "__main__":
+ _VERSION = get_distribution("mvcm").version
+ parser = argparse.ArgumentParser(prog="MVCM", description="")
+
+ parser.add_argument(
+ "--satellite",
+ help="satellite name, not case-sensitive. Available options: [snpp, ]",
+ choices=[
+ "snpp",
+ ],
+ )
+ parser.add_argument(
+ "--sensor",
+ help="sensor name, not case-sensitive. Available options: [viirs, ]",
+ choices=[
+ "viirs",
+ ],
+ )
+ parser.add_argument("--path", help="path where the data is stored")
+ parser.add_argument("--l1b", help="level 1b file name")
+ parser.add_argument("--geolocation", help="geolocation file name")
+ parser.add_argument("--hires-l1b", help="VIIRS IMG02 file name")
+ parser.add_argument("--hires-geo", help="VIIRS IMG03 file name")
+ parser.add_argument("-t", "--threshold", help="thresholds file name")
+ parser.add_argument(
+ "--atmos-1",
+ help="file name of the first GEOS-IT file for atmospheric parameters",
+ )
+ parser.add_argument(
+ "--atmos-2",
+ help="file name of the second GEOS-IT file for atmospheric parameters",
+ )
+ parser.add_argument("--land", help="file name of GEOS-IT land parameters")
+ parser.add_argument("--ocean", help="file name of GEOS-IT ocean parameters")
+ parser.add_argument("--constants", help="file name of GEOS-IT constants")
+ parser.add_argument("--ndvi", help="NDVI file name")
+ parser.add_argument("--sst", help="Sea surface temperature file name")
+ parser.add_argument("--eco", help="Ecosystem file")
+ parser.add_argument("-o", "--out", help="output file name")
+ parser.add_argument(
+ "-V",
+ "--version",
+ action="version",
+ version=_VERSION,
+ help="print version and exit",
+ )
+ 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")
+ parser.add_argument("--hires-only", action="store_true", help="run only high resolution code")
+ parser.add_argument("-x", "--mvcm-exe-path", help="path to mvcm executable")
+ parser.add_argument(
+ "--mvcm-thresholds", help="thresholds file name for operational MVCM", default=None
+ )
+
+ args = parser.parse_args()
+
+ satellite = args.satellite or "snpp"
+ sensor = args.sensor or "viirs"
+ data_path = args.path
+ mod02 = args.l1b
+ mod03 = args.geolocation
+ img02 = args.hires_l1b or None
+ img03 = args.hires_geo or None
+ threshold_file = args.threshold
+ geos_atm_1 = args.atmos_1
+ geos_atm_2 = args.atmos_2
+ geos_land = args.land
+ geos_ocean = args.ocean
+ constants = args.constants
+ ndvi_file = args.ndvi
+ sst_file = args.sst
+ eco_file = args.eco
+ out_file = args.out
+ verbose = args.verbose or False
+ debug = args.debug or False
+
+ main(
+ satellite=satellite,
+ sensor=sensor,
+ data_path=data_path,
+ mod02=mod02,
+ mod03=mod03,
+ img02=img02,
+ img03=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=constants,
+ ndvi_file=ndvi_file,
+ sst_file=sst_file,
+ eco_file=eco_file,
+ out_file=out_file,
+ verbose=verbose,
+ debug=debug,
+ exe_path=args.mvcm_exe_path,
+ hires_only=args.hires_only,
+ mvcm_threshold_file=args.mvcm_thresholds,
+ )
+
+# tracemalloc.stop()
diff --git a/mvcm/main_tests_only.py b/mvcm/main_tests_only.py
new file mode 100644
index 0000000000000000000000000000000000000000..f68660e621c97cf4fb062732f9c8bcfbc2a9ab4c
--- /dev/null
+++ b/mvcm/main_tests_only.py
@@ -0,0 +1,612 @@
+"""Main function for MVCM."""
+
+import argparse
+import gc
+import logging
+import os
+
+import numpy as np
+import psutil
+import xarray as xr
+from pkg_resources import get_distribution # type: ignore
+from rich.logging import RichHandler
+from ruamel.yaml import YAML
+
+import mvcm.spectral_tests as tst
+import mvcm.utility_functions as utils
+from mvcm.constants import SceneConstants, SensorConstants
+from mvcm.restoral import Restoral
+from mvcm.write_output import save_output
+
+# import mvcm.temp_spectral_tests as tst
+
+
+_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"]
+
+proc = psutil.Process(os.getpid())
+
+
+def main(
+ threshold_file: str = "",
+ pixel_data: str = "",
+ input_data: str = "",
+ out_file: str = "",
+ verbose: int = 0,
+ *,
+ debug: bool = False,
+ use_hires: bool = False,
+) -> None:
+ """MVCM main function.
+
+ 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: [optional] str
+ VIIRS IMG02 file name
+ img03: [optional] 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: str
+ Sea surface temperature file name
+ eco_file: str
+ Ecosystem File
+ out_file: str
+ Output file name
+ verbose: int
+ Verbosity level
+ exe_path: str
+ Path to the MVCM executable
+ debug: bool
+ Debug mode
+ hires_only: bool
+ Only run MVCM with high resolution
+
+ Returns
+ -------
+ None
+ """
+ verbose = 4
+ verbose_level = np.minimum(verbose + 1, 4)
+ logger.setLevel(LOG_LEVELS[verbose_level])
+
+ with open(threshold_file) as f:
+ text = f.read()
+ thresholds = YAML(typ="safe").load(text)
+
+ pixel_type = xr.open_dataset(pixel_data, chunks="auto")
+ viirs_data = xr.open_dataset(input_data, chunks="auto")
+
+ # cmin_3d = np.ones((18, viirs_data.M11.shape[0], viirs_data.M11.shape[1]), dtype=np.float32)
+ cmin = np.ones(viirs_data.latitude.shape, dtype=np.float32)
+ # cmin_3d = xr.ones_like(viirs_data.latitude, dtype=np.float32, chunks="auto")
+
+ logger.info(f"Memory usage #4: {proc.memory_info().rss / 1e6} MB")
+
+ ##########################################################
+ _bitlist = [
+ "01",
+ "02",
+ "03",
+ "04",
+ "05",
+ "06",
+ "07",
+ "08",
+ "09",
+ "10",
+ "11",
+ "12",
+ "13",
+ "14",
+ "15",
+ "16",
+ "r01",
+ "r02",
+ "r03",
+ "r04",
+ "r05",
+ "r06",
+ ]
+ bits = {}
+
+ for b in _bitlist:
+ bits[b] = {
+ "test": np.zeros(viirs_data.latitude.shape, dtype=np.ubyte),
+ "qa": np.zeros(viirs_data.latitude.shape, dtype=np.ubyte),
+ # "test": xr.zeros_like(viirs_data.latitude, dtype=np.ubyte, chunks="auto"),
+ # "qa": xr.zeros_like(viirs_data.latitude, dtype=np.ubyte, chunks="auto"),
+ }
+ # i = 0
+
+ # bits_xr = xr.Dataset(
+ # {
+ # "qa": np.zeros(viirs_data.latitude.shape, dtype=np.ubyte),
+ # "test": np.zeros(viirs_data.latitude.shape, dtype=np.ubyte),
+ # },
+ # # coords={"latitude": viirs_data.latitude, "longitude": viirs_data.longitude},
+ # )
+ viirs_data["locut"] = np.full_like(
+ viirs_data.latitude.shape, fill_value=np.nan, dtype=np.float32
+ )
+ viirs_data["midpt"] = np.full_like(
+ viirs_data.latitude.shape, fill_value=np.nan, dtype=np.float32
+ )
+ viirs_data["hicut"] = np.full_like(
+ viirs_data.latitude.shape, fill_value=np.nan, dtype=np.float32
+ )
+ # logger.info(f"viirs_data: {viirs_data}")
+ logger.info(f"Memory usage #5: {proc.memory_info().rss / 1e6} MB")
+ # scene_types = np.zeros(viirs_data.M11.shape, dtype=np.ubyte)
+ for scene_name in SceneConstants.SCENE_LIST:
+ # for scene_name in ["Land_Day"]:
+ # scene_types[viirs_data[scene_name].values == 1] = i
+
+ logger.info(f"Processing {scene_name}")
+
+ if np.all(viirs_data[scene_name].values == 0):
+ # if viirs_data[scene_name].all() is False:
+ logger.info("Skipping, no pixels in scene.")
+ continue
+
+ logger.debug("initializing CloudTests class")
+ my_scene = tst.CloudTests(
+ data=viirs_data, scene_name=scene_name, thresholds=thresholds, hires=use_hires
+ )
+
+ # Initialize the confidence arrays for the various test groups
+ # cmin_g1 = xr.ones_like(viirs_data.latitude, dtype=np.float32, chunks="auto")
+ # cmin_g2 = xr.ones_like(viirs_data.latitude, dtype=np.float32, chunks="auto")
+ # cmin_g3 = xr.ones_like(viirs_data.latitude, dtype=np.float32, chunks="auto")
+ # cmin_g4 = xr.ones_like(viirs_data.latitude, dtype=np.float32, chunks="auto")
+ # cmin_g5 = xr.ones_like(viirs_data.latitude, dtype=np.float32, chunks="auto")
+
+ cmin_g1 = np.ones(viirs_data.latitude.shape, dtype=np.float32)
+ cmin_g2 = np.ones(viirs_data.latitude.shape, dtype=np.float32)
+ cmin_g3 = np.ones(viirs_data.latitude.shape, dtype=np.float32)
+ cmin_g4 = np.ones(viirs_data.latitude.shape, dtype=np.float32)
+ cmin_g5 = np.ones(viirs_data.latitude.shape, dtype=np.float32)
+ # cmin_temp = np.ones(viirs_data.M11.shape)
+ # cmin_xr = xr.ones_like(viirs_data.latitude, dtype=np.float32, chunks="auto")
+ # cmin_xr2 = xr.ones_like(viirs_data.latitude, dtype=np.float32, chunks="auto")
+ # cmin_xr4 = xr.ones_like(viirs_data.latitude, dtype=np.float32, chunks="auto")
+ logger.debug("starting tests")
+
+ if use_hires is True:
+ m15_name = "M15hi"
+ else:
+ m15_name = "M15"
+
+ logger.info(f"Memory usage #6: {proc.memory_info().rss / 1e6} MB")
+ # Group 1
+ 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"] = my_scene.surface_temperature_test("", viirs_data, cmin_g1, bits["02"])
+ cmin_g1, bits["03"] = my_scene.sst_test("M15", "M16", cmin_g1, bits["03"])
+
+ # Group 2
+ 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_xr2, bits_xr = my_scene.xr_test_11_12um_diff(cmin_xr2)
+ # cmin_g2 = np.minimum(cmin_g2, cmin_xr2.values)
+ # bits["05"]["qa"] = bits_xr.qa.values
+ # bits["05"]["test"] = bits_xr.test.values
+ 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"] = my_scene.bt_difference_11_4um_test_land(
+ # "M15-M16", "M15-M12", cmin_g2, bits["08"]
+ # )
+ cmin_g2, bits["09"] = my_scene.oceanic_stratus_11_4um_test("M15-M12", cmin_g2, bits["09"])
+
+ # Group 3
+ 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"] = 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"] = 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"] = my_scene.test_1_38um_high_clouds("M09", cmin_g4, bits["15"])
+ # cmin_xr4, bits_xr = my_scene.xr_test_1_38um_high_clouds(cmin_xr4)
+ # bits["15"]["qa"] = bits_xr.qa.values
+ # bits["15"]["test"] = bits_xr.test.values
+ # cmin_g4 = np.minimum(cmin_g4, cmin_xr4.values)
+
+ # Group 5
+ cmin_g5, bits["16"] = my_scene.thin_cirrus_4_12um_BTD_test("M13-M16", cmin_g5, bits["16"])
+
+ # logger.debug(f"Memory: {tracemalloc.get_traced_memory()}")
+ bit = {}
+ qabit = {}
+ # restoral_bits = {}
+ for b in _bitlist:
+ bit[b] = bits[f"{b}"]["test"]
+ qabit[b] = bits[f"{b}"]["qa"]
+ # restoral_bits[b] = utils.restoral_flag(bits[f"{b}"])
+ # # if utils.group_count(bit) != 0:
+ # cmin_g4 = cmin_xr.values
+
+ # cmin_temp = np.power(
+ # cmin_g1 * cmin_g2 * cmin_g3 * cmin_g4 * cmin_g5,
+ # 1 / utils.group_count(qabit),
+ # )
+
+ cmin = np.minimum(
+ np.power(cmin_g1 * cmin_g2 * cmin_g3 * cmin_g4 * cmin_g5, 1 / utils.group_count(qabit)),
+ cmin,
+ )
+ # cmin_3d[i, :, :] = cmin
+ # cmin_3d = np.minimum(cmin_3d, cmin)
+ # cmin = np.minimum(cmin_temp, cmin)
+ # i += 1
+
+ logger.info(f"Memory usage #7: {proc.memory_info().rss / 1e6} MB")
+
+ # cmin = np.min(cmin_3d, axis=0)
+ # cmin = cmin_3d
+ # pixel_type = xr.open_dataset("pixels_data.nc")
+ restore = Restoral(data=viirs_data, thresholds=thresholds, scene_flags=pixel_type)
+ logger.debug("Instance of Restoral class created successfully.")
+ logger.info(f"Memory usage #8: {proc.memory_info().rss / 1e6} MB")
+
+ # logger.debug(f"Memory: {tracemalloc.get_traced_memory()}")
+ # bit = {}
+ restoral_bits = {}
+ for b in _bitlist:
+ # bit[b] = bits[f"{b}"]["test"]
+ restoral_bits[b] = utils.restoral_flag(bits[f"{b}"])
+ # if utils.group_count(bit) != 0:
+
+ # Restoral Tests
+
+ # Spatial Variability
+ idx = np.nonzero(
+ (pixel_type["uniform"].values == 1)
+ & (pixel_type["water"].values == 1)
+ & (pixel_type["ice"].values == 0)
+ & (cmin >= 0.05)
+ & (cmin <= 0.99)
+ )
+ cmin[idx], bits["r01"]["test"][idx] = restore.spatial_variability(cmin, idx, bits["r01"])
+
+ # logger.debug(f"Memory: {tracemalloc.get_traced_memory()}")
+ # Sun Glint
+ sunglint_bits = (
+ restoral_bits["01"] * restoral_bits["03"] * restoral_bits["05"] * restoral_bits["15"]
+ )
+
+ idx = np.nonzero(
+ (pixel_type["day"].values == 1)
+ & (pixel_type["water"].values == 1)
+ & (pixel_type["uniform"].values == 1)
+ & (pixel_type["sunglint"].values == 1)
+ & (cmin <= 0.95)
+ & (sunglint_bits == 1)
+ )
+ cmin[idx], bits["r02"]["test"][idx] = restore.sunglint(cmin, idx, bits["r02"])
+
+ # logger.debug(f"Memory: {tracemalloc.get_traced_memory()}")
+ # Shallow Water
+ # ...
+ sh_water_bits = {
+ "ir": restoral_bits["01"] * restoral_bits["05"],
+ "nir_1": restoral_bits["10"],
+ "nir_2": restoral_bits["15"],
+ "sst": restoral_bits["03"],
+ }
+ idx = np.nonzero(
+ (pixel_type["day"].values == 1)
+ & (pixel_type["water"].values == 1)
+ & (pixel_type["ice"].values == 0)
+ )
+ cmin[idx], bits["r03"]["test"][idx] = restore.shallow_water(
+ sh_water_bits, cmin, idx, bits["r03"]
+ )
+
+ desert_bits = restoral_bits["15"] * restoral_bits["05"]
+ # make sure that the land_bits are calculated properly
+ idx_desert = np.nonzero(pixel_type["desert"].values == 1)
+ # land_bits = np.zeros(desert_bits.shape)
+ temp_bit = restoral_bits["10"] + restoral_bits["13"]
+ temp_bit = np.clip(temp_bit, 0, 1)
+ land_bits = np.zeros(restoral_bits["15"].shape)
+ idx = np.nonzero(
+ (pixel_type["day"].values == 1)
+ & (pixel_type["land"].values == 1)
+ & (pixel_type["desert"].values == 0)
+ & (pixel_type["snow"].values == 0)
+ & (pixel_type["ice"].values == 0)
+ & (cmin <= 0.95)
+ )
+ idx_desert = np.nonzero(
+ (pixel_type["day"].values == 1)
+ & (pixel_type["land"].values == 1)
+ & (pixel_type["desert"].values == 1)
+ & (pixel_type["snow"].values == 0)
+ & (pixel_type["ice"].values == 0)
+ & (cmin <= 0.95)
+ )
+ tmp_bits = (
+ np.clip(restoral_bits["09"] + restoral_bits["13"], 0, 1)
+ * restoral_bits["15"]
+ * restoral_bits["05"]
+ * restoral_bits["11"]
+ )
+ land_bits[idx] = tmp_bits[idx]
+ land_bits[idx_desert] = desert_bits[idx_desert]
+ # make sure that the land_bits are calculated properly
+ idx = np.nonzero(
+ (pixel_type["day"].values == 1)
+ & (pixel_type["land"].values == 1)
+ # & (pixel_type["desert"].values == 0)
+ & (pixel_type["snow"].values == 0)
+ & (pixel_type["ice"].values == 0)
+ & (cmin <= 0.95)
+ & (land_bits == 1)
+ )
+ cmin[idx], bits["r04"]["test"][idx] = restore.land(land_bits, cmin, idx, bits["r04"])
+
+ # logger.debug(f"Memory: {tracemalloc.get_traced_memory()}")
+ # idx = np.nonzero(
+ # (pixel_type["day"] == 1)
+ # & (pixel_type["land"] == 1)
+ # & (pixel_type["desert"] == 1)
+ # & (pixel_type["snow"] == 0)
+ # & (pixel_type["ice"] == 0)
+ # & (cmin <= 0.95)
+ # & (desert_bits == 1)
+ # )
+ # 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"]
+
+ # logger.debug(f"Memory: {tracemalloc.get_traced_memory()}")
+ idx = np.nonzero(
+ (pixel_type["day"].values == 1)
+ & (pixel_type["land"].values == 1)
+ & (pixel_type["snow"].values == 0)
+ & (pixel_type["ice"].values == 0)
+ & (pixel_type["coast"].values == 1)
+ )
+ cmin[idx], bits["r05"]["test"][idx] = restore.coast(coast_bits, cmin, idx, bits["r05"])
+
+ # idx = np.nonzero((pixel_type['night'] == 1) & (pixel_type['land'] == 1) &
+ # (pixel_type['snow'] == 0) & (pixel_type['ice'] == 0) &
+ # (cmin <= 0.95))
+ # cmin[idx], bits['r06']['test'][idx] = restore.land_night(land_night_bits,
+ # cmin,
+ # bits['r06'])[idx]
+
+ cmin_final = cmin
+
+ # bits translation MVCM-python -> MVCM-C
+ # 01 13 test_11um
+ # 02 27 surface_temperature_test
+ # 03 27 sst_test
+ # 04 24 bt_diff_86_11um
+ # 05 18 test_11_12um_diff
+ # 07 19 oceanic_stratus_11_4um_test
+ # 10 20 nir_reflectance_test
+ # 11 21 vis_nir_ratio
+ # 12 23 test_16_21um_reflectance
+ # 15 16 test_1_38um_high_clouds
+ # 16 17 thin_cirrus_4_12um_BTD_test
+
+ csc = cmin_final
+ integer_cloud_mask = np.zeros(csc.shape, dtype=np.byte)
+ integer_cloud_mask[csc > 0.99] = 3
+ integer_cloud_mask[(csc <= 0.99) & (csc > 0.95)] = 2
+ integer_cloud_mask[(csc <= 0.95) & (csc > 0.66)] = 1
+ integer_cloud_mask[csc <= 0.66] = 0
+
+ # logger.debug(f"Memory: {tracemalloc.get_traced_memory()}")
+ cloud_mask = np.zeros((6, csc.shape[0], csc.shape[1]), dtype=np.byte)
+ quality_assurance = np.zeros((csc.shape[0], csc.shape[1], 10), dtype=np.ubyte)
+ # scn_id = scn.scene_id(pixel_type)
+
+ output_dict = {
+ "latitude": {"dims": ("x", "y"), "data": viirs_data.latitude.values},
+ "longitude": {"dims": ("x", "y"), "data": viirs_data.longitude.values},
+ "sensor_zenith": {"dims": ("x", "y"), "data": viirs_data.sensor_zenith.values},
+ "sensor_azimuth": {"dims": ("x", "y"), "data": viirs_data.sensor_azimuth.values},
+ "solar_zenith": {"dims": ("x", "y"), "data": viirs_data.solar_zenith.values},
+ "solar_azimuth": {"dims": ("x", "y"), "data": viirs_data.solar_azimuth.values},
+ "confidence": {"dims": ("x", "y"), "data": cmin_final},
+ "cloud_mask": {"dims": ("bits", "x", "y"), "data": cloud_mask},
+ "quality_assurance": {"dims": ("x", "y", "qa"), "data": quality_assurance},
+ "integer_cloud_mask": {"dims": ("x", "y"), "data": integer_cloud_mask},
+ }
+
+ old_debug_data = {
+ "ndvi": {"dims": ("x", "y"), "data": viirs_data.ndvi.values},
+ "sunglint_bits": {"dims": ("x", "y"), "data": sunglint_bits},
+ "eco": {"dims": ("x", "y"), "data": viirs_data.ecosystem.values},
+ "bit1": {"dims": ("x", "y"), "data": bit["01"]},
+ "bit2": {"dims": ("x", "y"), "data": bit["02"]},
+ "bit3": {"dims": ("x", "y"), "data": bit["03"]},
+ "bit4": {"dims": ("x", "y"), "data": bit["04"]},
+ "bit5": {"dims": ("x", "y"), "data": bit["05"]},
+ "bit6": {"dims": ("x", "y"), "data": bit["06"]},
+ "bit7": {"dims": ("x", "y"), "data": bit["07"]},
+ "bit8": {"dims": ("x", "y"), "data": bit["08"]},
+ "bit9": {"dims": ("x", "y"), "data": bit["09"]},
+ "bit10": {"dims": ("x", "y"), "data": bit["10"]},
+ "bit11": {"dims": ("x", "y"), "data": bit["11"]},
+ "bit12": {"dims": ("x", "y"), "data": bit["12"]},
+ "bit13": {"dims": ("x", "y"), "data": bit["13"]},
+ "bit14": {"dims": ("x", "y"), "data": bit["14"]},
+ "bit15": {"dims": ("x", "y"), "data": bit["15"]},
+ "bit16": {"dims": ("x", "y"), "data": bit["16"]},
+ "qabit1": {"dims": ("x", "y"), "data": qabit["01"]},
+ "qabit2": {"dims": ("x", "y"), "data": qabit["02"]},
+ "qabit3": {"dims": ("x", "y"), "data": qabit["03"]},
+ "qabit4": {"dims": ("x", "y"), "data": qabit["04"]},
+ "qabit5": {"dims": ("x", "y"), "data": qabit["05"]},
+ "qabit6": {"dims": ("x", "y"), "data": qabit["06"]},
+ "qabit7": {"dims": ("x", "y"), "data": qabit["07"]},
+ "qabit8": {"dims": ("x", "y"), "data": qabit["08"]},
+ "qabit9": {"dims": ("x", "y"), "data": qabit["09"]},
+ "qabit10": {"dims": ("x", "y"), "data": qabit["10"]},
+ "qabit11": {"dims": ("x", "y"), "data": qabit["11"]},
+ "qabit12": {"dims": ("x", "y"), "data": qabit["12"]},
+ "qabit13": {"dims": ("x", "y"), "data": qabit["13"]},
+ "qabit14": {"dims": ("x", "y"), "data": qabit["14"]},
+ "qabit15": {"dims": ("x", "y"), "data": qabit["15"]},
+ "qabit16": {"dims": ("x", "y"), "data": qabit["16"]},
+ "bit_r1": {"dims": ("x", "y"), "data": bit["r01"]},
+ "bit_r2": {"dims": ("x", "y"), "data": bit["r02"]},
+ "bit_r3": {"dims": ("x", "y"), "data": bit["r03"]},
+ "bit_r4": {"dims": ("x", "y"), "data": bit["r04"]},
+ "bit_r5": {"dims": ("x", "y"), "data": bit["r05"]},
+ "bit_r6": {"dims": ("x", "y"), "data": bit["r06"]},
+ # "M02": {"dims": ("x", "y"), "data": viirs_data.M02.values},
+ # "M03": {"dims": ("x", "y"), "data": viirs_data.M03.values},
+ # "M04": {"dims": ("x", "y"), "data": viirs_data.M04.values},
+ "M05": {"dims": ("x", "y"), "data": viirs_data.M05.values},
+ "M15": {"dims": ("x", "y"), "data": viirs_data.M15.values},
+ # "M15-M12": {"dims": ("x", "y"), "data": viirs_data["M15-M12"].values},
+ # "M15-M16": {"dims": ("x", "y"), "data": viirs_data["M15-M16"].values},
+ "Ocean_Day": {"dims": ("x", "y"), "data": viirs_data.Ocean_Day.values},
+ "Land_Day": {"dims": ("x", "y"), "data": viirs_data.Land_Day.values},
+ "Land_Night": {"dims": ("x", "y"), "data": viirs_data.Land_Night.values},
+ "Land_Day_Desert": {"dims": ("x", "y"), "data": viirs_data.Land_Day_Desert.values},
+ "Land_Day_Desert_Coast": {
+ "dims": ("x", "y"),
+ "data": viirs_data.Land_Day_Desert_Coast.values,
+ },
+ "Day_Snow": {"dims": ("x", "y"), "data": viirs_data.Day_Snow.values},
+ "Polar_Day_Ocean": {"dims": ("x", "y"), "data": viirs_data.Polar_Day_Ocean.values},
+ "Polar_Day_Land": {"dims": ("x", "y"), "data": viirs_data.Polar_Day_Land.values},
+ "Polar_Day_Desert": {"dims": ("x", "y"), "data": viirs_data.Polar_Day_Desert.values},
+ "Polar_Day_Desert_Coast": {
+ "dims": ("x", "y"),
+ "data": viirs_data.Polar_Day_Desert_Coast.values,
+ },
+ "elevation": {"dims": ("x", "y"), "data": viirs_data.height.values},
+ # "scene_type": {"dims": ("x", "y"), "data": scene_types},
+ }
+
+ # logger.debug(f"Memory: {tracemalloc.get_traced_memory()}")
+ if debug:
+ # debug_output.update(old_debug_data)
+ sf = {k: {"dims": ("x", "y"), "data": pixel_type[k]} for k in pixel_type.keys()}
+ old_debug_data = old_debug_data | sf
+
+ # debug_output = xr.Dataset.from_dict(old_debug_data)
+ # debug_output.to_netcdf("/ships19/hercules/pveglio/debug_data.nc", mode="w")
+
+ all_bits = np.zeros((bit["01"].shape[0], bit["01"].shape[1], 22), dtype=np.ubyte)
+ for i in range(len(_bitlist)):
+ all_bits[:, :, i] = bit[_bitlist[i]]
+
+ viirs_m_bands = SensorConstants.VIIRS_VIS_BANDS + SensorConstants.VIIRS_IR_BANDS
+
+ all_bands = np.zeros((viirs_data.M01.shape[0], viirs_data.M01.shape[1], 16), dtype=np.ubyte)
+ for j in range(len(viirs_m_bands)):
+ all_bands[:, :, j] = viirs_data[viirs_m_bands[j]].values
+
+ debug_data = {
+ "bits": {"dims": ("x", "y", "testbits"), "data": all_bits},
+ "bands": {"dims": ("x", "y", "nbands"), "data": all_bands},
+ "ecosystem": {"dims": ("x", "y"), "data": viirs_data.ecosystem.values},
+ "ndvi": {"dims": ("x", "y"), "data": viirs_data.ndvi.values},
+ # "sst": {"dims": ("x", "y"), "data": viirs_data.sst.values},
+ "tpw": {"dims": ("x", "y"), "data": viirs_data.geos_tpw.values},
+ "snow_fraction": {"dims": ("x", "y"), "data": viirs_data.geos_snow_fraction.values},
+ "ice_fraction": {"dims": ("x", "y"), "data": viirs_data.geos_ice_fraction.values},
+ "ocean_fraction": {"dims": ("x", "y"), "data": viirs_data.geos_ocean_fraction.values},
+ "land_ice_fraction": {"dims": ("x", "y"), "data": viirs_data.geos_land_ice_fraction.values},
+ "surface_temperature": {
+ "dims": ("x", "y"),
+ "data": viirs_data.geos_surface_temperature.values,
+ },
+ }
+
+ # logger.debug(f"Memory: {tracemalloc.get_traced_memory()}")
+ if debug:
+ output_dict.update(debug_data)
+
+ output_data = xr.Dataset.from_dict(output_dict)
+ attrs = {"version": _VERSION, "satellite": "VIIRS", "sensor": "VIIRS"}
+ save_output(output_data, attrs, out_file, compression=5, debug=debug)
+
+ # comp = {"zlib": True, "complevel": 5}
+ # if type(ds_out) is xr.Dataset:
+ # for var in ds_out:
+ # ds_out[var].encoding.update(comp)
+ # if type(ds_out) is xr.DataArray:
+ # ds_out.encoding.update(comp)
+ #
+ # for attr in attrs:
+ # ds_out.attrs[attr] = attrs[attr]
+ #
+ # ds_out.to_netcdf(out_file)
+
+
+if __name__ == "__main__":
+ _VERSION = get_distribution("mvcm").version
+ parser = argparse.ArgumentParser(prog="MVCM", description="")
+
+ parser.add_argument(
+ "--data",
+ help="file name and path containing satellite data and ancillary data interpolated to sensor resolution",
+ )
+ parser.add_argument(
+ "--pixel-type", help="file name and path containing all info on pixel and scene type"
+ )
+ parser.add_argument("-t", "--threshold", help="thresholds file name")
+ parser.add_argument(
+ "-v", "--verbose", action="count", default=1, help="print verbose information"
+ )
+ parser.add_argument("-o", "--out", help="output file name")
+ parser.add_argument("-d", "--debug", action="store_true", help="activate debug mode")
+ parser.add_argument("-H", "--hires", action="store_true", help="run using hires data")
+
+ args = parser.parse_args()
+
+ threshold_file = args.threshold
+ out_file = args.out
+ verbose = args.verbose or False
+ debug = args.debug or False
+ use_hires = args.hires or False
+
+ main(
+ threshold_file=threshold_file,
+ input_data=args.data,
+ pixel_data=args.pixel_type,
+ out_file=out_file,
+ verbose=verbose,
+ debug=debug,
+ use_hires=use_hires,
+ )
diff --git a/pyproject.toml b/pyproject.toml
index 0b148dbde895e3ce918789df6178c73ed59b8a59..519699ac8e535e2516697f2aabf62b5476e2b517 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -23,13 +23,14 @@ dependencies = [
'wheel',
'setuptools',
'numpy==1.26',
- 'xarray',
+ 'xarray==2023.6.0',
'attrs',
'netCDF4',
'rich',
'ruamel.yaml',
'pre-commit',
'dask',
+ 'psutil',
]
dynamic = ['version']