Skip to content
Snippets Groups Projects
Commit 0954b891 authored by Paolo Veglio's avatar Paolo Veglio
Browse files

fixed 11-12um BT test. Forgot to flip the res flag in the function

parent ea73ea41
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
hires_module.py
+ 39
− 17
View file @ 0954b891
......@@ -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
......
read_data.py
+ 6
− 2
View file @ 0954b891
......@@ -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)
......
restoral.py
+ 1
− 1
View file @ 0954b891
......@@ -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'),
......
scene.py
+ 6
− 6
View file @ 0954b891
......@@ -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'
......
sort_inputs.py
+ 53
− 14
View file @ 0954b891
......@@ -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()
spectral_tests.py
+ 8
− 6
View file @ 0954b891
......@@ -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) &
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment