The observation operators considered so far assume that the signal from the sensor only reflects the nature of the land surface. In practice, atmospheric particles will absorb and scatter radiation. This processes will contaminate the acquired land surface reflectance. Additionally, there is an interest in understanding the distribution and dynamics of atmospheric components. Understanding the effect of atmospheric processes on radiation gives insight on the distribution of aerosols, gaseous species, water vapour, etc. A number of physically based models have been developed to explain gaseous absorption and scattering. Without going into too much detail, it can be said that these models couple the atmosphere with the land surface, and that in general, there is a requirement to understand the BRDF of the land surface in order to correct for atmospheric effects. So one way of going about accounting for the effect of the atmosphere relies on having an estimate of the surface’s BRDF and then solving for the parameters that result in the measured signal at the sensor. The practicity of this is limited, as we are actually interested in the surface’s BRDF! In this case, simple models of the surface are prescribed (i.e., assume a Lambertian surface with a typical reflectance value for vegetation, bare soil or water). This “Lambertian assumption” only requires then some atmospheric composition parameters to account for the atmosphere.
A more complicated setup couples the atmosphere and the land surface, resulting in an extended state vector, where the parameters that define the atmospheric composition are now part of the state vector. This of course requires an strategy for the coupling.
In the frame of EOLDAS, we use the 6S atmospheric model. This model has been widely utilised to carry out atmospheric corrections. In order to simplify things, some decisions have been made to limit the complexity of the model:
The computational cost of the model makes the use of monochromatic wavebands necessary, although this can be modified. Using a Lambertian surface requires that the concentration of water vapour, ozone and the aerosol optical depth at 550nm are specified. The model then is able to predict what the “at sensor” radiances will be. The BRDF-coupled mode requires samples of the surface BRDF and hemispherical albedo. These are then fed into 6S, in addition to the atmospheric parameters detailed above, to provide a calculation of the “at sensor” radiance.
The coupling of 6S with the land surface is made through deriving the Observation_Operator class. Under the Lambertian assumption and assuming a monochromatic sensor, the only parameters that control the atmosphere as modelled by 6S are water content, ozone concentration and aerosol optical thickness at 550nm. These are specified in a default configuration file, and the state vector is enlarged by these three parameters for each date. This is shown in the configuration file semid_atcorr_default.conf:
# configuration file
# note that some words are protected: values, keys
# so dont use them
# other than that, you can set any level of hierarchy in the
# representation
# this next part serves no purpose other than
# to show how different options can be set
# It is a requirement that parameter.names be defined
[parameter]
location = ['time']
limits = [[1,365,1]]
names=gamma_time,uwc,uo3,aot,xlai, xhc, rpl, xkab, scen, xkw, xkm, xleafn, xs1,xs2,xs3,xs4,lad
solve = [1]*len($parameter.names)
help_solve='flags for which state vector elements to solve for'
datatypes=x
[parameter.result]
filename = 'test/semid/output/test.params'
help_filename="state vector results file"
format = 'PARAMETERS'
[parameter.assoc_solve]
gamma_time = 0
uo3 = 0
lad = 0
xs3 = 0
xs4 = 0
xhc = 0
rpl = 0
[parameter.x]
datatype = x
names = $parameter.names
default = [0.05]*len($parameter.names)
help_default = "Set the parameter default values"
apply_grid = True
sd = [1.]*len($parameter.names)
bounds = [[0.01,0.99]]*len($parameter.names)
####state = fullState.dat
[parameter.x.assoc_bounds]
gamma_time = 0.000001,100000
uwc = 0.05, 5.
uo3 = 0.2, 0.45
aot = 0.1, 0.99
xlai = 0.05,0.99
xhc = 0.05,10.0
rpl = 0.01,0.10
xkab = 0.01,0.99
scen = 0.001,1.0
xkw = 0.01,0.99
xkm = 0.3,0.9
xleafn = 0.9,2.5
xs1 = 0.01, 4.
xs2 = 0.01, 5.
xs3 = 0., 0.
xs4 = 0.,0.
lad = 1,5
[parameter.x.assoc_default]
####gamma_time = 0.01
uwc = 1.
uo3 = 0.390
aot = 0.25
xlai = 0.95
xhc = 5
rpl = 0.01
xkab = 0.95
scen = 0.001
xkw = 0.95
xkm = 0.35
xleafn = 1.5
xs1 = 1.0
xs2 = 0.001
xs3 = 0
xs4 = 0
lad = 5
help_lad = 'lad value'
[general]
datadir = test/semid/input/,.
is_spectral = True
grid=True
calc_posterior_unc=False
write_results=True
doplot=False
plotmod=0
plotmovie=False
[general.optimisation]
randomise=False
[operator]
modelt.name=DModel_Operator
modelt.datatypes = x
obs.name=Observation_Operator_AtCorr
obs.datatypes = x,y
[operator.obs.rt_model]
model=rtmodel_ad_trans2
use_median=True
help_use_median = "Flag to state whether full bandpass function should be used or not.\nIf True, then the median wavelength of the bandpass function is used"
bounds = [400,2500,1]
help_bounds = "The spectral bounds (min,max,step) for the operator'
ignore_derivative=False
help_ignore_derivative = "Set to True to override loading any defined derivative functions in the library and use numerical approximations instead"
[operator.modelt.x]
names = $parameter.names
sd = [1.0]*len($operator.modelt.x.names)
datatype = x
[operator.modelt.rt_model]
model_order=2
help_model_order='The differential model order'
wraparound=periodic,365
help_wraparound='Wraparound conditions (periodic,reflexive,none). If periodic, then the period can be specified as a second a
rgument (e.g. --operator.modelt.rt_model.wraparound=periodic,365)'
lag=1
help_lag='The grid lag spacing between samples for the differential operator. This can be a list, in which case lag weight is
also required'
lag_weight = 1
help_lag_weight='The weight for the lag. Only of value if multiple lags are specified'
[operator.obs.x]
names = $parameter.names[1:]
sd = [1.0]*len($operator.obs.x.names)
datatype = x
[operator.obs.y]
control = 'mask vza vaa sza saa'.split()
names = "465.6 553.6 645.5 856.5 1241.6 1629.1 2114.1".split()
help_names='Band names'
sd = "0.003 0.004 0.004 0.015 0.013 0.01 0.006".split()
help_sd='Observations StdDev'
datatype = y
state = test/data_type/input/test.brf
help_state='set the obs state file'
# angle_transform not yet implemented
#angle_transform = False
#help_angle_transform = "Transform the relative azimuth of the input data by adding 180 degrees. Useful if model and data are not consistently defined."
[operator.obs.y.result]
filename = 'test/semid/output/test-fwd.brf'
help_filename = 'forward modelling results file'
format = 'PARAMETERS'
The previous file does three main things:
When run, it will read parameters from filename = ‘test/semid/output/test.params’, and it will forward model the top of atmosphere radiance for the sensor specified in operator.obs.y (and its associated data file, state= test/data_type/input/test.brf.
Although the Lambertian assumption only requires a single estimate of the surface reflectance, setting up the atmosphere and performing the correction can be quite time consuming. The lack of an adjoint to the coupled atmosphere-surface continuum requires that if the atmospheric parameters need to be solved for, a numeric approximation to the observation operator’s cost function will need to be calculated, making the proceedings very slow. However, once the land surface estate is estimated, forward running the surface parameters to at-sensor radiances is possible and practical
Another option to do an atmospheric correction is to couple the surface’s BRDF with the atmosphere. In 6S, this is done by providing adequate sampling of the BRDF and spherical albedo, which in our case is realised by having multiple forward model runs of the land surface observation operator. The reflectances resulting from these multiple runs are fed into 6S, in addition to atmospheric composition parameters, and the at-sensor radiance is returned. The operation is analogous to the Lambertian case, and only the name of the Observation Operator is different: we set operator.obs.nameObservation_Operator_AtCorr_Cpled. A configuration file, semid_atcorr_cpld_default.conf is provided for this. This is essentially the same as semid_atcorr_default.conf, it only defines a different observation operator:
# configuration fle
# note that some words are protected: values, keys
# so dont use them
# other than that, you can set any level of hierarchy in the
# representation
# this next part serves no purpose other than
# to show how different options can be set
# It is a requirement that parameter.names be defined
[parameter]
location = ['time']
limits = [[1,365,1]]
names=gamma_time,uwc,uo3,aot,xlai, xhc, rpl, xkab, scen, xkw, xkm, xleafn, xs1,xs2,xs3,xs4,lad
solve = [1]*len($parameter.names)
help_solve='flags for which state vector elements to solve for'
[parameter.result]
filename = 'test/semid/output/test.params'
help_filename="state vector results file"
format = 'PARAMETERS'
[parameter.assoc_solve]
gamma_time = 0
lad = 0
xs3 = 0
xs4 = 0
xhc = 0
rpl = 0
[parameter.x]
datatype = x
names = $parameter.names
default = [0.05]*len($parameter.names)
help_default = "Set the parameter default values"
apply_grid = True
sd = [1.]*len($parameter.names)
bounds = [[0.01,0.99]]*len($parameter.names)
####state = fullState.dat
[parameter.x.assoc_bounds]
gamma_time = 0.000001,100000
uwc = 0.001, 0.9
uo3 = 0.001, 0.9
aot = 0.1, 0.99
xlai = 0.05,0.99
xhc = 0.05,10.0
rpl = 0.01,0.10
xkab = 0.01,0.99
scen = 0.001,1.0
xkw = 0.01,0.99
xkm = 0.3,0.9
xleafn = 0.9,2.5
xs1 = 0.01, 4.
xs2 = 0.01, 5.
xs3 = 0., 0.
xs4 = 0.,0.
lad = 1,5
[parameter.x.assoc_default]
####gamma_time = 0.01
uwc = 0.1
uo3 = 0.1
aot = 0.2
xlai = 0.95
xhc = 5
rpl = 0.01
xkab = 0.95
scen = 0.001
xkw = 0.95
xkm = 0.35
xleafn = 1.5
xs1 = 1.0
xs2 = 0.001
xs3 = 0
xs4 = 0
lad = 5
help_lad = 'lad value'
[general]
datadir = test/semid/input/,.
is_spectral = True
calc_posterior_unc=False
write_results=True
doplot=True
plotmod=10
plotmovie=True
[general.optimisation]
randomise=False
[operator]
modelt.name=DModel_Operator
modelt.datatypes = x
obs.name=Observation_Operator_AtCorrCpld
obs.datatypes = x,y
[operator.obs.rt_model]
model=rtmodel_ad_trans2
use_median=True
help_use_median = "Flag to state whether full bandpass function should be used or not.\nIf True, then the median wavelength of the bandpass function is used"
bounds = [400,2500,1]
help_bounds = "The spectral bounds (min,max,step) for the operator'
ignore_derivative=True
help_ignore_derivative = "Set to True to override loading any defined derivative functions in the library and use numerical approximations instead"
[operator.modelt.x]
names = $parameter.names
sd = [1.0]*len($operator.modelt.x.names)
datatype = x
[operator.modelt.rt_model]
model_order=2
help_model_order='The differential model order'
wraparound=periodic,365
help_wraparound='Wraparound conditions (periodic,reflexive,none). If periodic, then the period can be specified as a second a
rgument (e.g. --operator.modelt.rt_model.wraparound=periodic,365)'
lag=1
help_lag='The grid lag spacing between samples for the differential operator. This can be a list, in which case lag weight is
also required'
lag_weight = 1
help_lag_weight='The weight for the lag. Only of value if multiple lags are specified'
[operator.obs.x]
names = $parameter.names[1:]
sd = [1.0]*len($operator.obs.x.names)
datatype = x
[operator.obs.y]
control = 'mask vza vaa sza saa'.split()
names = "465.6 553.6 645.5 856.5 1241.6 1629.1 2114.1".split()
sd = "0.003 0.004 0.004 0.015 0.013 0.01 0.006".split()
datatype = y
state = test/data_type/input/test.brf
help_state='set the obs state file'
# angle_transform not yet implemented
#angle_transform = False
#help_angle_transform = "Transform the relative azimuth of the input data by adding 180 degrees. Useful if model and data are not consistently defined."
[operator.obs.y.result]
filename = 'test/semid/output/test-fwd.brf'
help_filename = 'forward modelling results file'
format = 'PARAMETERS'
Unfortunately, while providing a more detailed description of the surface will improve the atmospheric correction, the lack of an adjoint coupled with the large number of forward land surface model runs that are required to sample the surface BRDF, make this methodology computationally very expensive.