«  A tutorial guide of the EOLDAS prototype software   ::   Contents   ::   Radiative transfer modelling for optical remote sensing  »

A simple observation operator example, running in eoldas

Purpose of this section

This section of the user guide will take a simple example of combining two cost functions in a variational DA sense: an Identity (or Prior) operator and a smoothness (derivative) constraint. The examples used here are of filtering noisy time series of NDVI data.

As well as learning about these concepts if you are unfamiliar with them, this section will also take you through some examples of running the eoldas to solve this sort of problem. Towards the end of the section, we delve into writing some python code to make use of eoldas, and also mention other ways of interacting with it.

We learn that at the heart of eoldas is one or more configuration files that allow us to solve problems without the need for any code writing (though you can if you want to!).

Introducing the observation operator concept

In the previous sections, we have assumed that the state of the land surface can be observed directly, and that these observations are only limited by noise in their acquisition. In general, the state of the surface and the observations will be different entities. For example, the state of the surface may include parameters such as leaf area index, chlorophyll concentration or soil moisture, whereas the observations may be of directional surface reflectance, top of atmosphere radiance or microwave temperatures. The link between these magnitudes is the observation operator, which maps the land surface parameters into quantities that are measured by a sensor. Many types of observation operator are available: from the statistical to the physics based. The simplest case, however, is the identity operator. In fact, we have already used it, as we have assumed that the observations are just direct measurements of the state vector.

A simple data assimilation example using an identity observation operator

The simplest observation operator is the identity observation operator, where the observations are identical to the state vector components. We can see the use of this operator as a way of optimally smoothing univariate timeseries, for example NDVI. Additionally, the use of a DA system allows to interpolate where data points are not available. The following demonstrates how the EOLDAS prototype can be used for this task, and also allows us to explore the use of the weak assimilation paradigm conveniently.

The main way to run EOLDAS is via one or more configuration files. This is partly to make sure that a record exists of a particular experimental setup and partly to allow flexible running of the system without the need for further coiding by the user (unless he/she wants to add new classes or methods).

Here is the start of the configuration file that we are going to use:

# An EOLDAS confioguration file for
# a simple Identity operator and a regulariser

[parameter]
location = ['time']
limits = [[1,365,1]]
names = gamma_time,NDVI
solve = 0,1
datatypes = x

[parameter.result]
filename = output/Identity/NDVI_Identity.params
format = 'PARAMETERS'
help_filename='Set the output state filename'

[parameter.x]
datatype = x
names = $parameter.names
default = 200,0
help_default="Set the default values of the states"
apply_grid = True
sd = [1.]*len($parameter.names)
bounds =  [[0.000001,1000000],[-1,1]]

This sets up the section [parameter], which described the state variables within EOLDAS.

The section must contain fields describing the locational information of the data to be used as well as the names of the state variables we will consider. In this case, we have two state variables gamma_time and NDVI that we wish to estimate over the (time) range 1 to 365 (inclusive) in steps of 1 (day).

The subsection parameter.solve gives a list of codes indicating whether we wish to solve for each parameter or not.

A code of 0 tells the EOLDAS not to solve, i.e. just to use the data that are read in or any other default values.

A code of 1 tells EOLDAS to solve for that state variable at all locations (times here).

A code of 2 tells EOLDAS to solve for the state variable, but to maintain a single value over all locations.

The section parameter.result gives information on any output file name and format for the state.

Finally in this section, the field parameter.x sets up data and conditions for the state vector x. Remember that it is this state vector x that we will solve for in the EOLDAS. Here, we specify that it is of datatype x (i.e. the x state vector), that it has the same names as we set up in parameter.names, that default values to be assigned are 25 and 0 for the two variables respectively. If any data are read in, these override the default values, but in this case, we simply start with the defaults.

The text help_default allows the field parameter.x.default to be set from the command line with the option --parameter.x.default=a,b.

The flag apply_grid, which is the default, tells the EOLDAS to produce the state vector on a grid over the bounds defined in parameter.limits.

Finally, we define default uncertainty information for the state vector (this does not directly affect the running of the EOLDAS (parameter.x.sd), and define the bounds for each state vector (use None if no bound is to be used). Here, we set the lower bound as 0 and the upper bound as 1 for both parameters.

The next section sets up general conditions:

[general]
doplot=True
help_do_plot='do plotting'

In this case, we set a flag to do plotting when the results are written out. Plotting will use the filenames in any state variable section e.g. parameter.result.filename to generate a series of plots with filenames the same as this data filename. We will see some examples later.

The next section sets up the operators that we want to define here.

[operator]
modelt.name=DModel_Operator
modelt.datatypes = x
obs.name=Operator
obs.datatypes = x,y

Here, we define two operators, DModel_Operator and Operator. These names refer directly to python classes for the operators in EOLDAS. The base class is ‘Operator’ which implements the Identity operator. All other classes are derived from this. The differential operator works only on the x state vector, which is equlivalent to defining \(\mathbf{y_{obs}}=0\). The operator ‘Operator’ access both x and y data if it to act as a Prior constraint, so we set up x and y datatypes.

Next we set the details of these operators. First, the differential operator:

[operator.modelt.x]
names = $parameter.names
sd = [1.0]*len($operator.modelt.x.names)
datatype = x

[operator.modelt.rt_model]
model_order=1
help_model_order='The differential model order'
wraparound=periodic,365

where we specify which state vector elements this operator has access to (all of those in parameter.names here) and set up the default uncertainty and datatype.

We then set the parameters specific to the ‘model’ \(H(\mathbf{x})\), in this case the order of the differential model (2 here) and the edge conditions (periodic, with a period of 200 (days)).

Finally, we set up the operator operator.obs, specifically, parameters for its x and y state vectors.

[operator.obs.x]
names = $parameter.names[1:]
sd = [1.0]*len($operator.obs.x.names)
help_sd='Set the observation sd'
datatype = x

[operator.obs.y]
control = [mask]
names = $parameter.names[1:]
sd = [0.15]*len($operator.obs.x.names)
help_sd="set the sd for the observations"
datatype = y
state = data/Identity/random_ndvi1.dat
help_state = "Set the y state vector"

[operator.obs.y.result]
filename = output/Identity/NDVI_fwd.params

specifying default uncertainty information, data types and any required output files. Here, we wish to write out the results in operator.obs.y, so we specify a filename for this. Again, we see the use of a ‘help’ variable, which here allows operator.obs.y.result.filename to be set from the command line. This interfacing to the command line means that a single configuration file can generally serve for multiple experiments and the user does not need to keep generating new ones.

The main program can be accessed in various ways. One way is to write some front end code that calls the eoldas python code.

An example is solve_eoldas_identity.py that includes three sections:

First, some code to generate a synthetic dataset.

import pdb
import numpy as np

def create_data ( n_per=4, noise=0.15, obs_off=0.33, \
					window_size=0, order=4):
    """
    Create synthetic "NDVI-like" data for a fictitious time series. We return
    the original data, noisy data (using IID Gaussian noise), the QA flag as well
    as the time axis.
    
    Missing observations are simulated by drawing a random number between 0 and 1
    and checking against obs_off.
    
    Parameters
    ----------
    n_per : integer
    Observation periodicity. By default, assumes every 8 days

    noise : float
    The noise standard deviation. By default, 0.15

    obs_off : float
    The threshold to decide on missing observations in the time series.

    window_size : integer, odd
    window size for savitzky_golay filtering. A large window size will lead
    to larger data gaps by correlating the noise. Set to zero by default
    which applies no smoothing.

    order : integer
    order of the savitzky_golay filter. By default 4.
    
    """
    from savitzky_golay import savitzky_golay
    import numpy as np

    doys  = np.arange ( 1, 365+1, n_per)
    ndvi_clean =  np.clip(np.sin((doys-1)/72.), 0,1) 
    ndvi =  np.clip(np.sin(doys/72.), 0,1) 
    # add Gaussian noise of sd noise
    ndvi = np.random.normal(ndvi,noise,ndvi.shape[0])
     
    # set the qa flags for each sample to 1 (good data)
    qa_flag = np.ones_like ( ndvi).astype( np.int32 )
    passer = np.random.rand( ndvi.shape[0])
    if window_size >0:
        # force odd
        window_size = 2*(window_size/2)+1
        passer = savitzky_golay(passer, window_size=window_size, order=order)
    # assign a proportion of the qa to 0 from an ordering of the smoothed 
    # random numbers
    qa_flag[np.argsort(passer)[:passer.size * obs_off]]  = 0

    return ( doys, ndvi_clean, ndvi, qa_flag )

Here, we generate some NDVI data which has a trajectory of a sine wave for the first half of the year and is flat at zero for the second half. The sampling is controlled by n_per and obs_off. The parameter obs_off randomly removes a proportion of the data. If window_size and order are set then a savitzky_golay filter is used to induce correlation in the timing of the samples that are removed from thie dataset (qa=0). This mimics what we practically have in Optical Earth Observation with temporal correlation in cloud cover.

The next section calculates the ‘ideal’ value of \(\gamma\) by calculating the root mean squared deviation of the original dataset. We use this ideal gamma here for to demonstrate the physical meaning of the gamma value, though in practice this would be unknown. This section also writes the dataset to a temporary file in ‘BRDF’ format.

All of this so far standard python coding, though such datasets could be generated in many other ways. The final section interfaces to the top level of the eoldas code, which is what is of immediate concern in this tutorial.

    import sys,tempfile
    this = sys.argv[0]
    import  eoldas

    # SD of noise
    noise=0.15
    # nominal sampling period
    n_per=7
    # proportion of missing samples
    obs_off=0.33
    # order of differential model (integer)
    model_order=1
    # sgwindow is larger to create larger data gaps
    sgwindow=10

    # set up data for this experiment
    file, ideal_gamma,doys,ndvi_clean,ndvi,qa_flag  = \
		prepare(noise=noise,n_per=n_per,\
			obs_off=obs_off,model_order=model_order,\
			sgwindow=sgwindow)

    # set gamma to less than the the theoretical value
    gamma = ideal_gamma*0.33

    # set op file names
    xfile = 'output/Identity/NDVI_Identity.params'
    yfile = 'output/Identity/NDVI_fwd.params'

    # initialise optuions for DA overriding any in config files
    cmd = 'eoldas ' +  \
                ' --conf=config_files/eoldas_config.conf --conf=config_files/Identity.conf ' + \
                ' --logfile=mylogs/Identity.log ' + \
                ' --calc_posterior_unc ' + \
                ' --parameter.solve=0,1 ' + \
		' --parameter.result.filename=%s '%xfile +\
                ' --parameter.x.default=%f,0.0 '%(gamma) + \
		' --operator.obs.y.result.filename=%s'%yfile +\
                ' --operator.obs.y.state=%s'%file+\
                ' --operator.modelt.rt_model.model_order=%d '%model_order 

    # initialise eoldas
    self = eoldas.eoldas(cmd)
    # solve DA
    self.solve(write=True)
    

The first part of the code extends the system path for where it searches for libraries. This is done relative to where solve_eoldas_identity.py is (in the bin directory of the distribution). After that, we set up values for the parameters for generating the synthetic dataset.

The interface to the eoldas here is mainly to make a string with a number of flags. The most important flag is --conf=config_files/Identity.conf which specifies the configuration file for this experiment. In addition, --conf=eoldas_config.conf is given, which specifies a system default configuration file, eoldas_config.conf.

So, the simplest ‘top level’ interface to eoldas from python code involves:

    gamma = ideal_gamma*0.33

    # set op file names
    xfile = 'output/Identity/NDVI_Identity.params'
    yfile = 'output/Identity/NDVI_fwd.params'

    # initialise optuions for DA overriding any in config files
    cmd = 'eoldas ' +  \
                ' --conf=config_files/eoldas_config.conf --conf=config_files/Identity.conf ' + \
                ' --logfile=mylogs/Identity.log ' + \
                ' --calc_posterior_unc ' + \
                ' --parameter.solve=0,1 ' + \
		' --parameter.result.filename=%s '%xfile +\
                ' --parameter.x.default=%f,0.0 '%(gamma) + \
		' --operator.obs.y.result.filename=%s'%yfile +\
                ' --operator.obs.y.state=%s'%file+\
                ' --operator.modelt.rt_model.model_order=%d '%model_order 

and

                ' --operator.modelt.rt_model.model_order=%d '%model_order 

    # initialise eoldas
    self = eoldas.eoldas(cmd)
    # solve DA
    self.solve(write=True)

where we set up the text string, initiate the eoldas object (eoldas.eoldas) and then call the eoldas.solve() method.

This command string is clearly of some importance. The flags --logfile=mylogs/Identity.log and --calc_posterior_unc correspond to items in the [general] section of the configuration file. Here, eoldas_config.conf contains the lines:

[general]
datadir = ., data
help_datadir ="Specify where the data and or conf files are"
here = os.getcwdu()
grid = True
is_spectral = True
calc_posterior_unc=False
help_calc_posterior_unc ="Switch to calculate the posterior uncertainty"
write_results=True
help_write_results="Flag to make eoldas write its results to files"
init_test=False
help_init_test="Flag to make eoldas run a test of the cost functions before proceeding with DA"
doplot=False
plotmod=100000000
plotmovie=False

[general.optimisation]
# These are the default values
iprint=1
gtol=1e-3
help_gtol = 'set the relative tolerance for termination'
maxIter=1e4
maxFunEvals=2e4
plot=0
# see http://openopt.org/NLP#Box-bound_constrained
solverfn=scipy_lbfgsb 
randomise=False
help_randomise='Randomise the starting point'

To understand how e.g the flag --calc_posterior_unc can be used we can look at:

calc_posterior_unc=False
help_calc_posterior_unc ="Switch to calculate the posterior uncertainty"

where we set the default value of the item (False) and also have a ‘help’ statement which allows this value to be overridden. You shouldn’t normally need to change things in the system configuration file eoldas_config.conf. This flag controls whether we calculate the posterior iuncertainty or not. The default is False because it can be quite computationally expensive and is often best done in a post processing step.

The flag --operator.obs.y.state=filename refers specifically to the section operator.obs.y.state of the configuration, which is something we have set up in Identity.conf with a ‘help’ field, so we can override the default value set in the configuration file.

[operator.obs.y]
control = [mask]
names = $parameter.names[1:]
sd = [0.15]*len($operator.obs.x.names)
help_sd="set the sd for the observations"
datatype = y

If we run solve_eoldas_identity.py, it sends logging information to mylogs/Identity.log and reports (to the stderr) the progress of the optimisation, ‘f’ being the total of all of the \(J\) terms for this configuration. It should converge to a solution within some tens of iterations and result in a final value of \(J\) of around 1800. We set the name of the logfile in solve_eoldas_identity.py:

    if window_size >0:

There is a lot of detail in the log file about exactly what value terms are set to and the progress of the eoldas. It also contains information on the individual \(J\) terms:

2012-06-14 15:15:57,455 - eoldas.solver.eoldas.solver-obs-x - INFO -      J   = 293.213607
2012-06-14 15:15:57,456 - eoldas.solver.eoldas.solver-modelt-x - INFO -      J   = 34.354830
2012-06-14 15:15:57,457 - eoldas.solver.eoldas.solver-obs-x - INFO -      J   = 196.657803
2012-06-14 15:15:57,458 - eoldas.solver.eoldas.solver-modelt-x - INFO -      J   = 87.225910
2012-06-14 15:15:57,459 - eoldas.solver.eoldas.solver-obs-x - INFO -      J   = 112.026943
2012-06-14 15:15:57,461 - eoldas.solver.eoldas.solver-modelt-x - INFO -      J   = 113.479088
2012-06-14 15:15:57,462 - eoldas.solver.eoldas.solver-obs-x - INFO -      J   = 78.729997
2012-06-14 15:15:57,463 - eoldas.solver.eoldas.solver-modelt-x - INFO -      J   = 93.761052
2012-06-14 15:15:57,464 - eoldas.solver.eoldas.solver-obs-x - INFO -      J   = 43.172989
2012-06-14 15:15:57,466 - eoldas.solver.eoldas.solver-modelt-x - INFO -      J   = 99.078418
2012-06-14 15:15:57,467 - eoldas.solver.eoldas.solver-obs-x - INFO -      J   = 40.564246
2012-06-14 15:15:57,468 - eoldas.solver.eoldas.solver-modelt-x - INFO -      J   = 74.685441
2012-06-14 15:15:57,469 - eoldas.solver.eoldas.solver-obs-x - INFO -      J   = 39.829009
2012-06-14 15:15:57,470 - eoldas.solver.eoldas.solver-modelt-x - INFO -      J   = 63.881006
2012-06-14 15:15:57,471 - eoldas.solver.eoldas.solver-obs-x - INFO -      J   = 27.079490

A log file is important to the running of eoldas, as the processing can take quite some time for some problems.

The state vector results will be written to output/Identity/NDVI_Identity.params because we first specified this in Identity.conf:

[parameter.result]
filename = output/Identity/NDVI_Identity.params
format = 'PARAMETERS'

(actually, since we also set the comamnd ' --operator.obs.y.result.filename=%s'%yfile in the code we have written, this will override what is in the parameter file).

The output file format is a general ‘PARAMETERS’ format that should look something like this:

The first line of the file is a header statement that describes the location fields (‘time’ here) and state variables (gamma_time and NDVI here) and subsequent lines give the values for those fields. This format can be used for input data, but it cannot at present be used to define a full input covariance matrix, only the standard deviation for each observation as in this example.

A graph of the result is in output/Identity/NDVI_Identity.params.plot.x.png:

output/Identity/NDVI_Identity.params.plot.x.png

This result is quite interesting for understanding how our DA system works: The plot shows the mean of the estimate of the NDVI state vector in the lower panel (the upper panel shows the value of gamma_time which we did not solve for) as a red line as we solve for the state every day (of 365 days). The 95% confidence interval is shown shaded in grey. We will show this below in a more refined plot of the results.

We also specified the ‘y’ data to be written out (in Identity.conf):

help_state = "Set the y state vector"

[operator.obs.y.result]
filename = output/Identity/NDVI_fwd.params

so this goes to a file output/Identity/NDVI_fwd.params unless the flag --operator.obs.y.result.filename=somethingElse.dat is used. The format of this file is the same as above, but it shows the retrieved NDVI data since this reports \(H(x)\).

#PARAMETERS time mask NDVI sd-NDVI
   1.000000    1.000000    0.211816    0.000000
   8.000000    1.000000    0.275085    0.000000
  15.000000    1.000000    0.336462    0.000000
  22.000000    1.000000    0.401630    0.000000
  29.000000    1.000000    0.455163    0.000000
  36.000000    1.000000    0.542784    0.000000
  43.000000    1.000000    0.610468    0.000000
  50.000000    1.000000    0.669554    0.000000
  92.000000    1.000000    0.834135    0.000000

The original dataset is output form convenience in the same format as the y state, being output/Identity/NDVI_fwd.params_orig here:

#PARAMETERS time mask NDVI sd-NDVI
   1.000000    1.000000    0.315700    0.150000
   8.000000    1.000000    0.286300    0.150000
  15.000000    1.000000    0.313996    0.150000
  22.000000    1.000000    0.470529    0.150000
  29.000000    1.000000    0.253425    0.150000
  36.000000    1.000000    0.660812    0.150000
  43.000000    1.000000    0.661460    0.150000
  50.000000    1.000000    0.856904    0.150000
  92.000000    1.000000    0.894455    0.150000

Various other graphics are output for a ‘y’ state:

A plot at each location, in output/Identity/NDVI_fwd.params.plot.y.png (not of much relevance in this example):

output/Identity/NDVI_fwd.params.plot.y.png

A plot of the observations and modelled values as a function of location (the observations are the green dots) in output/Identity/NDVI_fwd.params.plot.y2.png:

output/Identity/NDVI_fwd.params.plot.y2.png

and a plot of the x state vector associated with this operator (NDVI here).

output/Identity/NDVI_fwd.params.plot.x.png

Example plotting data from the output files

The above plots are automatically generated by eoldas provided general.doplot is True but these are intended as quicklooks, and users are likely to want to form their own plots.

An example of this is implemented in example1plot.py:

#!/usr/bin/env python
import numpy as np
import pylab as plt
#
# Some plotting of the synthetic and retrieved data

# generate the clean data
doys  = np.arange ( 1, 365+1, 1)
ndvi_clean =  np.clip(np.sin((doys-1)/72.), 0,1)

# read the state file
state = open('output/Identity/NDVI_Identity.params').readlines()
sdata = np.array([np.array(i.split()).astype(float) for i in state[1:]])
sNdvi = sdata[:,2]
sdNdvi = sdata[:,4]

# noisy sample data
noisy = open('output/Identity/NDVI_fwd.params_orig').readlines()
ndata = np.array([np.array(i.split()).astype(float) for i in noisy[1:]])
nNdvi = ndata[:,2]
sdnNdvi = ndata[:,3]
ndoys = ndata[:,0]

# retrieved
plt.fill_between(doys,y1=sNdvi-1.96*sdNdvi,y2=sNdvi+1.96*sdNdvi,facecolor='0.8')
p1 = plt.plot(doys,sNdvi,'r')
# original
p2 = plt.plot(doys,ndvi_clean,'g')
# samples used
p3 = plt.errorbar(ndoys,nNdvi,yerr=sdnNdvi*1.96,fmt='bo')

plt.legend([p2,p3,p1],['Original state','Sampled noisy state','Retrieved state'])

#plt.show()
plt.savefig('images/example1plot.png')
example 1 plot

We can see that smoothers of this sort have some difficulty maintaining sudden changes, although this is quite challenging here given the level of noise amd the rather large data gaps.

More importantly, we have used the smoother to interpolate over the missing observations and to reduce uncertainty.

If we inspect the file output/Identity/NDVI_Identity.params <output/Identity/NDVI_Identity.params>:

#PARAMETERS time gamma_time NDVI sd-gamma_time sd-NDVI
   1.000000   42.912561    0.211816    0.000000    0.063006
   2.000000   42.912561    0.220854    0.000000    0.064164
   3.000000   42.912561    0.229893    0.000000    0.065042
   4.000000   42.912561    0.238931    0.000000    0.065651
   5.000000   42.912561    0.247970    0.000000    0.066000
   6.000000   42.912561    0.257008    0.000000    0.066092
   7.000000   42.912561    0.266046    0.000000    0.065928
   8.000000   42.912561    0.275085    0.000000    0.065507
   9.000000   42.912561    0.283853    0.000000    0.066408

alongside the original data output/Identity/NDVI_fwd.params_orig <output/Identity/NDVI_Identity.params>:

#PARAMETERS time mask NDVI sd-NDVI
   1.000000    1.000000    0.315700    0.150000
   8.000000    1.000000    0.286300    0.150000
  15.000000    1.000000    0.313996    0.150000
  22.000000    1.000000    0.470529    0.150000
  29.000000    1.000000    0.253425    0.150000
  36.000000    1.000000    0.660812    0.150000
  43.000000    1.000000    0.661460    0.150000
  50.000000    1.000000    0.856904    0.150000
  92.000000    1.000000    0.894455    0.150000

We note that the original state here (green line) lies entirely within the 95% CI. The error reduction has been of the order of 2.2 (compare sd-NDVI after the DA with that prior to it). We can see that the uncertainty at the datapoints has been reduced from 0.15 (that of the input data) to typically around 0.065. This grows slightly (to around 0.10) when the data gaps are large.

The impact of ‘filtering’ in this way (optimising a fit of the model to the observations) is to smooth the data and reduce uncertainty in the output. The reduction in uncertainty is related to the amount of smoothing that we apply. The second effect is to interpolate between observations, with uncertainty growing where we have no observations.

You might try changing the value of gamma used and seeing the effect on the results.

We could do this by modifying a few lines of solve_eoldas_identity.py to produce solve_eoldas_identity1.py:

if __name__ == "__main__":
    import sys,tempfile
    this = sys.argv[0]
    import eoldas

    # SD of noise
    noise=0.15
    # nominal sampling period
    n_per=7
    # proportion of missing samples
    obs_off=0.33
    # order of differential model (integer)
    model_order=1
    # sgwindow is larger to create larger data gaps
    sgwindow=10

    # set up data for this experiment
    file, ideal_gamma,doys,ndvi_clean,ndvi,qa_flag  = \
		prepare(noise=noise,n_per=n_per,\
			obs_off=obs_off,model_order=model_order,\
			sgwindow=sgwindow)

    # set gamma to less than the the theoretical value
    gamma = ideal_gamma*0.45

    # set op file names
    xfile = 'output/Identity/NDVI_Identity1.params'
    yfile = 'output/Identity/NDVI_fwd1.params'

    # initialise optuions for DA overriding any in config files

which writes out to output/Identity/NDVI_Identity1.params and output/Identity/NDVI_fwd1.params and has a gamma value that is 0.45/0.33 of that previously used.

Now, plotting this using example1plot1.py:

example 1 plot

This is possibly a better result, but in fact what we see is further limitation of the model that we have chosen here: we enforce wraparound (i.e. the NDVI at day 1 is expected to be the same as at day 365) and we enforce smoothness (so as we increase the gamma, the smoothness, we over-smooth at the sudden change that occurs half way through the year).

We could remove the wraparound condition, but in practice, it is better simply to weaken this constraint. We have done this in solve_eoldas_identity2.py:

if __name__ == "__main__":
    import sys,tempfile
    this = sys.argv[0]
    import eoldas

    # SD of noise
    noise=0.15
    # nominal sampling period
    n_per=7
    # proportion of missing samples
    obs_off=0.33
    # order of differential model (integer)
    model_order=1
    # sgwindow is larger to create larger data gaps
    sgwindow=10

    # set up data for this experiment
    file, ideal_gamma,doys,ndvi_clean,ndvi,qa_flag  = \
		prepare(noise=noise,n_per=n_per,\
			obs_off=obs_off,model_order=model_order,\
			sgwindow=sgwindow)

    # set gamma to less than the the theoretical value
    gamma = ideal_gamma*0.33

    # set op file names
    xfile = 'output/Identity/NDVI_Identity2.params'
    yfile = 'output/Identity/NDVI_fwd2.params'

    # initialise optuions for DA overriding any in config files
    cmd = 'eoldas ' +  \
                ' --conf=config_files/eoldas_config.conf --conf=config_files/Identity.conf ' + \
                ' --logfile=mylogs/Identity.log ' + \
                ' --calc_posterior_unc ' + \
                ' --parameter.solve=0,1 ' + \
		' --parameter.result.filename=%s '%xfile +\
                ' --parameter.x.default=%f,0.0 '%(gamma) + \

which write out to output/Identity/NDVI_Identity2.params and output/Identity/NDVI_fwd2.params and has the same (higher) gamma value used above.

Now, plotting this using example1plot2.py:

example 1 plot

If you re-run these scripts several times, so that you see different configurations for the temporal sampling, you will notice that the interpolation sometimes behaves well over the entire dataset (for some given gamma) and sometimes doesn’t. For example:

example 1 plot

Interfacing a little more deeply with the eoldas code

Whilst considering writing wrapper codes around eoldas functionality and outputs it is instructive to explore some of the data structure available.

We can re-use the example solve_eoldas_identity.py developed above and access some of the data structure as shown in solve_eoldas_identity_a.py.

#!/usr/bin/env python
import pdb
import numpy as np

if __name__ == "__main__":
    import sys,tempfile
    this = sys.argv[0]
    import  eoldas
    from solve_eoldas_identity import *
    # import the setup methods from solve_eoldas_identity
    import pylab as plt

    # SD of noise
    noise=0.15
    # nominal sampling period
    n_per=7
    # proportion of missing samples
    obs_off=0.33
    # order of differential model (integer)
    model_order=1
    # sgwindow is larger to create larger data gaps
    sgwindow=10

    # set up data for this experiment
    file, ideal_gamma,doys,ndvi_clean,ndvi,qa_flag  = \
		prepare(noise=noise,n_per=n_per,\
			obs_off=obs_off,model_order=model_order,\
			sgwindow=sgwindow)

    # set gamma to thge theoretical value
    gamma = ideal_gamma

    # set op file names
    xfile = 'output/Identity/NDVI_Identity_a.params'
    yfile = 'output/Identity/NDVI_fwd_a.params'

    # initialise options for DA overriding any in config files
    # make sure we use some different output file names to othe scripts
    cmd = 'eoldas ' +  \
                ' --conf=config_files/eoldas_config.conf --conf=config_files/Identity.conf ' + \
                ' --logfile=mylogs/Identity.log ' + \
                ' --calc_posterior_unc ' + \
                ' --parameter.solve=0,1 ' + \
		' --parameter.result.filename=%s '%xfile +\
                ' --parameter.x.default=%f,0.0 '%(gamma) + \
		' --operator.obs.y.result.filename=%s'%yfile +\
                ' --operator.obs.y.state=%s'%file+\
                ' --operator.modelt.rt_model.model_order=%d '%model_order 

    # initialise eoldas
    self = eoldas.eoldas(cmd)
    # solve DA
    self.solve(write=True)
    
    # now pull some data out of the eoldas

    # the 'root' of the DA data structure is in self.solver.root

    root = self.solver.root

    # The state vector data are stored in root.x
    # with ancillary information in root.x_meta
    # so the state vector is e.g. in root.x.state
    # and the names are in root.x_meta.state

    state_names = root.x_meta.state
    state = root.x.state

    # The sd representation of the posterior is in root.x.sd
    # This is all set up in eoldas_Solver.py
    # All storage is of type ParamStorage, an extended
    # dictionary structure. You can explore it interactively
    # with e.g. root.dict().keys() or self.keys() since
    # self here is a straight dictionary
    sd = root.x.sd

    # The full inverse Hessian is in self.solver.Ihessian
    Ihessian = self.solver.Ihessian

    # A mask to reduce this to only the state variables
    # being targeted (solve == 1) is through a call to:
    NDVI_Ihessian = self.solver.unloader(None,Ihessian,M=True)

    # This is of shape (365,365) here.
    # so now lets produce an image of it
    # to visualise the structure
    fig = plt.figure()
    ax = fig.add_subplot(111)
    cax = ax.imshow(NDVI_Ihessian,interpolation='nearest')
    ax.set_title('Posterior Uncertainty matrix for NDVI')
    # Add colorbar, make sure to specify 
    # tick locations to match desired ticklabels
    cbar = fig.colorbar(cax) #), ticks=[-1, 0, 1])
    #cbar.ax.set_yticklabels(['< -1', '0', '> 1'])# vertically oriented colorbar
    # see http://matplotlib.sourceforge.net/plot_directive/mpl_examples/pylab_examples/colorbar_tick_labelling_demo.py
    # save it   
    plt.savefig('output/IHessianNDVI_expt1.png')

 
   
 

    

The comments in the code should be self explanatory and anyone interested in delving much further into the eoldas codes should see the full class documentation. Here, we can see at least how to access the posterior estimate of the state and its uncertainty. We write out the uncertainty to an image using matplotlib (pylab):

images/example1plit, from ../solve_eoldas_identity.py

Now, this is a very interesting figure for understanding how these multiple constraints are interacting. The observation uncertainty is just described by standard deviation, so lies along the leading diagonal of the a priori uncertainty. Further, it only exists where there are data points. The impact of applying the (regularisation) model constraint is to reduce the uncertainty at the observation points as we would expect. We can see a ‘sausage’ pattern from above in this figure quite clearly. The ‘pinch points’ are when the sample points are dense. Where there are large data gaps (from our simulated cloud impacts here) the uncertainty is higher, but it is ‘spread out’ from the leading diagonal by applying temporal covariance. The impact of the filtering is large where the observation impact is low (or non existant). We will see these same effects in many DA experiments.

Running EOLDAS from the command line

An alternative to writing your own python code for the front end is to use eoldas.py, which can be direcly run from the command line.

Help on command line options is available by typing::

eoldas_run.py --help

As an aid to setting up the correct python variables etc, a front end script, eoldas can also be accessed (in the bin directory).

eoldas_run.py --conf=confs/Identity.conf

N.B. Make sure the eoldas_run.py script is in your path. If you install the python packages for a single user, in UNIX it will usually be under ~/.local/bin/. You may want to add that path to your users’ path.

Application to MODIS reflectance data

We can now apply the concepts demonstrated above to real EO data from the MODIS sensors.

This is actually quite trivial to achieve as it involves the same basic configuration file as previously. This time, we will access it from eoldas_run.py.

The MODIS data file format

An example data file then is data/modis_botswana.dat which is MODIS reflectance data for a site in Botswana. The format of this file is ‘BRDF’, which looks like:

BRDF 67 7 645 858.5 469 555 1240 1640 2130  0.003 0.004 0.004 0.015 0.013 0.01 0.006
181 1 62.830002 -83.040001 37.910000 23.219999 0.032900 0.068600 0.018000 0.028400 0.098100 0.097600 0.081900 
182 1 33.980000 99.540001 44.970001 39.290001 0.050500 0.091800 0.030000 0.045300 0.123200 0.137100 0.124800 
184 1 51.660000 99.300003 46.889999 42.349998 0.056200 0.106000 0.029500 0.053800 0.151600 0.154600 0.141400 
185 1 32.880001 -81.900002 40.540001 31.510000 0.039200 0.069600 0.022900 0.034400 0.089800 0.102500 0.092800 
186 1 63.049999 100.470001 48.910000 45.169998 0.058200 0.138600 0.038900 0.055200 0.184600 0.183800 0.143100 
187 1 6.430000 -79.459999 42.090000 35.259998 0.043700 0.074600 0.026400 0.038700 0.103100 0.118400 0.103600 
189 1 22.389999 97.720001 43.790001 38.750000 0.044500 0.088400 0.027100 0.040200 0.129500 0.135800 0.127600 
190 1 57.400002 -82.720001 38.029999 26.480000 0.050700 0.087700 0.024600 0.045200 0.119900 0.118000 0.098100 
191 1 44.040001 99.419998 45.619999 41.990002 0.062800 0.107800 0.035600 0.057100 0.141600 0.154800 0.137600 
192 1 42.689999 -82.820000 39.290001 30.719999 0.041300 0.078800 0.022700 0.036600 0.106400 0.109000 0.096100 
193 1 58.139999 99.800003 47.560001 45.000000 0.063500 0.124200 0.034900 0.060800 0.169700 0.170900 0.143400 
194 1 20.129999 -82.099998 40.730000 34.689999 0.038600 0.073200 0.024000 0.035100 0.106900 0.111700 0.105200 
196 1 8.940000 96.330002 42.330002 38.389999 0.052200 0.090000 0.030300 0.047200 0.118400 0.131700 0.123100 
197 1 62.480000 -83.120003 36.610001 25.559999 0.053800 0.092700 0.025000 0.049200 0.118100 0.120200 0.093200 
198 1 34.730000 99.230003 44.060001 41.840000 0.054500 0.103800 0.030300 0.051400 0.144800 0.147200 0.133300 
199 1 50.549999 -82.849998 37.750000 30.030001 0.044100 0.082200 0.025500 0.041300 0.104800 0.114500 0.093900 
200 1 52.070000 100.010002 45.910000 45.029999 0.074300 0.132000 0.037400 0.069800 0.174200 0.169700 0.143300 
201 1 32.130001 -82.889999 39.080002 34.240002 0.041500 0.077000 0.026200 0.038900 0.071200 0.108400 0.096300 
202 1 63.400002 100.080002 47.880001 48.009998 0.082600 0.155400 0.043600 0.073800 0.191000 0.187000 0.137200 
203 1 5.270000 -76.970001 40.590000 38.189999 0.053000 0.092500 0.026800 0.047500 0.122900 0.121400 0.109300 
205 1 23.219999 97.779999 42.220001 41.840000 0.067000 0.113900 0.034600 0.060000 0.143100 0.148500 0.131800 
206 1 56.959999 -83.139999 35.910000 29.340000 0.065800 0.102000 0.033800 0.062000 0.117800 0.113800 0.087400 

The first line is a header, as with the previous ‘PARAMETERS’ format we saw. The first word on the header must be BRDF (a # can be included for compatitbilty with some other formats).

The second item on the header line is the number of samples in the file (92 here). The third item is the number of wavebands represented (7 here) then this is followed by wavelengths, wavebands, or waveband identifiers.

If a single float number is given, it is assumed to represent a narrow waveband centred around that wavelength. Wavelength is assumed to be in nm. If it is two float numbers connected by a dash (e.g. 45-550) then it is taken to be a tophat function representation of a waveband, where the two number represent the minimum and maximum wavelength respectively. If some other code is found (e.g. B2), it is assumed to be a tag associated with a particular waveband. In this case, any spectral libraries loaded in the configuration files are searched to see if they contain a definition of this waveband. If nothing is found and it cannot be interpreted as a number, it is assigned an arbitrary wavelength (index starting from 0) and should not be interpreted as a wavelength identifier.

The 7 columns following the wavebands specify the standard deviation in the reflectance data. All subsequent lines contain reflectance, as a function of ‘locational’, ‘control’ and spectral information. The first column here defines location (time), columns 2 to 6 are ‘control’ variables: things the observation waries as a function of that are not spectral or locational. In this case, they are (all angles are in degrees):

# a data mask (1 for good), # view zenith angle, # view azimuth angle, # solar zenith angle, # solar azimuth angle.

The final 7 columns give the reflectance in the 7 wavebands for each location defined.

The data span days of year 181 to 272 inclusive. There is a lot of day to day variation in the data, but a clear underlying trend in the Near infrared. We will set ourselves the task of trying to use eoldas to filter the dataset so that the trend becomes more apparent. We will also try to extrapolate from the sample days to the whole year (which is clearly quite a challenge, but instructive for understanding uncertainty).

We will apply the same data assimilation components as previously, i.e. an Identity operator for the NIR reflectance (in other words, we take the reflectance as a prior estimate of what we wish to estimate) and a regularisation filter implemented as a constraint on first order derivatives at lag one day. We will suppose that we know the uncertainty in this model (the temporal constraint) to be expressed as a standard deviation of .002 (1/500), i.e. the root mean squared deviation from one day to the next is expected to be around .002. An examination of the dataset will show that the NIR reflectance roughly increases from around 0.1 to 0.2 over about 100 days, so we could use a gamma value of 1000. However, we don’t want to impose that constraint too strongly, so we will instead use a gamma of 500. We are not trying to fit a trendline here, just to smooth the dataset so that we can detect the underlying behaviour. The result should however not be very sensitive to this guess and can always be explored using approaches such as cross validation.

The configuration file

We have a configuration file set up in config_files/Identity2.conf. This provides a full description of the problem we wish to solve, which is to reduce noise in the NIR observation series.

# An EOLDAS configuration file for
# a simple Identity operator and a regulariser

[parameter]
location = ['time']
limits = [[1,365,1]]
help_limits="define the limits for the required locational information"
names = "gamma_time 858.5".split()
solve = 0,1
datatypes = x

[parameter.result]
filename = output/Identity/MODIS_botswana.params
format = 'PARAMETERS'
help_filename='Set the output state filename'

[parameter.x]
datatype = x
names = $parameter.names
default = 500,0.1
state= data/modis_botswana.dat
help_default="Set the default values of the states"
apply_grid = True
sd = [1.]*len($parameter.names)
bounds =  [[0.000001,1000000],[0,1]]

[general]
doplot=True
help_do_plot='do plotting'
calc_posterior_unc=False

[operator]
modelt.name=DModel_Operator
modelt.datatypes = x
obs.name=Operator
obs.datatypes = x,y

[operator.modelt.x]
names = $parameter.names
sd = [1.0]*len($operator.modelt.x.names)
datatype = x

[operator.modelt.rt_model]
model_order=1
help_model_order='The differential model order'
wraparound=periodic,365
#wraparound=None

[operator.obs.x]
names = $parameter.names[1:]
sd = [1.0]*len($operator.obs.x.names)
help_sd='Set the observation sd'
datatype = x

[operator.obs.y]
control = [mask,vza,vaa,sza,saa]
names = $parameter.names[1:]
sd = [0.015]*len($operator.obs.x.names)
#help_sd="set the sd for the observations"
datatype = y
state = data/modis_botswana.dat
help_state = "Set the y state vector"

[operator.obs.y.result]
filename = output/Identity/Botswana_fwd.params
format = 'PARAMETERS'
help_filename = 'Set the fwd modelling result filename'



In this, we can see the required location field defined in the [parameters] section (i.e. parameters.location=[[1,365,1]]). The control variables are associated with the observation operator, here given as [mask,vza,vaa,sza,saa].

To run eoldas with this configuration then all we need do is type:

eoldas_run.py --conf=config_files/eoldas_config.conf \
  --conf=config_files/Identity2.conf --calc_posterior_unc

This writes the files:

output/Identity/MODIS_botswana.params: state estimation (smoothed reflectance)

 181.000000    1.000000  500.000000    0.095014    0.000000    0.005513
 182.000000    1.000000  500.000000    0.095247    0.000000    0.005255
 183.000000    1.000000  500.000000    0.095541    0.000000    0.005063
 184.000000    1.000000  500.000000    0.095835    0.000000    0.004840
 185.000000    1.000000  500.000000    0.095949    0.000000    0.004669
 186.000000    1.000000  500.000000    0.096530    0.000000    0.004541
 187.000000    1.000000  500.000000    0.096364    0.000000    0.004451
 188.000000    1.000000  500.000000    0.096585    0.000000    0.004396
 189.000000    1.000000  500.000000    0.096806    0.000000    0.004293
 190.000000    1.000000  500.000000    0.097176    0.000000    0.004216
 191.000000    1.000000  500.000000    0.097714    0.000000    0.004161

output/Identity/Botswana_fwd.params: forward modelling of y

#PARAMETERS time mask,vza,vaa,sza,saa 858.5 sd-858.5
 181.000000    0.000000    0.095014    0.000000
 182.000000    0.000000    0.095247    0.000000
 184.000000    0.000000    0.095835    0.000000
 185.000000    0.000000    0.095949    0.000000
 186.000000    0.000000    0.096530    0.000000
 187.000000    0.000000    0.096364    0.000000
 189.000000    0.000000    0.096806    0.000000
 190.000000    0.000000    0.097176    0.000000
 191.000000    0.000000    0.097714    0.000000

as well as the orignal data in output/Identity/Botswana_fwd.params_orig

#PARAMETERS time mask,vza,vaa,sza,saa 858.5 sd-858.5
 181.000000    0.000000    0.068600    0.015000
 182.000000    0.000000    0.091800    0.015000
 184.000000    0.000000    0.106000    0.015000
 185.000000    0.000000    0.069600    0.015000
 186.000000    0.000000    0.138600    0.015000
 187.000000    0.000000    0.074600    0.015000
 189.000000    0.000000    0.088400    0.015000
 190.000000    0.000000    0.087700    0.015000
 191.000000    0.000000    0.107800    0.015000

with appropriate graphics for the state:

_images/MODIS_botswana.params.plot.x.png

and the y data:

_images/Botswana_fwd.params.plot.y2.png _images/Botswana_fwd.params.plot.x.png

The resultant state data are quite instructive: where we have observations, the uncertainty is reduced from 0.015 to around 0.004 (the actual degree of noise reduction depends on the value of gamma used). Where there are no data, the uncertainty grows to around 0.017. It is slightly reduced at the year start/end because of the wrparound condition used here.

In this case, we have used a first order differential constraint with a periodic boundary condition. These are quite important in this case: we only have observations in a limited time window, so using a periodic boundary condition is one way to place some form of constraint at what happens when we have no data. The first order differential model will in essence perform a linear interpolation where there are no data, which is probably appropriate for this case. You can try changing the model order to see what happens. Run e.g.:

eoldas_run.py --conf=config_files/eoldas_config.conf \
  --conf=config_files/Identity2.conf --calc_posterior_unc \
  --operator.modelt.rt_model.model_order=2 \
  --parameter.x.default=5000,0.1 \
  --operator.obs.y.result.filename=output/Identity/Botswana_fwd.params2 \
  --parameter.result.filename=output/Identity/MODIS_botswana.params2

and have a look at output/Identity/MODIS_botswana.params2.plot.x.png:

_images/MODIS_botswana.params2.plot.x.png

The assumed behaviour is quite different to the first order differential constraint outside of the observations. With a high value of gamma, the result is essentially a straight line where there are observations. The influence of the wraparound condition is also clear here.

If we used a lower gamma, we would see some features of using a second order model:

eoldas_run.py --conf=config_files/eoldas_config.conf \
  --conf=config_files/Identity2.conf --calc_posterior_unc \
  --operator.modelt.rt_model.model_order=2 \
  --parameter.x.default=200,0.1  \
  --operator.obs.y.result.filename=output/Identity/Botswana_fwd.params3  \
  --parameter.result.filename=output/Identity/MODIS_botswana.params3
_images/MODIS_botswana.params3.plot.x.png

Where there are observations, the second order difference constraint would have high frequency oscillations. This can be a positive feature or an annoyance: it all depends on how you expect the function to behave. It is not an arbitrary choice: the user has imposed a particular expectation of behaviour here through the model. There clearly needs to be some evidence for choosing one form of model over another, although that is not always straightforward.

One other interesting feature of this result is that the mean estimate is bounded (0.0,1.0) which is reasonable for reflectance data (theoretically, BRF can go above 1, but this is rarely observed). This means that the mean value (the red line) is forced to stay above 0 (days 100 to 150) even though the apparent trend from the observations might otherwise make it go below zero. This condition is imposed in the configuration by the line:

bounds =  [[0.000001,1000000],[0,1]]

One further thing to note about this result is that when the uncertainty is so large (as is the case here when we are extrapolating) that the confidence interval spans more than the entire data range (0,1) the mean and standard deviation described in the Gaussian statistics are a little meaningless or at least should be interpreted with caution. It might in this case be better to describe the data (in the extrapolated region) as simply being somewhere between the physical limits.

A final comment is that actually, a large amount of the departure of the signal from our smooth trajectory, the high frequency variation that we have treated as noise here, can most likely be described by considering the physics os scattering by the surface (there are, on the whole, BRDF effects). We will return to this later.

As well as raising a few issues with regard to model selection, we have demonstrated further use of the eoldas command line for quickly changing some terms in an experiment. However, it is all very well in showing that eoldas can ingest satellite data and apply constraints to provide an estimate of state variables, but the state variables here (reflectance) do not directly help us monitor the properties of the land surface or link to process models.

To do that using EO data, we generally need more complex obeseravtion operators. These are dealt with in the next section.

«  A tutorial guide of the EOLDAS prototype software   ::   Contents   ::   Radiative transfer modelling for optical remote sensing  »