-
Paolo Veglio authoredPaolo Veglio authored
main.py 26.13 KiB
"""Main function for MVCM."""
import argparse
import logging
import numpy as np
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
import mvcm.scene as scn
import mvcm.spectral_tests as tst
import mvcm.utility_functions as utils
from mvcm.restoral import Restoral
from mvcm.write_output import 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"]
# #################################################################### #
# TEST CASE
# data:
_datapath = "/ships19/hercules/pveglio/mvcm_viirs_hires"
# _fname_mod02 = glob(
# f'{_datapath}/VNP02MOD.A2022173.1312.001.*.uwssec_bowtie_corrected.nc')[0]
# _fname_mod03 = glob(
# f'{_datapath}/VNP03MOD.A2022173.1312.001.*.uwssec.nc')[0]
# _fname_img02 = glob(
# f'{_datapath}/VNP02IMG.A2022173.1312.001.*.uwssec_bowtie_corrected.nc')[0]
# _fname_img03 = glob(
# f'{_datapath}/VNP03IMG.A2022173.1312.001.*.uwssec.nc')[0]
# thresholds:
_threshold_file = "/home/pveglio/mvcm/thresholds/thresholds.mvcm.snpp.v0.0.1.yaml"
# ancillary files:
_geos_atm_1 = "GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1200.V01.nc4"
_geos_atm_2 = "GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1500.V01.nc4"
_geos_land = "GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20220622_1330.V01.nc4"
_geos_ocean = "GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_1330.V01.nc4"
_geos_constants = "GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4"
_ndvi_file = "NDVI.FM.c004.v2.0.WS.00-04.177.hdf"
_sst_file = "oisst.20220622"
_eco_file = "goge1_2_img.v1"
# #################################################################### #
_scene_list = [
"Land_Day",
"Land_Day_Coast",
"Land_Day_Desert",
"Land_Day_Desert_Coast",
"Ocean_Day",
"Ocean_Night",
"Polar_Day_Ocean",
"Polar_Night_Ocean",
"Polar_Day_Land",
"Polar_Day_Coast",
"Polar_Day_Desert",
"Polar_Day_Desert_Coast",
"Polar_Day_Snow",
"Land_Night",
"Polar_Night_Land",
"Polar_Night_Snow",
"Day_Snow",
"Night_Snow",
]
_mod_bands = [
"M01",
"M02",
"M03",
"M04",
"M05",
"M06",
"M07",
"M08",
"M09",
"M10",
"M11",
"M12",
"M13",
"M14",
"M15",
"M16",
]
_img_bands = ["I01", "I02", "I03", "I04", "I05"]
def main(
satellite: str = "snpp",
sensor: str = "viirs",
data_path: str = _datapath,
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,
geos_land: str = _geos_land,
geos_ocean: str = _geos_ocean,
geos_constants: str = _geos_constants,
ndvi_file: str = _ndvi_file,
sst_file: str = _sst_file,
eco_file: str = _eco_file,
out_file: str = "temp_out.nc",
verbose: int = 0,
debug: bool = False, # noqa
) -> 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
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}/ancillary",
}
with open(threshold_file) as f:
text = f.read()
thresholds = YAML(typ="safe").load(text)
# 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 _img_bands:
try:
vnp02[b]
except KeyError:
logger.info(f"Band {b} not found in file. No output will be written.")
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 _mod_bands:
try:
vnp02[b]
except KeyError:
logger.info(f"Band {b} not found in file. No output will be written.")
return
logger.info(f"All bands found in file {file_names['MOD02']}. The code will run.")
with Dataset(file_names["MOD03"]) as f:
attrs = {attr: f.getncattr(attr) for attr in f.ncattrs()}
viirs_data = rd.get_data(satellite, sensor, file_names, thresholds, hires=use_hires)
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)
cmin_3d = np.ones((18, viirs_data.M11.shape[0], viirs_data.M11.shape[1]), dtype=np.float32)
##########################################################
_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.M11.shape, dtype=np.int8),
"qa": np.zeros(viirs_data.M11.shape, dtype=np.int8),
}
i = 0
viirs_data["locut"] = xr.DataArray(
np.zeros(viirs_data.M11.shape) * np.nan, dims=viirs_data.M11.dims
)
viirs_data["midpt"] = xr.DataArray(
np.zeros(viirs_data.M11.shape) * np.nan, dims=viirs_data.M11.dims
)
viirs_data["hicut"] = xr.DataArray(
np.zeros(viirs_data.M11.shape) * np.nan, dims=viirs_data.M11.dims
)
scene_types = np.zeros(viirs_data.M11.shape, dtype=np.int8)
for _scene_i, scene_name in enumerate(_scene_list):
scene_types[viirs_data[scene_name].values == 1] = i
logger.info(f"Processing {scene_name}")
if np.all(viirs_data[scene_name].values == 0):
logger.info("Skipping, no pixels in scene.")
continue
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, dtype=np.float32)
cmin_g2 = np.ones(viirs_data.M11.shape, dtype=np.float32)
cmin_g3 = np.ones(viirs_data.M11.shape, dtype=np.float32)
cmin_g4 = np.ones(viirs_data.M11.shape, dtype=np.float32)
cmin_g5 = np.ones(viirs_data.M11.shape, dtype=np.float32)
# cmin_temp = np.ones(viirs_data.M11.shape)
if use_hires is True:
m15_name = "M15hi"
else:
m15_name = "M15"
# 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_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"])
# Group 5
cmin_g5, bits["16"] = my_scene.thin_cirrus_4_12um_BTD_test("M13-M16", cmin_g5, bits["16"])
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_3d[i, :, :] = np.fmin(cmin_temp,
np.power(cmin_g1 * cmin_g2 * cmin_g3 * cmin_g4 * cmin_g5,
1/utils.group_count(bit)))
"""
cmin = np.power(
cmin_g1 * cmin_g2 * cmin_g3 * cmin_g4 * cmin_g5, 1 / utils.group_count(qabit)
)
cmin_3d[i, :, :] = cmin
i += 1
cmin = np.min(cmin_3d, axis=0)
# 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(
(scene_flags["uniform"] == 1)
& (scene_flags["water"] == 1)
& (scene_flags["ice"] == 0)
& (cmin >= 0.05)
& (cmin <= 0.99)
)
cmin[idx], bits["r01"]["test"][idx] = restore.spatial_variability(cmin, idx, bits["r01"])
# Sun Glint
sunglint_bits = (
restoral_bits["01"] * restoral_bits["03"] * restoral_bits["05"] * restoral_bits["15"]
)
idx = np.nonzero(
(scene_flags["day"] == 1)
& (scene_flags["water"] == 1)
& (scene_flags["uniform"] == 1)
& (scene_flags["sunglint"] == 1)
& (cmin <= 0.95)
& (sunglint_bits == 1)
)
cmin[idx], bits["r02"]["test"][idx] = restore.sunglint(cmin, idx, bits["r02"])
# 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(
(scene_flags["day"] == 1) & (scene_flags["water"] == 1) & (scene_flags["ice"] == 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(scene_flags["desert"] == 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(
(scene_flags["day"] == 1)
& (scene_flags["land"] == 1)
& (scene_flags["desert"] == 0)
& (scene_flags["snow"] == 0)
& (scene_flags["ice"] == 0)
& (cmin <= 0.95)
)
idx_desert = np.nonzero(
(scene_flags["day"] == 1)
& (scene_flags["land"] == 1)
& (scene_flags["desert"] == 1)
& (scene_flags["snow"] == 0)
& (scene_flags["ice"] == 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(
(scene_flags["day"] == 1)
& (scene_flags["land"] == 1)
# & (scene_flags["desert"] == 0)
& (scene_flags["snow"] == 0)
& (scene_flags["ice"] == 0)
& (cmin <= 0.95)
& (land_bits == 1)
)
cmin[idx], bits["r04"]["test"][idx] = restore.land(land_bits, cmin, idx, bits["r04"])
# idx = np.nonzero(
# (scene_flags["day"] == 1)
# & (scene_flags["land"] == 1)
# & (scene_flags["desert"] == 1)
# & (scene_flags["snow"] == 0)
# & (scene_flags["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"]
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], 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,
# 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)
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
cloud_mask = np.zeros((6, csc.shape[0], csc.shape[1]))
quality_assurance = np.zeros((csc.shape[0], csc.shape[1], 10))
# scn_id = scn.scene_id(scene_flags)
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},
}
if debug is True:
# debug_output.update(old_debug_data)
sf = {k: {"dims": ("x", "y"), "data": scene_flags[k]} for k in scene_flags}
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))
for i in range(len(_bitlist)):
all_bits[:, :, i] = bit[_bitlist[i]]
all_bands = np.zeros((viirs_data.M01.shape[0], viirs_data.M01.shape[1], 16))
for j in range(len(_mod_bands)):
all_bands[:, :, j] = viirs_data[_mod_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,
},
}
if debug is True:
output_dict.update(debug_data)
output_data = xr.Dataset.from_dict(output_dict)
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(
"--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")
args = parser.parse_args()
satellite = args.satellite or "snpp"
sensor = args.sensor or "viirs"
data_path = args.path or _datapath
mod02 = args.l1b
mod03 = args.geolocation
img02 = args.hires_l1b or None
img03 = args.hires_geo or None
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
geos_land = args.land or _geos_land
geos_ocean = args.ocean or _geos_ocean
constants = args.constants or _geos_constants
ndvi_file = args.ndvi or _ndvi_file
sst_file = args.sst or _sst_file
eco_file = args.eco or _eco_file
out_file = args.out or "test_out.nc"
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,
)