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

started implementation of sunglint calculations.

parent f044f62a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
main.py
+ 27
− 8
View file @ 8f867ce8
import tests
import ruamel_yaml as yml
import numpy as np
import read_data as rd
import tests
def main():
datapath = '/ships19/hercules/pveglio/neige_data/snpp_test_input'
fname_l1b = 'VNP02MOD.A2014213.1548.001.2017301015346.uwssec.bowtie_restored_scaled.nc'
fname_geo = 'VNP03MOD.A2014213.1548.001.2017301015705.uwssec.nc'
thresh_file = '/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.yaml'
viirs_data = rd.read_data(f'{datapath}/{fname_l1b}', f'{datapath}/{fname_geo}')
with open(thresh_file) as f:
text = f.read()
thresholds = yml.safe_load(text)
confidence = np.zeros(2, viirs_data['M01'].shape[0], viirs_data['M01'].shape[1])
confidence[0, :, :] = tests.test_11um(viirs_data.M15.values, thresholds['Daytime_Ocean'])
confidence[1, :, :] = tests.test_11_4diff(viirs_data.M15.values, viirs_data.M13.values,
thresholds['Daytime_Ocean'])
return confidence
def test_main():
rad1 = [[255, 260, 265, 248, 223],
[278, 285, 270, 268, 256],
......@@ -13,7 +36,6 @@ def main():
[280, 281, 272, 270, 267]]
thresh_file = '/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.yaml'
with open(thresh_file) as f:
text = f.read()
......@@ -23,12 +45,9 @@ def main():
confidence = np.zeros((2, rad1.shape[0], rad1.shape[1]))
thresholds = yml.safe_load(text)
thr_11um = thresholds['Daytime_Ocean']
thr_11_4diff = thresholds['Daytime_Ocean']
confidence[0, :, :] = tests.test_11um(rad1, thr_11um['bt11'])
confidence[1, :, :] = tests.test_11_4diff(rad1, rad2, thr_11_4diff['test11_4lo'])
confidence[0, :, :] = tests.test_11um(rad1, thresholds['Daytime_Ocean'])
confidence[1, :, :] = tests.test_11_4diff(rad1, rad2, thresholds['Daytime_Ocean'])
print(f'Confidence[0,:,:]: \n {confidence[0, :, :]}')
print(f'Confidence[1,:,:]: \n {confidence[1, :, :]}')
......@@ -37,4 +56,4 @@ def main():
if __name__ == "__main__":
main()
test_main()
read_data.py
+ 24
− 1
View file @ 8f867ce8
from netCDF4 import Dataset
# from netCDF4 import Dataset
import xarray as xr
import numpy as np
_DTR = np.pi/180.
_RTD = 180./np.pi
def read_data(sensor: str, l1b_filename: str, geo_filename: str):
......@@ -28,6 +31,26 @@ def read_data(sensor: str, l1b_filename: str, geo_filename: str):
in_data = in_data.merge(data)
relazi = relative_azimuth_angle(data['sensor_azimuth'].values, data['solar_azimuth'].values)
sunglint = sun_glint_angle(data['sensor_zenith'].values, data['solar_zenith'].values, relazi)
in_data['relative_azimuth'] = (('number_of_lines', 'number_of_pixels'), relazi)
in_data['sunglint_angle'] = (('number_of_lines', 'number_of_pixels'), sunglint)
return in_data
def relative_azimuth_angle(sensor_azimuth: float, solar_azimuth: float) -> float:
rel_azimuth = np.abs(180. - np.abs(sensor_azimuth - solar_azimuth))
return rel_azimuth
def sun_glint_angle(sensor_zenith: float, solar_zenith: float, rel_azimuth: float) -> float:
cossna = (np.sin(sensor_zenith*_DTR) * np.sin(solar_zenith*_DTR) * np.cos(rel_azimuth*_DTR) +
np.cos(sensor_zenith*_DTR) * np.cos(solar_zenith*_DTR))
sunglint_angle = np.arccos(cossna) * _RTD
return sunglint_angle
def correct_reflectances():
pass
tests.py
+ 70
− 178
View file @ 8f867ce8
import numpy as np
_test_rad = np.random.randint(25, size=[6, 8])
#_test_thr = [15, 10, 5, 1, 1]
_test_thr = [5, 10, 15, 1, 1]
import utils
def test_11um(rad, threshold):
#if (~np.isnan(rad) or threshold[4] == 1.0):
confidence = conf_test(rad, threshold)
thr = threshold['bt11']
if thr[4] == 1:
# the C code has the line below that I don't quite understand the purpose of.
# It seems to be setting the bit to 0 if the BT value is greater than the midpoint
#
# if (m31 >= dobt11[1]) (void) set_bit(13, pxout.testbits);
confidence = utils.conf_test(rad, thr)
return confidence
def test_11_4diff(rad1, rad2, threshold):
raddiff = rad1 - rad2;
confidence = conf_test(raddiff, threshold)
day = True
sunglint = False
thr = threshold['test11_4lo']
raddiff = rad1 - rad2
if day is True and sunglint is False:
confidence = utils.conf_test(raddiff, thr)
return confidence
def nir_refl_test(rad, threshold):
sunglint = False
sza = 0
refang = 0
vza = 0
dtr = 0
band_n = 2
vzcpow = 0
sunglint_thr = np.zeros((4,))
# ref2 [5]
# b2coeffs [4]
# b2mid [1]
# b2bias_adj [1]
# b2lo [1]
# vzcpow [3] (in different place)
coeffs = threshold['b2coeffs']
hicut0 = coeffs[0] + coeffs[1]*sza + coeffs[2]*np.power(sza, 2) + coeffs[3]*np.power(sza, 3)
hicut0 = (hicut0 * 0.01) + threshold['b2adj']
hicut0 = hicut0 * threshold['b2bias_adj']
midpt0 = hicut0 + (threshold['b2mid'] * threshold['b2bias_adj'])
locut0 = midpt0 + (threshold['b2lo'] * threshold['b2bias_adj'])
if sunglint is True:
# Keep in mind that band_n uses MODIS band numbers (2=0.86um and 7=2.1um)
# For VIIRS would be 2=M7 (0.865um) and 7=M11 (2.25um)
sunglint_thr = utils.get_sunglint_thresholds(refang, threshold['Sun_Glint'],
band_n, sunglint_thr)
locut1 = sunglint[0]
midpt1 = sunglint[1]
hicut1 = sunglint[2]
else:
locut1 = locut0
midpt1 = midpt0
hicut1 = hicut0
cosvza = np.cos(vza*dtr)
locut2 = locut1 * (1./np.power(cosvza, vzcpow[0]))
midpt2 = midpt1 * (1./np.power(cosvza, vzcpow[0]))
hicut2 = hicut1 * (1./np.power(cosvza, vzcpow[0]))
corr_thr = [locut2, hicut2, 1.0, midpt2]
confidence = utils.conf_test(rad, corr_thr)
return confidence
def simple_threshold_test(rad, threshold):
return conf_test(rad, threshold)
def vis_nir_ratio_test(rad1, rad2, threshold):
pass
class CloudMaskTests(object):
......@@ -51,177 +108,12 @@ class CloudMaskTests(object):
pass
def conf_test(rad=_test_rad, thr=_test_thr):
'''
Assuming a linear function between min and max confidence level, the plot below shows
how the confidence (y axis) is computed as function of radiance (x axis).
This case illustrates alpha < gamma, obviously in case alpha > gamma, the plot would be
flipped.
gamma
c 1 ________
o | /
n | /
f | /
i | beta /
d 1/2 |....../
e | /
n | /
c | /
e 0________/
| alpha
--------- radiance ---------->
'''
coeff = np.power(2, (thr[3] - 1))
hicut = thr[0]
beta = thr[1]
locut = thr[2]
power = thr[3]
radshape = rad.shape
rad = rad.reshape((rad.shape[0]*rad.shape[1]))
c = np.zeros(rad.shape)
if hicut > locut:
gamma = thr[0]
alpha = thr[2]
flipped = False
else:
gamma = thr[2]
alpha = thr[0]
flipped = True
# Rad between alpha and beta
range_ = 2. * (beta - alpha)
s1 = (rad[rad <= beta] - alpha) / range_
if flipped is False:
c[rad <= beta] = coeff * np.power(s1, power)
if flipped is True:
c[rad <= beta] = 1. - (coeff * np.power(s1, power))
# Rad between beta and gamma
range_ = 2. * (beta - gamma)
s1 = (rad[rad > beta] - gamma) / range_
if flipped is False:
c[rad > beta] = 1. - (coeff * np.power(s1, power))
if flipped is True:
c[rad > beta] = coeff * np.power(s1, power)
# Rad outside alpha-gamma interval
if flipped is False:
c[rad > gamma] = 1
c[rad < alpha] = 0
if flipped is True:
c[rad > gamma] = 0
c[rad < alpha] = 1
c[c > 1] = 1
c[c < 0] = 0
confidence = c.reshape(radshape)
return confidence
def conf_test_dble(rad, coeffs):
'''
gamma1 gamma2
c 1_______ ________
o | \ /
n | \ /
f | \ /
i | \ beta1 beta2 /
d 1/2 \....| |...../
e | \ /
n | \ /
c | \ /
e 0 \_____________/
| alpha1 alpha2
--------------------- radiance ------------------------->
'''
hicut = [coeffs[0], coeffs[1]]
locut = [coeffs[2], coeffs[3]]
midpt = [coeffs[4], coeffs[5]]
power = coeffs[6]
gamma1 = hicut[0]
gamma2 = hicut[1]
alpha1 = locut[0]
alpha2 = locut[1]
beta1 = midpt[0]
beta2 = midpt[1]
coeff = np.power(2, (power - 1))
radshape = rad.shape
rad = rad.reshape((rad.shape[0]*rad.shape[1]))
c = np.zeros(rad.shape)
# Find if interval between inner cutoffs passes or fails test
if (alpha1 - gamma1 > 0):
# Value is within range of lower set of limits
range_ = 2 * (beta1 - alpha1)
s1 = (rad[(rad <= alpha1) & (rad >= beta1)] - alpha1) / range_
c[(rad <= alpha1) & (rad >= beta1)] = coeff * np.power(s1, power)
range_ = 2 * (beta1 - gamma1)
s1 = (rad[(rad <= alpha1) & (rad < beta1)] - gamma1) / range_
c[(rad <= alpha1) & (rad < beta1)] = coeff * np.power(s1, power)
# Value is within range of upper set of limits
range_ = 2 * (beta2 - alpha2)
s1 = (rad[(rad > alpha1) & (rad <= beta2)] - alpha2) / range_
c[(rad > alpha1) & (rad <= beta2)] = coeff * np.power(s1, power)
range_ = 2 * (beta2 - gamma2)
s1 = (rad[(rad > alpha1) & (rad > beta2)] - gamma2) / range_
c[(rad > alpha1) & (rad > beta2)] = coeff * np.power(s1, power)
# Inner region fails test
# Check for value beyond function range
c[(rad > alpha1) & (rad < alpha2)] = 0
c[(rad < gamma1) | (rad > gamma2)] = 1
else:
# Value is withing range of lower set of limits
range_ = 2 * (beta1 - alpha1)
s1 = (rad[(rad <= gamma1) & (rad <= beta1)] - alpha1) / range_
c[(rad <= gamma1) & (rad <= beta1)] = coeff * np.power(s1, power)
range_ = 2 * (beta1 - gamma1)
s1 = (rad[(rad <= gamma1) & (rad > beta1)] - gamma1) / range_
c[(rad <= gamma1) & (rad > beta1)] = coeff * np.power(s1, power)
# Value is within range of upper set of limits
range_ = 2 * (beta2 - alpha2)
s1 = (rad[(rad > gamma1) & (rad >= beta2)] - alpha2) / range_
c[(rad > gamma1) & (rad >= beta2)] = coeff * np.power(s1, power)
range_ = 2 * (beta2 - gamma2)
s1 = (rad[(rad > gamma1) & (rad < beta2)] - gamma2) / range_
c[(rad > gamma1) & (rad < beta2)] = coeff * np.power(s1, power)
# Inner region passes test
# Check for value beyond function range
c[(rad > gamma1) & (rad < gamma2)] = 1
c[(rad < alpha1) | (rad > alpha2)] = 0
c[c>1] = 1
c[c<0] = 0
confidence = c.reshape(radshape)
return confidence
def test():
rad = np.random.randint(50, size=[4, 8])
#coeffs = [5, 42, 20, 28, 15, 35, 1]
#coeffs = [20, 28, 5, 42, 15, 35, 1]
# coeffs = [5, 42, 20, 28, 15, 35, 1]
# coeffs = [20, 28, 5, 42, 15, 35, 1]
coeffs = [35, 15, 20, 1, 1]
#confidence = conf_test_dble(rad, coeffs)
# confidence = conf_test_dble(rad, coeffs)
confidence = test_11um(rad, coeffs)
print(rad)
print('\n')
......
utils.py 0 → 100644
+ 225
− 0
View file @ 8f867ce8
import numpy as np
_test_rad = np.random.randint(25, size=[6, 8])
# _test_thr = [15, 10, 5, 1, 1]
_test_thr = [5, 10, 15, 1, 1]
# this function creates a map of sunglint areas, based on the different angles set in the
# threshold file. The goal is to create an array of indices that I can use to quickly assign
# different coefficients depending on the angle interval. This will be mostly used in the
# function get_sunglint_thresholds().
# All of this is because we want to be able to process the whole array, instead of iterating
# over all pixels one by one.
def sunglint_scene(refang, sunglint_thr):
sunglint_flag = np.zeros(refang.shape)
sunglint_flag[refang <= sunglint_thr['bounds'][3]] = 1
sunglint_flag[refang <= sunglint_thr['bounds'][2]] = 2
sunglint_flag[refang <= sunglint_thr['bounds'][1]] = 3
return sunglint_flag
def get_sunglint_thresholds(refang, thresholds, band_n, sunglint):
band = f'band{band_n}'
# if refang > thresholds['bounds'][3]:
# sunglint = sunglint
# # dosgref[2] = hicnf
# # dosgref[0] = locnf
# # dosgref[1] = mdcnf
# # sunglint[3] = doref2[3]
if refang <= thresholds['bounds'][1]:
sunglint = thresholds[f'{band}_0deg']
else:
if (refang > thresholds['bounds'][1] and refang <= thresholds['bounds'][2]):
lo_ang = thresholds['bounds'][1]
hi_ang = thresholds['bounds'][2]
lo_ang_val = thresholds[f'{band}_10deg'][0]
hi_ang_val = thresholds[f'{band}_10deg'][1]
power = thresholds[f'{band}_10deg'][3]
conf_range = thresholds[f'{band}_10deg'][2]
elif (refang > thresholds['bounds'][2] and refang <= thresholds['bounds'][3]):
lo_ang = thresholds['bounds'][2]
hi_ang = thresholds['bounds'][3]
lo_ang_val = thresholds[f'{band}_20deg'][0]
hi_ang_val = sunglint[1]
power = thresholds[f'{band}_20deg'][3]
conf_range = thresholds[f'{band}_20deg'][2]
a = (refang - lo_ang) / (hi_ang - lo_ang)
midpt = lo_ang_val + a*(hi_ang_val - lo_ang_val)
sunglint[1] = midpt
sunglint[2] = midpt - conf_range
sunglint[0] = midpt + conf_range
sunglint[3] = power
return sunglint
def conf_test(rad=_test_rad, thr=_test_thr):
'''
Assuming a linear function between min and max confidence level, the plot below shows
how the confidence (y axis) is computed as function of radiance (x axis).
This case illustrates alpha < gamma, obviously in case alpha > gamma, the plot would be
flipped.
gamma
c 1 ________
o | /
n | /
f | /
i | beta /
d 1/2 |....../
e | /
n | /
c | /
e 0________/
| alpha
--------- radiance ---------->
'''
coeff = np.power(2, (thr[3] - 1))
hicut = thr[0]
beta = thr[1]
locut = thr[2]
power = thr[3]
radshape = rad.shape
rad = rad.reshape((rad.shape[0]*rad.shape[1]))
c = np.zeros(rad.shape)
if hicut > locut:
gamma = thr[0]
alpha = thr[2]
flipped = False
else:
gamma = thr[2]
alpha = thr[0]
flipped = True
# Rad between alpha and beta
range_ = 2. * (beta - alpha)
s1 = (rad[rad <= beta] - alpha) / range_
if flipped is False:
c[rad <= beta] = coeff * np.power(s1, power)
if flipped is True:
c[rad <= beta] = 1. - (coeff * np.power(s1, power))
# Rad between beta and gamma
range_ = 2. * (beta - gamma)
s1 = (rad[rad > beta] - gamma) / range_
if flipped is False:
c[rad > beta] = 1. - (coeff * np.power(s1, power))
if flipped is True:
c[rad > beta] = coeff * np.power(s1, power)
# Rad outside alpha-gamma interval
if flipped is False:
c[rad > gamma] = 1
c[rad < alpha] = 0
if flipped is True:
c[rad > gamma] = 0
c[rad < alpha] = 1
c[c > 1] = 1
c[c < 0] = 0
confidence = c.reshape(radshape)
return confidence
def conf_test_dble(rad, coeffs):
# '''
# gamma1 gamma2
# c 1_______ ________
# o | \ /
# n | \ /
# f | \ /
# i | \ beta1 beta2 /
# d 1/2 \....| |...../
# e | \ /
# n | \ /
# c | \ /
# e 0 \_____________/
# | alpha1 alpha2
# --------------------- radiance ------------------------->
# '''
hicut = [coeffs[0], coeffs[1]]
locut = [coeffs[2], coeffs[3]]
midpt = [coeffs[4], coeffs[5]]
power = coeffs[6]
gamma1 = hicut[0]
gamma2 = hicut[1]
alpha1 = locut[0]
alpha2 = locut[1]
beta1 = midpt[0]
beta2 = midpt[1]
coeff = np.power(2, (power - 1))
radshape = rad.shape
rad = rad.reshape((rad.shape[0]*rad.shape[1]))
c = np.zeros(rad.shape)
# Find if interval between inner cutoffs passes or fails test
if (alpha1 - gamma1 > 0):
# Value is within range of lower set of limits
range_ = 2 * (beta1 - alpha1)
s1 = (rad[(rad <= alpha1) & (rad >= beta1)] - alpha1) / range_
c[(rad <= alpha1) & (rad >= beta1)] = coeff * np.power(s1, power)
range_ = 2 * (beta1 - gamma1)
s1 = (rad[(rad <= alpha1) & (rad < beta1)] - gamma1) / range_
c[(rad <= alpha1) & (rad < beta1)] = coeff * np.power(s1, power)
# Value is within range of upper set of limits
range_ = 2 * (beta2 - alpha2)
s1 = (rad[(rad > alpha1) & (rad <= beta2)] - alpha2) / range_
c[(rad > alpha1) & (rad <= beta2)] = coeff * np.power(s1, power)
range_ = 2 * (beta2 - gamma2)
s1 = (rad[(rad > alpha1) & (rad > beta2)] - gamma2) / range_
c[(rad > alpha1) & (rad > beta2)] = coeff * np.power(s1, power)
# Inner region fails test
# Check for value beyond function range
c[(rad > alpha1) & (rad < alpha2)] = 0
c[(rad < gamma1) | (rad > gamma2)] = 1
else:
# Value is withing range of lower set of limits
range_ = 2 * (beta1 - alpha1)
s1 = (rad[(rad <= gamma1) & (rad <= beta1)] - alpha1) / range_
c[(rad <= gamma1) & (rad <= beta1)] = coeff * np.power(s1, power)
range_ = 2 * (beta1 - gamma1)
s1 = (rad[(rad <= gamma1) & (rad > beta1)] - gamma1) / range_
c[(rad <= gamma1) & (rad > beta1)] = coeff * np.power(s1, power)
# Value is within range of upper set of limits
range_ = 2 * (beta2 - alpha2)
s1 = (rad[(rad > gamma1) & (rad >= beta2)] - alpha2) / range_
c[(rad > gamma1) & (rad >= beta2)] = coeff * np.power(s1, power)
range_ = 2 * (beta2 - gamma2)
s1 = (rad[(rad > gamma1) & (rad < beta2)] - gamma2) / range_
c[(rad > gamma1) & (rad < beta2)] = coeff * np.power(s1, power)
# Inner region passes test
# Check for value beyond function range
c[(rad > gamma1) & (rad < gamma2)] = 1
c[(rad < alpha1) | (rad > alpha2)] = 0
c[c > 1] = 1
c[c < 0] = 0
confidence = c.reshape(radshape)
return confidence
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