diff --git a/mvcm/main.py b/mvcm/main.py
index ce38334d63ad73cea97209ae6069d0f224209730..db019316868e2286de58e949ffe81abecd558256 100644
--- a/mvcm/main.py
+++ b/mvcm/main.py
@@ -3,6 +3,9 @@ try:
except ImportError:
import ruamel_yaml as yml
+import functools
+import time
+
import numpy as np
import xarray as xr
@@ -41,6 +44,19 @@ _eco_file = 'goge1_2_img.v1'
# #################################################################### #
+def timer(func):
+ @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
+ print(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',
data_path: str = _datapath,
@@ -79,6 +95,7 @@ def main(satellite: str = 'snpp',
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
@@ -135,6 +152,8 @@ def main(satellite: str = 'snpp',
# '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']:
+ temp_bit = np.zeros((18, viirs_data.M11.shape[0], viirs_data.M11.shape[1]))
+ thr_temp = np.ones(viirs_data.M11.shape) * -12345.6
for scene_name in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert', 'Land_Day_Desert_Coast',
'Ocean_Day']:
print(f'Processing {scene_name}')
@@ -160,7 +179,7 @@ def main(satellite: str = 'snpp',
cmin_G1, bit['03'] = MyScene.sst_test('M15', 'M16', cmin_G1)
cmin_G2, bit['04'] = MyScene.bt_diff_86_11um('M14-M15', cmin_G2)
- cmin_G2, bit['05'] = MyScene.test_11_12um_diff('M15-M16', cmin_G2)
+ cmin_G2, bit['05'], thr_temp = MyScene.test_11_12um_diff('M15-M16', cmin_G2, thr_temp)
cmin_G2, bit['06'] = MyScene.variability_11um_test('M15', cmin_G2)
cmin_G2, bit['07'] = MyScene.bt_difference_11_4um_test_ocean('M15-M12', cmin_G2)
# # cmin_G2, bit['08'] = MyScene.bt_difference_11_4um_test_land('M15-M12', cmin_G2)
@@ -181,8 +200,10 @@ def main(satellite: str = 'snpp',
cmin_G5, bit['16'] = MyScene.thin_cirrus_4_12um_BTD_test('M13-M16', cmin_G5)
# 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_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(bit))
sunglint_bits = bit['01'] * bit['03'] * bit['05'] * bit['15']
@@ -200,37 +221,39 @@ def main(satellite: str = 'snpp',
idx = np.nonzero((scene_flags['uniform'] == 1) & (scene_flags['water'] == 1) &
(scene_flags['ice'] == 0) & (cmin >= 0.05) & (cmin <= 0.99))
- cmin[idx] = Restore.spatial_variability(cmin)[idx]
+ cmin[idx] = Restore.spatial_variability(copy(cmin))[idx]
idx = np.nonzero((scene_flags['day'] == 1) & (scene_flags['water'] == 1) &
(scene_flags['uniform'] == 1) & (scene_flags['sunglint'] == 1) &
(cmin <= 0.95))
- # cmin[idx] = Restore.sunglint(sunglint_bits, copy(cmin))[idx]
+ var_threshold = thresholds['Daytime_Ocean_Spatial_Variability']['var_11um']
+ cmin[idx] = Restore.sunglint(sunglint_bits, var_threshold, copy(cmin))[idx]
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, copy(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, copy(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, copy(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, copy(cmin))[idx]
cmin_3d[i, :, :] = cmin
# cmin_3d[i, :, :] = np.power(cmin_G1 * cmin_G2 * cmin_G3 * cmin_G4 * cmin_G5, 1/utils.group_count(bit))
i += 1
cmin_final = np.min(cmin_3d, axis=0)
+ # bit['15'] = np.max(temp_bit, axis=0)
'''
sunglint_bits = bit['01'] * bit['03'] * bit['05'] * bit['15']
sh_water_bits = {'ir': bit['01'] * bit['05'],
@@ -322,8 +345,14 @@ def main(satellite: str = 'snpp',
'data': bit['03']},
'bit5': {'dims': ('x', 'y'),
'data': bit['05']},
+ 'bit10': {'dims': ('x', 'y'),
+ 'data': bit['10']},
'bit15': {'dims': ('x', 'y'),
'data': bit['15']},
+ 'thr': {'dims': ('x', 'y'),
+ 'data': thr_temp},
+ # 'shallow_water_bits': {'dims': ('x', 'y'),
+ # 'data': sh_water_bits},
}
sf = {k: {'dims': ('x', 'y'), 'data': scene_flags[k]} for k in scene_flags}
diff --git a/mvcm/preprocess_thresholds.py b/mvcm/preprocess_thresholds.py
index d99fa92fe7a9e293be044d012ef23695aef0e86c..ed78068c5229680584624e9424f98d5599fa32ea 100644
--- a/mvcm/preprocess_thresholds.py
+++ b/mvcm/preprocess_thresholds.py
@@ -37,31 +37,33 @@ def prepare_11_12um_thresholds(thresholds: np.ndarray,
return thr_dict
-# function was called preproc
def thresholds_11_12um(data: xr.Dataset,
thresholds: np.ndarray,
scene: str,
scene_idx: np.ndarray) -> np.ndarray:
cosvza = np.cos(data.sensor_zenith.values[scene_idx].ravel() * _DTR)
schi = np.full(cosvza.shape, 99.0)
- schi[cosvza > 0] = 1/cosvza
+ schi[np.abs(cosvza) > 0] = 1.0/cosvza
schi = np.array(schi, dtype=np.float32) # this is because the C function expects a float
- m15 = data.M15.values[scene_idx].ravel()
+ m15 = np.array(data.M15.values[scene_idx].ravel(), dtype=np.float32)
latitude = data.latitude.values[scene_idx].ravel()
thr = anc.py_cithr(1, schi, m15)
thr_dict = prepare_11_12um_thresholds(thresholds, m15.shape[0])
midpt = np.full(m15.shape[0], thr)
- idx = np.nonzero((thr < 0.1) | (np.abs(schi-99) < 0.0001))
+ idx = np.nonzero((thr < 0.1) | (np.abs(schi-99.0) < 0.0001))
midpt[idx] = thr_dict['coeffs'][idx, 0]
locut = midpt + (thr_dict['cmult'] * midpt)
+
if scene in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert', 'Land_Day_Desert_Coast',
'Ocean_Day', 'Ocean_Night', 'Polar_Day_Ocean', 'Polar_Night_Ocean']:
hicut = midpt - thr_dict['adj']
+
elif scene in ['Polar_Day_Land', 'Polar_Day_Coast', 'Polar_Day_Desert',
'Polar_Day_Desert_Coast', 'Polar_Day_Snow']:
hicut = midpt - (thr_dict['adj'] * midpt)
+
elif scene in ['Land_Night', 'Polar_Night_Land', 'Polar_Night_Snow', 'Day_Snow', 'Night_Snow']:
_coeffs = {'Land_Night': 0.3, 'Polar_Night_Land': 0.3, 'Polar_Night_Snow': 0.3,
'Day_Snow': 0.0, 'Night_Snow': 0.3}
diff --git a/mvcm/utility_functions.py b/mvcm/utility_functions.py
index f5091897165fee63eec27cb2a3a8f607e9d04fe6..c53f5a89a772f3e5fe6b168b1ba89a97dfae0e46 100644
--- a/mvcm/utility_functions.py
+++ b/mvcm/utility_functions.py
@@ -25,7 +25,7 @@ def local_nxn_standard_deviation(arr: np.ndarray,
def spatial_var(rad: np.ndarray,
threshold: Dict) -> np.ndarray:
- """Compute 3x3 spatial variability and returns number of pixels within certain threshold"""
+ """Compute 3x3 spatial variability and returns fractional number of pixels within the given threshold"""
test = sliding_window_view(np.pad(rad, [1, 1], mode='constant'), (3, 3)) - np.expand_dims(rad, (2, 3))
test = np.abs(test.reshape(rad.shape[0], rad.shape[1], 9))