From 0954b89153feb67832ca484dc2220e7152b05fe0 Mon Sep 17 00:00:00 2001
From: Paolo Veglio <paolo.veglio@ssec.wisc.edu>
Date: Fri, 16 Dec 2022 22:21:31 +0000
Subject: [PATCH] fixed 11-12um BT test. Forgot to flip the res flag in the
function
---
hires_module.py | 56 +++++++++++++++++++++++++++------------
read_data.py | 8 ++++--
restoral.py | 2 +-
scene.py | 12 ++++-----
sort_inputs.py | 67 +++++++++++++++++++++++++++++++++++++----------
spectral_tests.py | 14 +++++-----
6 files changed, 113 insertions(+), 46 deletions(-)
diff --git a/hires_module.py b/hires_module.py
index 2f2193b..4b6f815 100644
--- a/hires_module.py
+++ b/hires_module.py
@@ -14,13 +14,17 @@ import restoral
import read_data as rd
import scene as scn
+import importlib
+
+importlib.reload(tst)
+
# #################################################################### #
# TEST CASE
# data:
_datapath = '/ships19/hercules/pveglio/mvcm_viirs_hires'
-_fname_mod02 = glob(f'{_datapath}/VNP02MOD.A2022173.1454.001.*.uwssec_bowtie_restored.nc')[0]
+_fname_mod02 = glob(f'{_datapath}/VNP02MOD.A2022173.1454.001.*.uwssec_bowtie_corrected.nc')[0]
_fname_mod03 = glob(f'{_datapath}/VNP03MOD.A2022173.1454.001.*.uwssec.nc')[0]
-_fname_img02 = glob(f'{_datapath}/VNP02IMG.A2022173.1454.001.*.uwssec_bowtie_restored.nc')[0]
+_fname_img02 = glob(f'{_datapath}/VNP02IMG.A2022173.1454.001.*.uwssec_bowtie_corrected.nc')[0]
_fname_img03 = glob(f'{_datapath}/VNP03IMG.A2022173.1454.001.*.uwssec.nc')[0]
# MVCM output
@@ -63,15 +67,22 @@ def main(*,
cm = f['geophysical_data/Cloud_Mask'][:]
f.close()
+ check = nc.Dataset(mod02)
+ if 'M01' not in list(check['observation_data'].variables):
+ check.close()
+ return
+
+ if check.isopen():
+ check.close()
with open(threshold_file) as f:
text = f.read()
thresholds = yml.safe_load(text)
- file_names = {'MOD02': f'{mod02}',
- 'MOD03': f'{mod03}',
- 'IMG02': f'{img02}',
- 'IMG03': f'{img03}',
+ file_names = {'MOD02': f'{data_path}/{mod02}',
+ 'MOD03': f'{data_path}/{mod03}',
+ 'IMG02': f'{data_path}/{img02}',
+ 'IMG03': f'{data_path}/{img03}',
'GEOS_atm_1': f'{geos_atm_1}',
'GEOS_atm_2': f'{geos_atm_2}',
'GEOS_land': f'{geos_land}',
@@ -85,6 +96,8 @@ def main(*,
sunglint_angle = thresholds['Sun_Glint']['bounds'][3]
viirs_data = rd.get_data(file_names, sunglint_angle, hires=True)
+ if 'I01' not in viirs_data:
+ return
cm_flag = np.array((np.array(cm[0, :, :], np.byte) &
(pow(2, 1) + pow(2, 2)))/pow(2, 1))
@@ -95,7 +108,7 @@ def main(*,
mvcm_conf[cm_flag == 1] = 0.66
mvcm_conf[cm_flag == 0] = 0
- viirs_lores = rd.read_data('viirs', mod02, mod03)
+ viirs_lores = rd.read_data('viirs', mod02, data_path + '/' + mod03)
viirs_data['M02'] = (('number_of_lines_2', 'number_of_pixels_2'),
viirs_lores.M02.values.repeat(2, 0).repeat(2, 1))
@@ -104,7 +117,10 @@ def main(*,
viirs_data['M01'] = (('number_of_lines_2', 'number_of_pixels_2'), np.zeros(viirs_data.I01.shape))
viirs_data['M10'] = viirs_data.I03
viirs_data['M12'] = viirs_data.I04
- viirs_data['M15'] = viirs_data.I05
+ viirs_data['M15'] = (('number_of_lines_2', 'number_of_pixels_2'),
+ viirs_lores.M15.values.repeat(2, 0).repeat(2, 1))
+ viirs_data['M16'] = (('number_of_lines_2', 'number_of_pixels_2'),
+ viirs_lores.M16.values.repeat(2, 0).repeat(2, 1))
cmin_G1 = np.ones(viirs_data.M15.shape)
cmin_G2 = np.ones(viirs_data.M15.shape)
@@ -112,18 +128,21 @@ def main(*,
MyScene = tst.CloudTests(viirs_data, 'Ocean_Day', thresholds)
- cmin_G1, bit1 = MyScene.test_11um('M15', cmin_G1, test_name='11um_Test')
+ cmin_G1, bit1 = MyScene.test_11um('I05', cmin_G1, test_name='11um_Test')
+ cmin_G1, bit2 = MyScene.sst_test('M15', 'M16', cmin_G1, test_name='SST_Test')
- cmin_G2, bit2 = MyScene.oceanic_stratus_11_4um_test('I05-I04', cmin_G2,
+ cmin_G2, bit3 = MyScene.oceanic_stratus_11_4um_test('I05-I04', cmin_G2,
test_name='11-4um_Oceanic_Stratus_Test')
- cmin_G3, bit3 = MyScene.nir_reflectance_test('M07', cmin_G3, test_name='NIR_Reflectance_Test')
- cmin_G3, bit4 = MyScene.vis_nir_ratio_test('I02-I01ratio', cmin_G3, test_name='Vis/NIR_Ratio_Test')
- cmin_G3, bit5 = MyScene.nir_reflectance_test('M10', cmin_G3, test_name='1.6_2.1um_NIR_Reflectance_Test')
+ cmin_G3, bit4 = MyScene.nir_reflectance_test('M07', cmin_G3, test_name='NIR_Reflectance_Test')
+ cmin_G3, bit5 = MyScene.vis_nir_ratio_test('I02-I01ratio', cmin_G3, test_name='Vis/NIR_Ratio_Test')
+ visnir, bit5 = MyScene.vis_nir_ratio_test('I02-I01ratio', np.ones(viirs_data.M15.shape), test_name='Vis/NIR_Ratio_Test')
+ cmin_G3, bit6 = MyScene.nir_reflectance_test('M10', cmin_G3, test_name='1.6_2.1um_NIR_Reflectance_Test')
confidence = np.power(cmin_G1 * cmin_G2 * cmin_G3, 1/3)
temp_confidence = np.power(cmin_G1 * cmin_G2 * cmin_G3, 1/3)
- total_bit = np.ones(confidence.shape)
+ total_bit = np.zeros(confidence.shape)
+ total_bit = bit1 + bit2
sunglint_angle = thresholds['Sun_Glint']['bounds'][3]
scene_flags = scn.find_scene(viirs_data, sunglint_angle)
@@ -143,9 +162,12 @@ def main(*,
date = mod02.split('.')[1]
time = mod02.split('.')[2]
- np.savez(f'{data_path}/output/pyMVCM_{date}.{time}',
- confidence=total_confidence, lat=viirs_data.latitude.values, lon=viirs_data.longitude.values,
- ocean_day_mask=viirs_data.Ocean_Day.values)
+ # np.savez(f'{data_path}/outputs/pyMVCM_{date}.{time}',
+ np.savez(f'{data_path}/pyMVCM_{date}.{time}',
+ confidence=total_confidence, py_conf=confidence, mvcm_conf=mvcm_conf,
+ lat=viirs_data.latitude.values, lon=viirs_data.longitude.values,
+ ocean_day_mask=viirs_data.Ocean_Day.values, visnir=visnir,
+ c1=cmin_G1, c2=cmin_G2, c3=cmin_G3)
# return confidence # , conf, total_bit
diff --git a/read_data.py b/read_data.py
index 1d89d4c..70a0dff 100644
--- a/read_data.py
+++ b/read_data.py
@@ -4,6 +4,8 @@ import numpy as np
import ancillary_data as anc
import scene as scn
+import os
+from datetime import datetime as dt
from typing import Dict
_DTR = np.pi/180.
@@ -80,7 +82,6 @@ def read_ancillary_data(file_names: str,
resolution: int = 1) -> xr.Dataset:
# Ancillary files temporarily defined here. Eventually we will find a better way to pass these
- start_time = '2022-06-22 14:54:00.000'
anc_dir = file_names['ANC_DIR']
sst_file = file_names['SST']
ndvi_file = file_names['NDVI']
@@ -89,6 +90,9 @@ def read_ancillary_data(file_names: str,
geos_ocn = file_names['GEOS_ocean']
geos1 = file_names['GEOS_atm_1']
geos2 = file_names['GEOS_atm_2']
+ vnptime = '.'.join(os.path.basename(file_names['MOD02']).split('.')[1:3])
+ start_time = dt.strftime(dt.strptime(vnptime, 'A%Y%j.%H%M'), '%Y-%m-%d %H:%M:%S.000')
+ # start_time = '2022-06-22 14:54:00.000'
dim1 = viirs_data.latitude.shape[0]
dim2 = viirs_data.latitude.shape[1]
@@ -179,7 +183,7 @@ def get_data(file_names: Dict[str, str],
viirs_data = read_data('viirs', f'{img02}', f'{img03}')
viirs_data['M05'] = viirs_data.I01
viirs_data['M07'] = viirs_data.I02
- viirs_data = read_ancillary_data(file_names, viirs_data)
+ viirs_data = read_ancillary_data(file_names, viirs_data, resolution=2)
viirs_data['I05-I04'] = (('number_of_lines_2', 'number_of_pixels_2'),
viirs_data.I05.values - viirs_data.I04.values)
diff --git a/restoral.py b/restoral.py
index f01d7e8..df3f404 100644
--- a/restoral.py
+++ b/restoral.py
@@ -24,7 +24,7 @@ def sunglint(viirs_data, threshold, bit, conf):
m02 = viirs_data.M07.values
irclr = np.ones(viirs_data.M02.shape)
- # irclr[bit < 1] = 0
+ irclr[bit < 2] = 0
var = spatial_var(m31, 0.4)
reg_var_mean = sliding_window_view(np.pad(m02, [1, 1], mode='constant'),
diff --git a/scene.py b/scene.py
index 4c78603..0e37c57 100644
--- a/scene.py
+++ b/scene.py
@@ -39,19 +39,19 @@ night = 200
# TEST CASE
# data:
datapath = '/ships19/hercules/pveglio/mvcm_viirs_hires'
-fname_mod02 = glob(f'{datapath}/VNP02MOD.A2022173.1448.001.*.uwssec.nc')[0]
+fname_mod02 = glob(f'{datapath}/VNP02MOD.A2022173.1448.001.*.uwssec*.nc')[0]
fname_mod03 = glob(f'{datapath}/VNP03MOD.A2022173.1448.001.*.uwssec.nc')[0]
-fname_img02 = glob(f'{datapath}/VNP02IMG.A2022173.1448.001.*.uwssec.nc')[0]
+fname_img02 = glob(f'{datapath}/VNP02IMG.A2022173.1448.001.*.uwssec*.nc')[0]
fname_img03 = glob(f'{datapath}/VNP03IMG.A2022173.1448.001.*.uwssec.nc')[0]
# thresholds:
threshold_file = '/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.yaml'
# ancillary files:
-geos_atm_1 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1500.V01.nc4'
-geos_atm_2 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1800.V01.nc4'
-geos_land = 'GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20220622_1630.V01.nc4'
-geos_ocean = 'GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_1630.V01.nc4'
+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_1430.V01.nc4'
+geos_ocean = 'GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_1430.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'
diff --git a/sort_inputs.py b/sort_inputs.py
index 5f547d2..282fb91 100644
--- a/sort_inputs.py
+++ b/sort_inputs.py
@@ -5,8 +5,10 @@ import numpy as np
from datetime import datetime
from glob import glob
+import hires_module as hm
-def main(vnp02mod_fname):
+
+def create_filenames(vnp02mod_fname):
datapath = os.path.dirname(vnp02mod_fname)
filename = os.path.basename(vnp02mod_fname)
@@ -17,7 +19,7 @@ def main(vnp02mod_fname):
geos_times = ['0000', '0300', '0600', '0900', '1200', '1500', '1800', '2100']
geos_date_from_viirs = datetime.strftime(datetime.strptime(vnpdate, 'A%Y%j'), '%Y%m%d')
geos_flist = glob(f'{datapath}/ancillary/GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.{geos_date_from_viirs}*.nc4')
-
+ geos_flist = np.sort(geos_flist)
fmt = '%H%M'
diff_times = [(datetime.strptime(gt, fmt) - datetime.strptime(vnptime, fmt)).total_seconds()
for gt in geos_times]
@@ -37,21 +39,58 @@ def main(vnp02mod_fname):
for gt in land_ocean_timelist]
land_ocean_index = np.argmin(np.abs(diff_times))
-
- geos_land_file = glob(f'{datapath}/ancillary/GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.{geos_date_from_viirs}*.nc4')[land_ocean_index]
- geos_ocean_file = glob(f'{datapath}/ancillary/GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.{geos_date_from_viirs}*.nc4')[land_ocean_index]
+ geos_land_flist = glob(f'{datapath}/ancillary/GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.{geos_date_from_viirs}*.nc4')
+ geos_land_file = np.sort(geos_land_flist)[land_ocean_index]
+ geos_ocean_flist = glob(f'{datapath}/ancillary/GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.{geos_date_from_viirs}*.nc4')
+ geos_ocean_file = np.sort(geos_ocean_flist)[land_ocean_index]
vnp03mod = glob(f'{datapath}/VNP03MOD.{vnpdate}.{vnptime}*.nc')[0]
- vnp02img = glob(f'{datapath}/VNP02IMG.{vnpdate}.{vnptime}*_bowtie_restored.nc')[0]
+ vnp02img = glob(f'{datapath}/VNP02IMG.{vnpdate}.{vnptime}*_bowtie_corrected.nc')[0]
vnp03img = glob(f'{datapath}/VNP03IMG.{vnpdate}.{vnptime}*.nc')[0]
-
- out_fnames = {'VNP03MOD': vnp03mod,
- 'VNP02IMG': vnp02img,
- 'VNP03IMG': vnp03img,
- 'GEOS_atm_1': geos_file1,
- 'GEOS_atm_2': geos_file2,
- 'GEOS_land': geos_land_file,
- 'GEOS_ocean': geos_ocean_file,
+ cld_msk = glob(f'{datapath}/CLDMSK_L2_VIIRS_SNPP.{vnpdate}.{vnptime}*.nc')[0]
+
+ out_fnames = {'VNP03MOD': os.path.basename(vnp03mod),
+ 'VNP02IMG': os.path.basename(vnp02img),
+ 'VNP03IMG': os.path.basename(vnp03img),
+ 'GEOS_atm_1': os.path.basename(geos_file1),
+ 'GEOS_atm_2': os.path.basename(geos_file2),
+ 'GEOS_land': os.path.basename(geos_land_file),
+ 'GEOS_ocean': os.path.basename(geos_ocean_file),
+ 'GEOS_constants': 'GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4',
+ 'SST_file': 'oisst.20220622',
+ 'ECO_file': 'goge1_2_img.v1',
+ 'NDVI_file': 'NDVI.FM.c004.v2.0.WS.00-04.177.hdf',
+ 'data_path': datapath,
+ 'cloud_mask_file': cld_msk
}
return out_fnames
+
+
+def call_mvcm():
+
+ # flist = np.sort(glob('/ships19/hercules/pveglio/mvcm_viirs_hires/VNP02MOD*.nc'))
+ flist = np.sort(glob('/ships19/hercules/pveglio/mvcm_viirs_hires/VNP02MOD.A2022173.1454*.nc'))
+
+ for fname in flist:
+ print(f'Processing {fname}...')
+ fnames_for_mvcm = create_filenames(fname)
+
+ hm.main(data_path=fnames_for_mvcm['data_path'],
+ mod02=fname,
+ mod03=fnames_for_mvcm['VNP03MOD'],
+ img02=fnames_for_mvcm['VNP02IMG'],
+ img03=fnames_for_mvcm['VNP03IMG'],
+ threshold_file=hm._threshold_file,
+ geos_atm_1=fnames_for_mvcm['GEOS_atm_1'],
+ geos_atm_2=fnames_for_mvcm['GEOS_atm_2'],
+ geos_land=fnames_for_mvcm['GEOS_land'],
+ geos_ocean=fnames_for_mvcm['GEOS_ocean'],
+ geos_constants=fnames_for_mvcm['GEOS_constants'],
+ ndvi_file=fnames_for_mvcm['NDVI_file'],
+ sst_file=fnames_for_mvcm['SST_file'],
+ cloud_mask_file=fnames_for_mvcm['cloud_mask_file'])
+
+
+if __name__ == "__main__":
+ call_mvcm()
diff --git a/spectral_tests.py b/spectral_tests.py
index 4c2962b..b09130a 100644
--- a/spectral_tests.py
+++ b/spectral_tests.py
@@ -89,7 +89,7 @@ class CloudTests(object):
cmin: np.ndarray,
test_name: str = 'SST_Test') -> np.ndarray:
- confidence = np.ones(self.data.M01.shape)
+ confidence = np.ones(self.data[band31].shape)
qa_bit = np.zeros(self.data[band31].shape)
test_bit = np.zeros(self.data[band31].shape)
threshold = self.thresholds[self.scene_name][test_name]
@@ -98,13 +98,13 @@ class CloudTests(object):
qa_bit[self.scene_idx] = 1
m31 = self.data[band31].values - 273.16
bt_diff = self.data[band31].values - self.data[band32].values
- sst = self.data.geos_sfct.values - 273.16
+ sst = self.data.sst.values - 273.16
cosvza = np.cos(self.data.sensor_zenith.values*_DTR)
c = threshold['coeffs']
modsst = 273.16 + c[0] + c[1]*m31 + c[2]*bt_diff*sst + c[3]*bt_diff*((1.0/cosvza) - 1.0)
- sfcdif = self.data.geos_sfct.values - modsst
+ sfcdif = self.data.sst.values - modsst
idx = np.nonzero((sfcdif < threshold['thr'][1]) &
(self.data[self.scene_name] == 1))
@@ -122,7 +122,7 @@ class CloudTests(object):
cmin: np.ndarray,
test_name: str = '8.6-11um_Test') -> np.ndarray:
- confidence = np.ones(self.data.M01.shape)
+ confidence = np.ones(self.data[band].shape)
qa_bit = np.zeros(self.data[band].shape)
test_bit = np.zeros(self.data[band].shape)
threshold = self.thresholds[self.scene_name][test_name]
@@ -566,13 +566,15 @@ class ComputeTests(CloudTests):
sunglint_angle = self.thresholds['Sun_Glint']['bounds'][3]
scene_flags = scn.find_scene(self.data, sunglint_angle)
+ cmin_tmp = cmin + 0
idx = np.nonzero((scene_flags['water'] == 1) & (scene_flags['ice'] == 0) &
(scene_flags['uniform'] == 1) & (cmin <= 0.99) & (cmin >= 0.05))
- cmin[idx] = restoral.spatial(self.data, self.thresholds['Sun_Glint'], scene_flags, cmin)[idx]
+ cmin[idx] = restoral.spatial(self.data, self.thresholds['Sun_Glint'], scene_flags, cmin_tmp)[idx]
+ cmin_tmp = cmin + 0
idx = np.nonzero((scene_flags['water'] == 1) & (scene_flags['sunglint'] == 1) &
(scene_flags['uniform'] == 1) & (cmin <= 0.95))
- cmin[idx] = restoral.sunglint(self.data, self.thresholds['Sun_Glint'], total_bit, cmin)[idx]
+ cmin[idx] = restoral.sunglint(self.data, self.thresholds['Sun_Glint'], total_bit, cmin_tmp)[idx]
"""
idx = np.nonzero((scene_flags['day'] == 1) & (scene_flags['land'] == 1) &
--
GitLab