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

added more spectral tests. Selected new scene and performed first combined test

parent 08600616
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
preprocess_thresholds.py
+ 53
− 90
View file @ 64486b66
......@@ -42,15 +42,14 @@ def thresholds_11_12um(data: xr.Dataset,
thresholds: np.ndarray,
scene: str,
scene_idx: np.ndarray) -> np.ndarray:
cosvza = np.cos(data.sensor_zenith[scene_idx].values.ravel() * _DTR)
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.array(schi, dtype=np.float32) # this is because the C function expects a float
m15 = data.M15.values[scene_idx].ravel()
latitude = data.latitude[scene_idx].values.ravel
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)
......@@ -128,46 +127,6 @@ def thresholds_NIR(data, thresholds, scene, test_name, scene_idx):
return corr_thr
def nir_refl(data, thresholds, scene_name):
sza = data.solar_zenith.values
band_n = 2
thr_nir = thresholds[scene_name]['1.6_2.1um_NIR_Reflectance_Test']
coeffs = thr_nir['coeffs']
vzcpow = thresholds['VZA_correction']['vzcpow'][0]
refang = data.sunglint_angle.values.reshape(np.prod(data.sunglint_angle.shape))
sunglint_thresholds = thresholds['Sun_Glint']
sunglint_flag = utils.sunglint_scene(refang, sunglint_thresholds).reshape(refang.shape)
hicut0 = np.array(coeffs[0] + coeffs[1]*sza + coeffs[2]*np.power(sza, 2) +
coeffs[3]*np.power(sza, 3) + coeffs[4]*np.power(sza, 4))
hicut0 = (hicut0*0.01) + thr_nir['adj']
hicut0 = (hicut0*thr_nir['bias']).reshape(refang.shape)
midpt0 = hicut0 + (thr_nir['midpt_coeff']*thr_nir['bias'])
locut0 = midpt0 + (thr_nir['locut_coeff']*thr_nir['bias'])
thr = np.array([locut0, midpt0, hicut0, np.ones(refang.shape)])
cosvza = np.cos(data.sensor_zenith*_dtr).values.reshape(refang.shape)
corr_thr = np.zeros((4, refang.shape[0]))
corr_thr[:3, sunglint_flag == 0] = thr[:3, sunglint_flag == 0] * (1./np.power(cosvza[sunglint_flag == 0],
vzcpow))
corr_thr[3, sunglint_flag == 0] = thr[3, sunglint_flag == 0]
for flag in range(1, 4):
if len(refang[sunglint_flag == flag]) > 0:
dosgref = utils.get_sunglint_thresholds(refang, sunglint_thresholds, band_n, flag, thr)
corr_thr[:3, sunglint_flag == flag] = dosgref[:3, sunglint_flag == flag] * \
(1./np.power(cosvza[sunglint_flag == flag], vzcpow))
corr_thr[3, sunglint_flag == flag] = dosgref[3, sunglint_flag == flag]
corr_thr = np.transpose(corr_thr.reshape((4, sza.shape[0], sza.shape[1])), (1, 2, 0))
return corr_thr
def preproc_surf_temp(data, thresholds):
thr_sfc1 = thresholds['Surface_Temperature_Test_1']
thr_sfc2 = thresholds['Surface_Temperature_Test_2']
......@@ -230,9 +189,9 @@ def var_11um(data, thresholds):
return var
def get_b1_thresholds(data, thresholds):
ndvi = data.ndvi.values.reshape(data.ndvi.shape[0]*data.ndvi.shape[1])
sctang = data.scattering_angle.values.reshape(data.ndvi.shape[0]*data.ndvi.shape[1])
def get_b1_thresholds(data, thresholds, scene_idx):
ndvi = data.ndvi.values[scene_idx].ravel()
sctang = data.scattering_angle.values[scene_idx].ravel()
# this is hardcoded in the function
delta_ndvi_bin = 0.1
......@@ -425,14 +384,14 @@ def get_nl_thresholds(data, threshold):
return out_thr
def vis_refl_thresholds(data, thresholds, scene):
def vis_refl_thresholds(data, thresholds, scene, scene_idx):
locut, midpt, hicut = get_b1_thresholds(data, thresholds)
locut, midpt, hicut = get_b1_thresholds(data, thresholds, scene_idx)
bias_adj = thresholds[scene]['Visible_Reflectance_Test']['adj']
ndvi = data.ndvi.values.reshape(data.ndvi.shape[0]*data.ndvi.shape[1])
m01 = data.M05.values.reshape(data.ndvi.shape[0]*data.ndvi.shape[1])
m02 = data.M07.values.reshape(data.ndvi.shape[0]*data.ndvi.shape[1])
m08 = data.M01.values.reshape(data.ndvi.shape[0]*data.ndvi.shape[1])
ndvi = data.ndvi.values[scene_idx].ravel()
m01 = data.M05.values[scene_idx].ravel()
m02 = data.M07.values[scene_idx].ravel()
m08 = data.M01.values[scene_idx].ravel()
m128 = m01
......@@ -458,38 +417,40 @@ def vis_refl_thresholds(data, thresholds, scene):
b1_power[idx] = thresholds[scene]['Visible_Reflectance_Test']['thr'][3]
m128[idx] = m02[idx]
cosvza = np.cos(data.sensor_zenith.values * _dtr).reshape(ndvi.shape)
cosvza = np.cos(data.sensor_zenith.values[scene_idx] * _DTR).ravel()
vzcpow = thresholds['VZA_correction']['vzcpow'][0]
b1_locut = (b1_locut * (1.0 / np.power(cosvza, vzcpow))).reshape(data.ndvi.shape)
b1_midpt = (b1_midpt * (1.0 / np.power(cosvza, vzcpow))).reshape(data.ndvi.shape)
b1_hicut = (b1_hicut * (1.0 / np.power(cosvza, vzcpow))).reshape(data.ndvi.shape)
out_thr = xr.DataArray(data=np.dstack((b1_locut, b1_midpt, b1_hicut, np.ones(data.ndvi.shape),
b1_power.reshape(data.ndvi.shape))),
dims=('number_of_lines', 'number_of_pixels', 'z'))
out_rad = xr.DataArray(data=m128.reshape(data.M01.shape), dims=('number_of_lines', 'number_of_pixels'))
b1_locut = (b1_locut * (1.0 / np.power(cosvza, vzcpow)))
b1_midpt = (b1_midpt * (1.0 / np.power(cosvza, vzcpow)))
b1_hicut = (b1_hicut * (1.0 / np.power(cosvza, vzcpow)))
return out_thr, out_rad
# out_thr = xr.DataArray(data=np.dstack((b1_locut, b1_midpt, b1_hicut, np.ones(data.ndvi.shape),
# b1_power.reshape(data.ndvi.shape))),
# dims=('number_of_lines', 'number_of_pixels', 'z'))
# out_rad = xr.DataArray(data=m128.reshape(data.M01.shape), dims=('number_of_lines', 'number_of_pixels'))
out_thr = np.dstack((b1_locut, b1_midpt, b1_hicut, b1_power))
out_rad = m128
return np.squeeze(out_thr.T), out_rad
def GEMI_test(data, thresholds, scene_name):
def GEMI_thresholds(data, thresholds, scene_name, scene_idx):
thresh = thresholds[scene_name]['GEMI_Test']
ndvi = data.ndvi.values[scene_idx].ravel()
gemi_thr = np.ones((ndvi.shape[0], 5))
gemi_thr = np.ones((data.M01.shape[0], data.M01.shape[1], 5))
idx = np.nonzero(ndvi < 0.1)
gemi_thr[idx, :3] = thresh['gemi0'][:3]
idx = np.nonzero((ndvi >= 0.1) & (data.ndvi < 0.2))
gemi_thr[idx, :3] = thresh['gemi1'][:3]
idx = np.nonzero((ndvi >= 0.2) & (ndvi < 0.3))
gemi_thr[idx, :3] = thresh['gemi2'][:3]
idx = np.nonzero(data.ndvi < 0.1)
gemi_thr[idx[0], idx[1], :3] = thresh['gemi0'][:3]
idx = np.nonzero((data.ndvi >= 0.1) & (data.ndvi < 0.2))
gemi_thr[idx[0], idx[1], :3] = thresh['gemi1'][:3]
idx = np.nonzero((data.ndvi >= 0.2) & (data.ndvi < 0.3))
gemi_thr[idx[0], idx[1], :3] = thresh['gemi2'][:3]
# thr_out = xr.DataArray(data=np.dstack((gemi_thr[:, :, 0], gemi_thr[:, :, 1], gemi_thr[:, :, 2],
# np.ones(gemi_thr[:, :, 0].shape),
# np.ones(gemi_thr[:, :, 0].shape))),
# dims=('number_of_lines', 'number_of_pixels', 'z'))
thr_out = xr.DataArray(data=np.dstack((gemi_thr[:, :, 0], gemi_thr[:, :, 1], gemi_thr[:, :, 2],
np.ones(gemi_thr[:, :, 0].shape),
np.ones(gemi_thr[:, :, 0].shape))),
dims=('number_of_lines', 'number_of_pixels', 'z'))
return thr_out
return gemi_thr
def bt11_4um_preproc(data, thresholds, scene_name):
......@@ -510,29 +471,31 @@ def bt11_4um_preproc(data, thresholds, scene_name):
return thr_out
def test_1_38um_preproc(data, thresholds, scene_name):
sza = data.solar_zenith.values
vza = data.sensor_zenith.values
def thresholds_1_38um_test(data, thresholds, scene_name, scene_idx):
sza = data.solar_zenith.values[scene_idx]
vza = data.sensor_zenith.values[scene_idx]
thresh = thresholds[scene_name]['1.38um_High_Cloud_Test']
c = thresh['coeffs']
vzcpow = thresholds['VZA_correction']['vzcpow'][1]
hicut0 = c[0] + c[1]*sza + c[2]*np.power(sza, 2) + c[3]*np.power(sza, 3) + \
hicut = c[0] + c[1]*sza + c[2]*np.power(sza, 2) + c[3]*np.power(sza, 3) + \
c[4]*np.power(sza, 4) + c[5]*np.power(sza, 5)
hicut0 = hicut0*0.01 + (np.maximum((sza/90.)*thresh['szafac']*thresh['adj'], thresh['adj']))
midpt0 = hicut0 + 0.001
locut0 = midpt0 + 0.001
hicut = hicut*0.01 + (np.maximum((sza/90.)*thresh['szafac']*thresh['adj'], thresh['adj']))
midpt = hicut + 0.001
locut = midpt + 0.001
cosvza = np.cos(vza*_dtr)
cosvza = np.cos(vza*_DTR)
locut0 = locut0 * (1/np.power(cosvza, vzcpow))
midpt0 = midpt0 * (1/np.power(cosvza, vzcpow))
hicut0 = hicut0 * (1/np.power(cosvza, vzcpow))
locut = locut * (1/np.power(cosvza, vzcpow))
midpt = midpt * (1/np.power(cosvza, vzcpow))
hicut = hicut * (1/np.power(cosvza, vzcpow))
out_thr = xr.DataArray(data=np.dstack((locut0, midpt0, hicut0, np.ones(data.ndvi.shape),
np.ones(data.ndvi.shape))),
dims=('number_of_lines', 'number_of_pixels', 'z'))
return out_thr
# out_thr = xr.DataArray(data=np.dstack((locut, midpt, hicut, np.ones(data.ndvi.shape),
# np.ones(data.ndvi.shape))),
# dims=('number_of_lines', 'number_of_pixels', 'z'))
out_thr = np.dstack((locut, midpt, hicut, np.ones(locut.shape), np.ones(locut.shape)))
return np.squeeze(out_thr.T)
# NOTE: 11-12um Cirrus Test
# hicut is computed in different ways depending on the scene
......
read_data.py
+ 8
− 0
View file @ 64486b66
......@@ -148,6 +148,13 @@ def get_data(file_names: Dict[str, str],
viirs_data = read_data('viirs', f'{mod02}', f'{mod03}')
viirs_data = read_ancillary_data(file_names, viirs_data)
m01 = viirs_data.M01.values
m02 = viirs_data.M02.values
r1 = 2.0 * (np.power(m02, 2.0) - np.power(m01, 2.0)) + (1.5 * m02) + (0.5 * m01)
r2 = m02 + m01 + 0.5
r3 = r1 / r2
gemi = r3 * (1.0 - 0.25*r3) - ((m01 - 0.125) / (1.0 - m01))
# Compute channel differences and ratios that are used in the tests
viirs_data['M14-M15'] = (('number_of_lines', 'number_of_pixels'),
viirs_data.M14.values - viirs_data.M15.values)
......@@ -159,6 +166,7 @@ def get_data(file_names: Dict[str, str],
viirs_data.M13.values - viirs_data.M16.values)
viirs_data['M07-M05ratio'] = (('number_of_lines', 'number_of_pixels'),
viirs_data.M07.values / viirs_data.M05.values)
viirs_data['GEMI'] = (('number_of_lines', 'number_of_pixels'), gemi)
scene_flags = scn.find_scene(viirs_data, sunglint_angle)
scene = scn.scene_id(scene_flags)
......
scene.py
+ 2
− 2
View file @ 64486b66
......@@ -11,7 +11,7 @@ import ancillary_data as anc
# lsf: land sea flag
_scene_list = ['Ocean_Day', 'Ocean_Night', 'Land_Day', 'Land_Night', 'Day_Snow', 'Night_Snow', 'Coast_Day',
'Desert_Day', 'Antarctic_Day', 'Polar_Day_Snow', 'Polar_Day_Desert', 'Polar_Day_Ocean',
'Land_Day_Desert', 'Antarctic_Day', 'Polar_Day_Snow', 'Polar_Day_Desert', 'Polar_Day_Ocean',
'Polar_Day_Desert_Coast', 'Polar_Day_Coast', 'Polar_Day_Land', 'Polar_Night_Snow',
'Polar_Night_Land', 'Polar_Night_Ocean', 'Land_Day_Desert_Coast', 'Land_Day_Coast', 'Desert']
_flags = ['day', 'night', 'land', 'coast', 'sh_lake', 'sh_ocean', 'water', 'polar', 'sunglint',
......@@ -364,7 +364,7 @@ def scene_id(scene_flag):
idx = np.nonzero((scene_flag['land'] == 1) & (scene_flag['day'] == 1) &
(scene_flag['desert'] == 1) & (scene_flag['coast'] == 0) &
(scene_flag['polar'] == 0) & (scene_flag['antarctica'] == 0))
scene['Desert_Day'][idx] = 1
scene['Land_Day_Desert'][idx] = 1
# Land Day Desert Coast
idx = np.nonzero((scene_flag['land'] == 1) & (scene_flag['day'] == 1) &
......
thresholds.mvcm.snpp.v0.0.1.yaml
+ 39
− 11
View file @ 64486b66
......@@ -16,7 +16,10 @@ Land_Day:
Visible_Reflectance_Test:
thr: [0.207, 0.169, 0.132, 1.0, 1.0] # used to be dlref1
adj: 0.94 # used to be dl_b1bias_adj
1.38um_High_Cloud_Test: [0.0375, 0.0250, 0.0125, 1.0, 1.0]
perform: True
1.38um_High_Cloud_Test:
thr: [0.0375, 0.0250, 0.0125, 1.0, 1.0]
perform: True
dlvrat : [1.80, 1.85, 1.90, 1.0, 1.0]
# dl11_12lcmult : 0.3
# dl11_12hcadj : 1.25
......@@ -68,7 +71,10 @@ Land_Day_Coast:
Visible_Reflectance_Test:
thr: [0.207, 0.169, 0.132, 1.0, 1.0]
adj: 0.94
1.38um_High_Cloud_Test: [0.0375, 0.0250, 0.0125, 1.0, 1.0]
perform: True
1.38um_High_Cloud_Test:
thr: [0.0375, 0.0250, 0.0125, 1.0, 1.0]
perform: True
# dl11_12lcmult_t2 : 0.3
# dl11_12hcadj_t2 : 1.25
......@@ -94,7 +100,10 @@ Land_Day_Desert:
thr: [0.326, 0.288, 0.250, 1.00, 1.0]
adj: 0.94
ndvi_thr: 0.25 # this used to be lds_ndvi
1.38um_High_Cloud_Test: [0.0375, 0.0250, 0.0125, 1.00, 1.0]
perform: True
1.38um_High_Cloud_Test:
thr: [0.0375, 0.0250, 0.0125, 1.0, 1.0]
perform: True
GEMI_Test:
gemi0: [0.085, 0.095, 0.115, 1.00, 1.0]
gemi1: [0.145, 0.170, 0.220, 1.00]
......@@ -115,10 +124,13 @@ Land_Day_Desert_Coast:
CO2_High_Clouds_Test: [222.0, 224.0, 226.0, 1.0, 1.0]
Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
Visible_Reflectance_Test:
thr: [0.326, 0.288, 0.250, 1.00, 1.0]
thr: [0.326, 0.288, 0.250, 1.0, 1.0]
adj: 0.94
ndvi_thr: 0.25
1.38um_High_Cloud_Test: [0.0375, 0.0250, 0.0125, 1.00, 1.0]
perform: True
1.38um_High_Cloud_Test:
thr: [0.0375, 0.0250, 0.0125, 1.0, 1.0]
perform: True
lds_ref3_tpw_c : -10.00
ldsbt1_c : 320.0
# lds11_12lcmult_c : 0.3
......@@ -218,7 +230,10 @@ Polar_Day_Land:
Visible_Reflectance_Test:
thr: [0.207, 0.169, 0.132, 1.0, 1.0]
adj: 0.94
1.38um_High_Cloud_Test: [0.0375, 0.0250, 0.0125, 1.0, 1.0]
perform: True
1.38um_High_Cloud_Test:
thr: [0.0375, 0.0250, 0.0125, 1.0, 1.0]
perform: True
Polar_Night_Land:
Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
......@@ -266,7 +281,10 @@ Polar_Day_Coast:
Visible_Reflectance_Test:
thr: [0.207, 0.169, 0.132, 1.0, 1.0]
adj: 0.94
1.38um_High_Cloud_Test: [0.0375, 0.0250, 0.0125, 1.0, 1.0]
perform: True
1.38um_High_Cloud_Test:
thr: [0.0375, 0.0250, 0.0125, 1.0, 1.0]
perform: True
Polar_Day_Desert:
11-12um_Cirrus_Test:
......@@ -279,7 +297,10 @@ Polar_Day_Desert:
Visible_Reflectance_Test:
thr: [0.326, 0.288, 0.250, 1.00, 1.0]
adj: 0.94
1.38um_High_Cloud_Test: [0.0375, 0.0250, 0.0125, 1.00, 1.0]
perform: True
1.38um_High_Cloud_Test:
thr: [0.0375, 0.0250, 0.0125, 1.00, 1.0]
perform: True
GEMI_Test:
gemi0: [0.085, 0.110, 0.135, 1.00, 1.0]
gemi1: [0.170, 0.220, 0.270, 1.00]
......@@ -298,7 +319,10 @@ Polar_Day_Desert_Coast:
Visible_Reflectance_Test:
thr: [0.326, 0.288, 0.250, 1.00, 1.0]
adj: 0.94
1.38um_High_Cloud_Test: [0.0375, 0.0250, 0.0125, 1.00, 1.0]
perform: True
1.38um_High_Cloud_Test:
thr: [0.0375, 0.0250, 0.0125, 1.00, 1.0]
perform: True
pds_ref3_tpw_c : -10.00
pdsbt1_c : 320.0
......@@ -311,7 +335,9 @@ Polar_Day_Snow:
dps_ref3_tpw : 0.75
dpsbt1 : 230.0
Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
1.38um_High_Cloud_Test: [0.060, 0.0525, 0.045, 1.0, 1.0]
1.38um_High_Cloud_Test:
thr: [0.060, 0.0525, 0.045, 1.0, 1.0]
perform: True
dps11_12adj : 0.8
dps4_11l : [20.0, 18.0, 16.0, 1.0]
dps4_11h : [18.0, 16.0, 14.0, 1.0]
......@@ -450,7 +476,9 @@ Day_Snow:
ds_ref3_tpw : 0.75
CO2_High_Clouds_Test: [222.0, 224.0, 226.0, 1.0, 1.0]
Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
1.38um_High_Cloud_Test: [0.060, 0.0525, 0.045, 1.0, 1.0]
1.38um_High_Cloud_Test:
thr: [0.060, 0.0525, 0.045, 1.0, 1.0]
perform: True
ds11_12adj : 0.8
ds11_12lcmult : 0.3
ds11_12hcmult : 0.3
......
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