User:Jbeezley

From CCM
Jump to: navigation, search

Contents

Using osimorph classes

This page is meant to be a basic walkthrough of the features of the osimorph classes. It assumes knowledge of basic Matlab programming. For more detailed information see the code and comments therein.

variable

  • describes a model scalar field (such as temperature)
  • does not contain any actual data, just meta information
  • properties:
    • char type[]: the data type of the variable (i.e. 'double', 'single', 'int32')
    • dimension dims{}: the grid over which the variable is defined
    • containers.Map attrs: attributes which control how the variable is treated (perturbation size, whether to morph the variable, how to plot, etc.)
  • static properties, which define attributes that can be customized (i.e. var.setattr(variable.c_pert_add,10)):
    • c_pert_add: magnitude of additive perturbations
    • c_pert_mul: magnitude of multiplicative perturbations
    • c_pert_alpha: smoothness of perturbations
    • pert_morphing: logical, do we perturb this variable spatially?
    • plot_type: how to display the variable ('mesh', 'image', 'none',etc.) see plot methods
    • vector: if this variable is a component of a vector, set this attribute to the name of the vector, otherwise, don't use
    • ivector: integer, component index of the vector (1 for 'x', 2 for 'y', 3 for 'z')
    • do_assim: if this variable is to be included in the assimilation state
    • do_morph: if this variable is to be transformed by morphing
    • datavar: if this variable is used as a synthetic observation
    • registration_var: if this variable is used for image registration

Examples:

% set up search path
cd osimorph/matlab
startup

% define a 10 x 11 grid
dims={dimension('x',10),dimension('y',11)};

% create the variable with type 'double' and default attributes
var=variable('varA',dims,'double');

% fill an array the size of the variable with random data
a=randn(var.size);

% write the data to a file, 'example.nc'
var.write('example.nc',a);

% plot the variable from the file
var.plot('example.nc');

% plot same figure, but this time save the image to a file
% called wiki_example1_var.png as well
var.plot('example.nc',[],'wiki_example1_var.png');

% change the plot type, and plot again
var.setattr(variable.plot_type,'image');
var.plot('example.nc',[],'wiki_example1_img.png');

% tag this as a registration variable
var.setattr(variable.registration_var,attribute.true);

% get all attributes
attrs=var.getattrs();

% randomly perturb the data in 'example.nc' by twice the standard magnitude
var.randomPert('example.nc',2)

% read the new data and plot the difference from the old data
% change the default title to describe the content of the figure
b=var.read('example.nc');
var.plot_data(a-b,'variable perturbation (a-b)','wiki_example1_diff.png');

variable_list

  • describes any list of variable objects
  • contains helper methods for acting on several variables at once

Examples:

% create a variable_list object from scratch
dims={dimension('x',10),dimension('y',11)};
vlist=variable_list({variable('varA',dims,'double'),variable('varB',dims,'double')});

% create a variable_list object from all variables in a netcdf file (for example a restart file from wrf-fire;
% copy such file in the current directory first as file tmp.nc)
vlist=variable_list.getFromFile('tmp.nc');

% list the names of all variables
vlist.getVarNames()

% extract variable objects from the list
tign_g=vlist.getVar('TIGN_G');
grnhfx=vlist.getVar('GRNHFX');
u=vlist.getVar('UF');
v=vlist.getVar('VF');

% adjust properties of the variables
grnhfx.setattr(variable.registration_var,attribute.true);
grnhfx.setattr(variable.plot_type,'image');
u.setattr(variable.vector,'wind');
u.setattr(variable.ivector,1);
v.setattr(variable.vector,'wind');
v.setattr(variable.ivector,2);

% create a new variable_list from these variables
vlist=variable_list({grnhfx,tign_g,u,v});

% read all of these variables from the file into a cell array
c=vlist.read('tmp.nc');

% plot all variables
vlist.plot('tmp.nc');

% write data to a new file
vlist.write('tmp1.nc',c);

% perturb new file
vlist.randomPert('tmp1.nc'); % standard perturbation
vlist.randomMorph('tmp1.nc'); % morphing perturbation

ensemble

  • a subclass of variable_list
  • contains a list of files representing the ensemble
  • most methods are overloaded from variable_list to eliminate the filename arguments

Example

% create a variable list from a wrf-fire output file, limited to the 4 variables given
vlist=variable_list.getFromFile('tmp.nc',{'TIGN_G','GRNHFX','UF','VF'});

% create an ensemble object with file names formatted as 'file_%03i.nc', with 3 members
ens=ensemble(3,vlist.variables,'file_%03i.nc');

% read the data from the initial file 'tmp.nc' and write to ensemble files unmodified to initialize
c=vlist.read('tmp.nc');
c=repmat(c,1,3);
ens.write(c);

% randomly perturb the ensemble members
ens.randomPert();
ens.randomMorph();

% get an ensemble matrix useful for input into enkf
c=ens.read();
X=ens.cell2matrix(c);

% convert an ensemble matrix back to a cell array for writing back to files
c=ens.matrix2cell(X);
ens.write(c);

models

A "model" in osimorph is the component that manages the initializing and advancing of the ensemble members. All model classes should be subclassed from the abstract class model, which defines a four standard methods that the model implementations must override. Because every model has its own list of variables that it manipulates, the class variable_list is a subclass of the abstract model; all variable_list methods are available in the implementation of the model.

Methods required of model implementations are as follows:

  • initialize(files): Initialize the model in the files given by the cell array of strings. After initialization the files should exist and contain all meta information and a model state at the initial time.
  • getInitialTime(): Get the starting time of the model. At the moment, the return class of this function is not necessarily numeric. It can be any class that supports logical operators `==' and `~='.
  • advance(time,files): Advance the model in `files' to the prescribed time. This files will already exist. It is up to the model to figure out what the current time of the files is.
  • getVariables(): Return a cell array containing variable objects for each variable in the model state. This method can simply return it own variables property.

Examples (see implementations in the classes directory):

  • model_identity: A "dummy" model that uses a single netCDF file as a template model state. Initialization of a new simulation simply copies the original file and advancement does nothing.
  • model_bump: Similar to model_identity except it constructs a single variable model state from scratch. The variable is a single gaussian "bump". The construct takes an additional optional argument isdata. If this argument is passed (it doesn't matter what it contains), then the bump will be in a different location. This class is design to be a simple test case for morphing.
  • model_epi: A more complex example of a model that calls an external program (the epidemic model) to initialize and advance the state.
% create a model_bump object
m=model_bump();

% initialize two files
m.initialize({'bump1.nc','bump2.nc'});

% perturb one of the files
m.randomPert({'bump1.nc'});
m.randomMorph({'bump1.nc'});

% plot the contents of each file
m.plot('bump1.nc');
m.plot('bump2.nc');

morphing

  • Matlab interface to morphing transformation binaries
  • needs some improvement to allow for customization of namelist parameters, to be done in future

Example:

% create a dummy model class from the wrf-fire restart file
mod=model_identity('tmp.nc');

% modify one of the variables to tag it as the registration variable
mod.getVar('GRNHFX').setattr(variable.registration_var,attribute.true);

% create an ensemble object from the model variables
ens=ensemble(3,mod.getVariables(),'file_%03i.nc');

% initialize the files and perturb
mod.initialize(ens.files());
ens.randomPert();
ens.randomMorph();

% create a morphing object
morph=morphing(mod);

% create a transformed ensemble
trans=morph.forward_trans(ens);

% invert the transformation
ens=morph.inverse_trans(trans);

observations

Observations in osimorph handle all kinds of data both real and synthetic. In general, observation classes can only be used with specific model classes. As with models, observation implementations should be subclassed from the observation abstract class. The data contained in an observation are split into multiple variables much like in a model. The data in each variable is an array of scalars representing a single kind of measurement, like an array of temperature sensors, or an image of heat flux. While observations are subclasses of variable_list, the components of the list are a special variables called data_variable. A data variable differs from a standard model variable in that it also contains variance information.

Implementations must override the following methods:

  • getVariance(): Return a cell array of data variances for each data variable contained in the observation.
  • obsFunction(): Return a function handle that computes synthetic data from a model state.
  • obsEnsemble(ens): Generates synthetic data from the given ensemble.
  • obsMatrix(): Returns an observation function in sparse matrix form as used in the EnKF formulation.
  • getNextObsTime(): Returns the time of the next observation.
  • deleteObs(t): Deletes the observation at time t.

Simplified observation subclasses exist to make it easier to develop basic simulations using synthetically generated data. This are described below.

observations on a subset of the model

The class simple_observation defines direct observations on one or more model variables. The constructor for this class has the following calling syntax.

simple_observation(model,obsvarnames,variance,obstime,obsfile)
  • model: a model object
  • obsvarnames: a cell array of variable names that are observed from the model
  • variance: a cell array of variance matrices for each observed variable
  • obstime: a vector containing the times at which the observations are made
  • obsfile: the name of a netcdf file containing observed data (typically a model output file)

synthetic observations generated from a model

The class synthetic_observation is much like simple_observation except it also handles the advancement of the synthetic data in time. The calling syntax of this class is as follows.

synthetic_observation(model,obsvarnames,variance,obstime,cpert,cpert_morph)
  • model: a model object
  • obsvarnames: a cell array of variable names that are observed from the model
  • variance: a cell array of variance matrices for each observed variable
  • obstime: a vector containing the times at which the observations are made
  • cpert: magnitude of additive perturbation for the initial model state
  • cpert_morph: magnitude of the spatial perturbation for the initial model state

Data assimilation with osimorph classes

Utilizing the model and observation classes and associated methods, implementing data assimilation is relatively easy. The abstract class assimilation defines the basic interface. The assimilation constructor takes an ensemble object and an observation object and stores them as properties. Data assimilation implementations override the update method, which takes no arguments and returns then analysis ensemble. The following assimilation algorithms are implemented.

  • enkf: classical perturbed data enkf algorithm
  • fft_enkf: enkf with diagonal fft covariance approximation
  • wavelet_enkf: enkf with diagonal wavelet covariance approximation

For fft_enkf and wavelet_enkf, the majority of the implementation is identical so they are both subclassed from spectral_enkf, which contains the common code.

Single cycle examples

Following are some examples of using prebuilt classes to create some simple data assimilation experiments. These examples build on all of the classes introduced above to produce a working simulation from start to finish.

Model bump

The bump model consists of a single variable that does not change in time. The variable consists of a single gaussian "bump" near the middle of the domain. The constructor allows the user to flag the object as data in which case the bump will be in a different spot on the domain. This feature is used to test assimilation of coherent features with large spatial error.

% model object for ensemble members
m_ens=model_bump();

% model object for data
m_data=model_bump(1);

% create and initialize 5 ensemble members
ens=ensemble(5,m_ens.getVariables());
m_ens.initialize(ens.files());

% randomly perturb ensemble with 2x standard morphing size 4x standard additive size
ens.randomPert(4)
ens.randomMorph(2)

% create data with uniform variance of 0.1
data=synthetic_observation(m_data,m_data.getVarNames(),0.1,[0]);

% plot ensemble and data (these are saved to files automatically)
ens.plot()
data.plot()

% set up morphing (use variance .001 for morphing function)
morph=morphing(m_ens);
morph.mxvariance=.001;
morph.myvariance=.001;

% do the morphing transformation of the ensemble and data
[m_ens,m_data]=morph.forward_trans(ens,data);

% plot transformed ensemble and data, set up format string to save to a different file name
m_ens.plot('morphed_mem%i_%%s','morphed_mem%i_%%s.png');
m_data.plot('morphed_data_%s','morphed_data_%s.png');

% do enkf on both the standard ensemble and the morphed ensemble
enkf_assim=enkf(ens,data);
enkf_assim.update()
m_enkf_assim=enkf(m_ens,m_data);
m_enkf_assim.update();

% do inverse transformation on the morphed analysis
m_ens=morph.inverse_trans(m_ens);

% plot analysis ensemble for both assimilations
ens.plot('analysis_mem%i_%%s','analysis_mem%i_%%s.png');
m_ens.plot('morphed_analysis_mem%i_%%s','morphed_analysis_mem%i_%%s.png');

Automated multi-cycle assimilation

Because the interfaces for assimilation methods and models are identical, a class is available that automates the creation of a simulation and generation of plots. This class is called simulation. This classes constructor has the following syntax. simulation(num_mem,model,obs,trans,assim)

  • num_mem: the number of ensemble members to use
  • model: an instantiated model object (subclassed from model)
  • obs: an instantiated observation object (subclassed from observation)
  • trans: an instantiation of a subclass of transformation (i.e. morphing)
  • assim: a function handle for the constructor of a subclass of assimilation (i.e. @enkf or @wavelet_enkf)

Calling the run method of the simulation will automatically initialize and perturb an ensemble, run the assimilation, and create png images of all variables at every step. Control over the parameters of the simulation can be modified by adjusting the properties of the input class objects and attributes of the model variables. An example of using this class is given below.

% create an bump model instance
mbump=model_bump();
mdata=model_bump(1);

% variable to use for observations
data_var='RAIN_DIF';

% variance of observation
variance=.1;

% variance of morphing function
morph_variance=.001;

% simulation times at which observations occur
observation_times=[10 20];

% construct observation instance
obs=synthetic_observation(mdata,{data_var},variance,observation_times);

% morphing class object
trans=morphing(mbump);
trans.mxvariance=morph_variance;
trans.myvariance=morph_variance;

% trans=@transformation; % for no morphing

% create simulation instance
sim=simulation(3,mbump,obs,trans,@enkf);

% run simulation
sim.run();
Personal tools
Namespaces

Variants
Actions
Navigation
Toolbox