-
Paolo Veglio authored11-12um Thin Cirrus Test has been implemented for all scene. Results still need to be verified, but look reasonable
Paolo Veglio authored11-12um Thin Cirrus Test has been implemented for all scene. Results still need to be verified, but look reasonable
preprocess_thresholds.py 8.29 KiB
import numpy as np
import xarray as xr
import ancillary_data as anc
import utils
_dtr = np.pi/180
# this case is written for the 11-12um Cirrus Test for scenes that follow pattern 1 (see note below)
def prepare_thresholds(data, thresholds):
coeff_values = np.empty((data.M01.shape[0], data.M01.shape[1], 2))
coeff_values[:, :, 0] = np.full(data.M01.shape, thresholds['11-12um_Cirrus_Test']['coeffs'][0])
coeff_values[:, :, 1] = np.full(data.M01.shape, thresholds['11-12um_Cirrus_Test']['coeffs'][1])
cmult_values = np.full(data.M01.shape, thresholds['11-12um_Cirrus_Test']['cmult'])
adj_values = np.full(data.M01.shape, thresholds['11-12um_Cirrus_Test']['adj'])
if 'bt1' in list(thresholds['11-12um_Cirrus_Test']):
bt1 = np.full(data.M01.shape, thresholds['11-12um_Cirrus_Test']['bt1'])
else:
bt1 = np.full(data.M01.shape, -999)
if 'lat' in list(thresholds['11-12um_Cirrus_Test']):
lat = np.full(data.M01.shape, thresholds['11-12um_Cirrus_Test']['lat'])
else:
lat = np.full(data.M01.shape, -999)
thr_dict = {'coeffs': (['number_of_lines', 'number_of_pixels', 'z'], coeff_values),
'cmult': (['number_of_lines', 'number_of_pixels'], cmult_values),
'adj': (['number_of_lines', 'number_of_pixels'], adj_values),
'bt1': (['number_of_lines', 'number_of_pixels'], bt1),
'lat': (['number_of_lines', 'number_of_pixels'], lat),
}
return xr.Dataset(data_vars=thr_dict)
# 3. Others
# - Land_Night
# - Polar_Night_Land
# - Polar_Night_Snow
# - Day_Snow
# - Night_Snow
def preproc(data, thresholds, scene):
cosvza = np.cos(data.sensor_zenith * _dtr)
schi = (1/cosvza).where(cosvza > 0, 99.0)
schi = schi.values.reshape(np.prod(schi.shape))
m15 = data.M15.values.reshape(np.prod(data.M15.shape))
thr = anc.py_cithr(1, schi, m15)
thr = thr.reshape(data.M15.shape)
schi = schi.reshape(data.M15.shape)
# thr_xr = xr.Dataset(np.full(data.sensor_zenith.shape, thresholds['coeffs']),
# dims=('number_of_lines', 'number_of_pixels'))
thr_xr = prepare_thresholds(data, thresholds)
midpt = thr_xr.coeffs[:, :, 0].where((thr < 0.1) | (np.abs(schi-99) < 0.0001), thr)
locut = midpt + (thr_xr.cmult * midpt)
if scene in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert', 'Land_Day_Desert_Coast',
'Ocean_Day', 'Ocean_Night', 'Polar_Ocean_Day', 'Polar_Ocean_Night']:
hicut = midpt - thr_xr.adj
elif scene in ['Polar_Day_Land', 'Polar_Day_Coast', 'Polar_Day_Desert',
'Polar_Day_Desert_Coast', 'Polar_Day_Snow']:
hicut = midpt - (thr_xr.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}
midpt = midpt - (_coeffs[scene] * locut)
if scene in ['Polar_Night_Land', 'Polar_Night_Snow', 'Night_Snow']:
hicut = (midpt - (0.2 * locut)).where(data.M15 < thr_xr.bt1, midpt - 1.25)
elif scene in ['Land_Night']:
hicut = -0.1 - np.power(90.0 - np.abs(data.latitude)/60, 4) * 1.15
hicut = hicut.where((data.M15 < thr_xr.bt1) & (data.latitude > thr_xr.lat), 1.25)
elif scene in ['Day_Snow']:
hicut = locut - (thr_xr.cmult * locut)
else:
print('Scene not recognized\n')
else:
print('Scene not recognized\n')
thr_out = xr.DataArray(data=np.dstack((locut, midpt, hicut, np.ones(locut.shape), np.ones(locut.shape))),
dims=('number_of_lines', 'number_of_pixels', 'z'))
return thr_out
# return locut, hicut, midpt
def preproc_sst(data, thresholds):
m31c = data.M15 - 273.16
m32c = data.M16 - 273.16
m31c_m32c = m31c - m32c
sstc = data.geos_sfct - 273.16
cosvza = np.cos(data.sensor_zenith*_dtr)
a = thresholds['coeffs']
modsst = 273.16 + a[0] + a[1]*m31c + a[2]*m31c_m32c*sstc + a[3]*m31c_m32c*((1/cosvza) - 1)
sfcdif = data.geos_sfct - modsst
return sfcdif
def preproc_nir(data, thresholds, scene):
sza = data.solar_zenith
band_n = 2
# NOTE: the visud condition in the C code is equivalent to having sza <= 85
# For the time being the visud filtering is not implemented
a = np.array(thresholds[scene]['NIR_Reflectance_Test']['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)
nir_thresh = thresholds[scene]['NIR_Reflectance_Test']
hicut0 = a[0] + a[1]*sza + a[2]*np.power(sza, 2) + a[3]*np.power(sza, 3)
hicut0 = (hicut0 * 0.01) + nir_thresh['adj']
hicut0 = (hicut0 + nir_thresh['bias']).values.reshape(refang.shape)
midpt0 = hicut0 + (nir_thresh['midpt_coeff'] * nir_thresh['bias'])
locut0 = midpt0 + (nir_thresh['locut_coeff'] * nir_thresh['bias'])
thr = np.array([locut0, midpt0, hicut0, nir_thresh['thr'][3]*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']
thr_df1 = thresholds['Surface_Temperature_Test_df1']
thr_df2 = thresholds['Surface_Temperature_Test_df2']
max_vza = 70.13 # This values is set based on sensor. Check mask_processing_constants.h for MODIS value
rs = np.prod(data.M15.shape)
df1 = (data.M15 - data.M16).values.reshape(rs)
df2 = (data.M15 - data.M13).values.reshape(rs)
desert_flag = data.Desert.values.reshape(rs)
thresh = np.ones((rs, )) * thr_sfc1
idx = np.where((df1 >= thr_df1[0]) | ((df1 < thr_df1[0]) & ((df2 <= thr_df2[0]) | (df2 >= thr_df2[1]))))
thresh[idx] = thr_sfc2
idx = np.where(desert_flag == 1)
thresh[idx] == thr_sfc1
midpt = thresh
idx = np.where(df1 >= thr_df1[1])
midpt[idx] = thresh[idx] + 2.0*df1[idx]
corr = np.power(data.sensor_zenith.values/max_vza, 4) * 3.0
midpt = midpt.reshape(corr.shape) + corr
locut = midpt + 2.0
hicut = midpt - 2.0
thr_out = xr.DataArray(data=np.dstack((locut, midpt, hicut, np.ones(locut.shape), np.ones(locut.shape))),
dims=('number_of_lines', 'number_of_pixels', 'z'))
return thr_out
def get_b1_thresholds():
# fill_ndvi[0] is for fill_ndvi1
# fill_ndvi[1] is for fill_ndvi2
pass
# NOTE: 11-12um Cirrus Test
# hicut is computed in different ways depending on the scene
# 1. midpt - adj
# - Land_Day
# - Land_Day_Coast
# - Land_Day_Desert
# - Land_Day_Desert_Coast
# - Ocean_Day
# - Ocean_Night
# - Polar_Day_Ocean
# - Polar_Night_Ocean
#
# 2. midpt - (btd_thr * adj)
# - Polar_Day_Land
# - Polar_Day_Coast
# - Polar_Day_Desert
# - Polar_Day_Desert_Coast
# - Polar_Day_Snow
#
# 3. Others
# - Land_Night
# - Polar_Night_Land
# - Polar_Night_Snow
# - Day_Snow
# - Night_Snow
# NOTE: 1.38um High Cloud Test
# thresholds are not always computed the same way. In group 1 there's no preprocessing required,
# in group 2 some calcuations are needed
# 1.
# - Land_Day
# - Land_Day_Coast
# - Land_Day_Desert
# - Land_Day_Desert_Coast
# - Polar_Day_Land
# - Polar_Day_Coast
# - Polar_Day_Desert
# - Polar_Day_Desert_Coast
# - Polar_Day_Snow
# - Day_Snow
#
# 2.
# - Ocean_Day
# - Polar_Ocean_Day