Taylor diagram

 


Taylor diagrams plot

SUB-PROGRAM




Need the Python source code ? Click here.

Description of sub-program code:
<*>=
<Module Documentation>
<Import library modules>
<Constant definitions>
<Sub-routines of this module>
<Main routine>


#=============================================================================
#=============================================================================

<Module Documentation>= (<-U)

"""
    Taylor diagrams set for MODEL-DATA comparison:
    ==============================================

    Plot one or multiple Taylor diagrams based on comparison statistics
    between a list of model output and a data set.

    This program has three modes depending on a command line option:
        Mode 1: for one variable, plot multiple Taylor diagrams,
                one per statistical component (the default mode)
        Mode 2: for one variable, plot all statistical components
                on one Taylor diagram.
        Mode 3: for multiple variables, plot the total space-time components
                on one Taylor diagram.

    For one variable (e.g. Sea Surface Temperature), comparison statistics
    between observation data and models are read from a netCDF file.
    This file contains, for each model, comparison statistics for all
    statistical components.

    The component list is the following:

        - 1. Total space-time field
        - 2. Zonal annual mean
        - 3. Zonal monthly mean
        - 4. Zonal monthly anomaly: Monthly anomaly to zonal annual mean
        - 5. Annual mean map
        - 6. Monthly map anomaly:   Monthly anomaly to annual mean map
        - 7. Longitudinal anomaly:  Spatial anomaly to zonal monthly mean
        - 8. Annual mean longitudinal anomaly:  Annual mean longitudinal anomaly to zonal mean

    Additional points may be shown on the diagrams along with models, like
        - all models mean or median
        - additionnal data-set

    Each model is represented on the diagrams as
        - either a one-or-two letter code, optionally with a specific color.
        - or a coloured dot with a one-or-two letter code. The dot color is
          taken from a color scale varying from purple to blue and indicates
          the model bias.

    By default, model names are extracted from input NetCDF file.
    Code-letters are derived from these names and predefined colors are
    automatically assigned to them.

    An optional parameter file allows the user to change the model list
    and specify their attributes. (see  section below)

    Output files:
    ------------
    Each Taylor diagram is saved in a file. Default format is PostScript.
    Alternate option is GIF. Respective file name extension is '.ps' or '.gif'

    When running in mode 1, each file name is after the component short name
    (first listed above) with all white space replaced by underscores.
    For example:  Total_space-time_field

    Directory where these files are created is by default the current working
    directory. This may be specified differently, using a command line option.
    If the specified directory does not exist, it is created. If output files
    already exist in that directory, they are overwritten.

    When running in mode 2, output image file name is after the variable name
    from input NetCDF file and is in current working directory.
    Path and name may be specified differently using a command line option.

    When running in mode 3, output image file name is after the FIRST variable name
    from (FIRST) input NetCDF file and is in current working directory.
    Path and name may be specified differently using a command line option.
"""

#=============================================================================
#=============================================================================

<Import library modules>= (<-U)

import Numeric
import cdms, vcs
import string, sys, os
import getopt


#=============================================================================
#=============================================================================

<Constant definitions>= (<-U)

model_list = []
nb_models = 0
all_dictionary = {}

component_names = (                             \
'Total space-time field',                       \
'Zonal annual mean',                            \
'Zonal monthly mean',                           \
'Monthly anomaly to zonal annual mean',         \
'Annual mean map',                              \
'Monthly anomaly to annual mean map',           \
'Spatial anomaly to zonal monthly mean',        \
'Annual mean spatial anomaly to zonal mean')

component_short_names = (            \
'Total space-time field',            \
'Zonal annual mean',                 \
'Zonal monthly mean',                \
'Zonal monthly anomaly',             \
'Annual mean map',                   \
'Monthly map anomaly',               \
'Longitudinal anomaly',              \
'Annual mean longitudinal anomaly')

nb_components = len (component_names)

# Color map varying from dark blue to light grey to dark red
# defined as a dictionnary (only entries in the range [6, 239] are defined)

red_blue_map_dictionary = { \
6 : [0, 0, 44] , \
7 : [0, 0, 44] , \
8 : [0, 0, 44] , \
9 : [0, 0, 44] , \
10 : [0, 0, 47] , \
11 : [0, 0, 47] , \
12 : [0, 0, 47] , \
13 : [0, 0, 52] , \
14 : [0, 0, 52] , \
15 : [0, 0, 52] , \
16 : [0, 0, 55] , \
17 : [0, 0, 55] , \
18 : [0, 0, 55] , \
19 : [0, 0, 55] , \
20 : [0, 0, 55] , \
21 : [0, 0, 58] , \
22 : [0, 0, 58] , \
23 : [0, 0, 62] , \
24 : [0, 0, 62] , \
25 : [0, 0, 62] , \
26 : [0, 0, 62] , \
27 : [0, 0, 62] , \
28 : [0, 0, 65] , \
29 : [0, 0, 68] , \
30 : [0, 0, 68] , \
31 : [0, 0, 68] , \
32 : [0, 0, 68] , \
33 : [0, 0, 68] , \
34 : [0, 0, 71] , \
35 : [0, 0, 71] , \
36 : [0, 0, 71] , \
37 : [0, 0, 74] , \
38 : [0, 0, 74] , \
39 : [0, 0, 74] , \
40 : [0, 0, 74] , \
41 : [0, 0, 78] , \
42 : [0, 0, 78] , \
43 : [0, 0, 78] , \
44 : [0, 0, 81] , \
45 : [0, 0, 81] , \
46 : [0, 0, 81] , \
47 : [0, 0, 84] , \
48 : [0, 0, 84] , \
49 : [0, 0, 84] , \
50 : [0, 0, 87] , \
51 : [0, 0, 87] , \
52 : [0, 0, 87] , \
53 : [0, 0, 87] , \
54 : [0, 0, 90] , \
55 : [0, 0, 90] , \
56 : [0, 0, 90] , \
57 : [0, 0, 94] , \
58 : [0, 0, 94] , \
59 : [0, 0, 94] , \
60 : [0, 0, 97] , \
61 : [0, 0, 97] , \
62 : [0, 0, 97] , \
63 : [0, 0, 100] , \
64 : [0, 0, 100] , \
65 : [0, 0, 100] , \
66 : [0, 1, 100] , \
67 : [3, 3, 100] , \
68 : [3, 4, 100] , \
69 : [6, 6, 100] , \
70 : [9, 7, 97] , \
71 : [9, 9, 97] , \
72 : [9, 11, 97] , \
73 : [12, 12, 97] , \
74 : [15, 14, 97] , \
75 : [15, 15, 97] , \
76 : [15, 17, 97] , \
77 : [19, 19, 97] , \
78 : [22, 20, 97] , \
79 : [22, 22, 97] , \
80 : [22, 23, 97] , \
81 : [22, 23, 97] , \
82 : [25, 25, 97] , \
83 : [28, 27, 97] , \
84 : [28, 28, 94] , \
85 : [28, 30, 94] , \
86 : [31, 31, 94] , \
87 : [35, 33, 94] , \
88 : [35, 35, 94] , \
89 : [35, 36, 94] , \
90 : [38, 38, 94] , \
91 : [38, 39, 94] , \
92 : [41, 41, 94] , \
93 : [41, 43, 94] , \
94 : [44, 44, 94] , \
95 : [44, 46, 94] , \
96 : [47, 47, 94] , \
97 : [47, 49, 94] , \
98 : [47, 49, 90] , \
99 : [52, 50, 90] , \
100 : [52, 52, 90] , \
101 : [55, 54, 90] , \
102 : [55, 55, 90] , \
103 : [55, 57, 90] , \
104 : [58, 58, 90] , \
105 : [62, 60, 90] , \
106 : [62, 62, 90] , \
107 : [62, 63, 90] , \
108 : [65, 65, 90] , \
109 : [68, 66, 90] , \
110 : [68, 68, 90] , \
111 : [68, 70, 90] , \
112 : [71, 71, 90] , \
113 : [71, 73, 87] , \
114 : [74, 74, 87] , \
115 : [74, 74, 87] , \
116 : [74, 76, 87] , \
117 : [78, 78, 87] , \
118 : [78, 79, 87] , \
119 : [81, 81, 87] , \
120 : [81, 82, 87] , \
121 : [84, 84, 87] , \
122 : [84, 86, 87] , \
123 : [87, 87, 87] , \
124 : [87, 86, 84] , \
125 : [87, 84, 84] , \
126 : [87, 82, 81] , \
127 : [87, 81, 81] , \
128 : [87, 79, 78] , \
129 : [87, 78, 78] , \
130 : [87, 76, 74] , \
131 : [87, 74, 74] , \
132 : [87, 73, 71] , \
133 : [87, 71, 71] , \
134 : [90, 70, 68] , \
135 : [90, 68, 68] , \
136 : [90, 68, 68] , \
137 : [90, 66, 68] , \
138 : [90, 65, 65] , \
139 : [90, 63, 62] , \
140 : [90, 62, 62] , \
141 : [90, 60, 62] , \
142 : [90, 58, 58] , \
143 : [90, 57, 55] , \
144 : [90, 55, 55] , \
145 : [90, 54, 55] , \
146 : [90, 52, 52] , \
147 : [90, 50, 52] , \
148 : [94, 49, 47] , \
149 : [94, 47, 47] , \
150 : [94, 46, 44] , \
151 : [94, 44, 44] , \
152 : [94, 43, 41] , \
153 : [94, 41, 41] , \
154 : [94, 39, 38] , \
155 : [94, 38, 38] , \
156 : [94, 36, 35] , \
157 : [94, 35, 35] , \
158 : [94, 33, 35] , \
159 : [94, 33, 35] , \
160 : [94, 31, 31] , \
161 : [94, 30, 28] , \
162 : [97, 28, 28] , \
163 : [97, 27, 28] , \
164 : [97, 25, 25] , \
165 : [97, 23, 22] , \
166 : [97, 22, 22] , \
167 : [97, 20, 22] , \
168 : [97, 19, 19] , \
169 : [97, 17, 15] , \
170 : [97, 15, 15] , \
171 : [97, 14, 15] , \
172 : [97, 12, 12] , \
173 : [97, 11, 9] , \
174 : [97, 9, 9] , \
175 : [97, 7, 9] , \
176 : [100, 6, 6] , \
177 : [100, 4, 3] , \
178 : [100, 3, 3] , \
179 : [100, 1, 0] , \
180 : [100, 0, 0] , \
181 : [100, 0, 0] , \
182 : [100, 0, 0] , \
183 : [97, 0, 0] , \
184 : [97, 0, 0] , \
185 : [97, 0, 0] , \
186 : [94, 0, 0] , \
187 : [94, 0, 0] , \
188 : [94, 0, 0] , \
189 : [94, 0, 0] , \
190 : [90, 0, 0] , \
191 : [90, 0, 0] , \
192 : [90, 0, 0] , \
193 : [87, 0, 0] , \
194 : [87, 0, 0] , \
195 : [87, 0, 0] , \
196 : [87, 0, 0] , \
197 : [84, 0, 0] , \
198 : [84, 0, 0] , \
199 : [84, 0, 0] , \
200 : [81, 0, 0] , \
201 : [81, 0, 0] , \
202 : [81, 0, 0] , \
203 : [78, 0, 0] , \
204 : [78, 0, 0] , \
205 : [78, 0, 0] , \
206 : [78, 0, 0] , \
207 : [74, 0, 0] , \
208 : [74, 0, 0] , \
209 : [74, 0, 0] , \
210 : [71, 0, 0] , \
211 : [71, 0, 0] , \
212 : [71, 0, 0] , \
213 : [71, 0, 0] , \
214 : [68, 0, 0] , \
215 : [68, 0, 0] , \
216 : [68, 0, 0] , \
217 : [68, 0, 0] , \
218 : [68, 0, 0] , \
219 : [65, 0, 0] , \
220 : [65, 0, 0] , \
221 : [62, 0, 0] , \
222 : [62, 0, 0] , \
223 : [62, 0, 0] , \
224 : [62, 0, 0] , \
225 : [62, 0, 0] , \
226 : [58, 0, 0] , \
227 : [55, 0, 0] , \
228 : [55, 0, 0] , \
229 : [55, 0, 0] , \
230 : [55, 0, 0] , \
231 : [55, 0, 0] , \
232 : [55, 0, 0] , \
233 : [52, 0, 0] , \
234 : [52, 0, 0] , \
235 : [52, 0, 0] , \
236 : [47, 0, 0] , \
237 : [47, 0, 0] , \
238 : [47, 0, 0] , \
239 : [44, 0, 0] }


#=============================================================================
#=============================================================================

<Sub-routines of this module>= (<-U)

<Function usage>
<Function plot_taylor_diagram>
<Function plot_color_bar>
<Function determine_model_list_from_file>
<Function read_needed_statistics_from_file>
<Class reference_std_dev>

#=============================================================================
#=============================================================================

<Function usage>= (<-U)

def usage():
    """
    Function usage
    --------------

    Print out program usage.
    Return nothing
    """
    exec_name = os.path.basename (sys.argv[0])
    print __doc__
    print """
    Usage: %s [Options] [stats-filename]

       or: %s (--mix_variables|-v) [stats-filename stats-filename2 ...]

        <stats-filename> and <stats-filename2> are NetCDF input file names.
        If none is given, <stats-filename> defaults to 'model_vs_data_comp_stats.nc'
        in current working directory.

    Options are:

    *** General options ***

     --help
     -h     print out this help

     --output=<pathname>
     -o <pathname>
            Depending on whether single or multiple image file to be output,
            specify output filename or output directory pathname.
            If multiple output files, file names are fixed (not an option).

            In any case, default directory is current working directory.
            For single output file, default filename depends on chosen mode
            option (see below).

     --gif
     -g     Output file format is GIF. File extension is '.gif'

     --portrait
     -p     Output image orientation is 'portrait'. (Default is 'landscape')

     --models=<model-list>
     -m <model-list>
            Comma-separated list of all models/composite-var to show
            on the diagram. Composite variables may be 'Median' or 'Mean',
            for example.
            By default, all variables present in the input NetCDF file are shown.

            Listed names must be consistent with variable ids in input NetCDF file.
            That is, when <name> appears in the list, <name>_Correlation and
            name>_StandardDeviation must exist as variables in the input file.

            This option is incompatible with --param-file or -f option below.

     --param-file=<filename>
     -f <filename>
            Specify a model parameter file that contains python statements
            defining model list and plot parameters for each models.
            Listed models must be consistent with variable ids in input file.
            See paragraph 'Model parameter file' below for more details.

            This option is incompatible with --models or -m option above.

            If neither model list nor model parameter file is specified on the
            command line, the program looks for a parameter file named
            'model_dictionary.py' in current working directory.

     --symbol-size=<size>
     -i=<size>
            Adjust size of points and associated symbols on the diagrams.
            <size> is an integer value in the range [1, 100].  Default is 50.

     --title=<title>
     -t <title>
            Add a title to the top of each Taylor diagram

     --units=<units>
     -u <units>
            Add the variable units on the standard deviation axis.


     *** Mode options ***

     --superimpose
     -s     Selector for mode 2. All components are superimposed on the same
            Taylor diagram. One different color of standard palette is used
            for each component. Only one image file is output.

            Default is mode 1 where each component is plot separately
            on a Taylor diagram.

     --mix-variables
     -v     Selector for mode 3. Main components (Total space-time field)
            of all variables are on the same Taylor diagram.
            One different color of standard palette is used for each variable.
            Only one image file is output.


     *** Miscelaneous (mode dependendant) options ***

     --components=<component-list>                 (valid only for mode 1 and 2)
     -C <component-list>
            Choose a subset of components to plot.
            <component-list> is a list of numbers corresponding to ranks
            in available component list above (from 1 to 8).
            List order is meaningless.

     --bias                                              (valid only for mode 1)
     -b     For the main statistical component, models are represented on the
            diagram with a dot whose color is function of model bias.
            By default, models are represented with only a coloured letter-code.

     --color                                             (valid only for mode 1)
     -c     If no option --bias or -b is used, represent models on the diagrams
            with coloured code-letter instead of black code-letter.

     --discrete                                          (valid only for mode 1)
     -d     When using --bias or -b option, bias is represented using a limited
            number of colors instead of a continuous color scale.

     --red-blue                                          (valid only for mode 1)
     -r     When using --bias or -b option, bias is represented using a standard
            blue to red palette centered on value 0.


    Model parameter file:

    A model parameter file is an optional python script that lists all models and
        composite vars used for comparison. It describes also, for each model/
        composite-var, the letter-code and the color that represents this model/var
        on the Taylor diagram.

    By default, the program looks for that file under the name
        'model_dictionary.py' in current working directory. User may specify
        a different one using option '-f'.

    File structure:
        It defines two Python variables named: 'model_dictionary' (mandatory, for models)
        and 'other_dictionary' (optional, for additionnal variables).

        These variables are Python dictionnaries with,
        - as keys: model or other variable names
        - as values: inner two-item dictionnaries defining
                     code-letter and color for each model/variable

        Code-letter is a string of one or more letter(s), typically 1 or 2.
        Color is a number ranging from 241 to 255. It is an index
        into the default color palette.
            color: 241   Black
            color: 242   Red
            color: 243   Green
            color: 244   Blue
            color: 245   Yellow
            color: 246   Cyan
            color: 247   Magenta
            color: 248   Orange
            color: 249   Marroon
            color: 250   Navy
            color: 251   Brown
            color: 252   Grey
            color: 253   Pale green
            color: 254   Rose
            color: 255   Light blue

        The color parameter is relevant only when using mode 1 (default)
        whith -c or --color option.

    Example:
         Let's say we deal with two models 'Model_one' and 'Model_two'
         and one additionnal data-set named 'Data'

         'model_dictionary' variable will be defined as the following:
            model_dictionary  = {\\
                "Model_one":      { "code_letter": "A",  "color": 241 },  \\
                "Model_two":      { "code_letter": "2",  "color": 242 },  \\
                }

         'other_dictionary' variable will be defined as the following:
            other_dictionary  = {\\
                "Data":     { "code_letter": "D",  "color": 244 },  \\
                }

    """ % (exec_name, exec_name)


#=============================================================================
#=============================================================================

<Function plot_taylor_diagram>= (<-U)

def plot_taylor_diagram (reference, data, data_name_list, color_list, bias_list, discrete_option, red_blue_option, \
    skill_score_option, max_std_dev, quadrans, ref_color, component_name):

    """
    Function plot_taylor_diagram
    ----------------------------

    Plot one Taylor diagram for several models and one statistical component.
    Models are represented with predefined letter-codes and colors.

    Call:
        plot_taylor_diagram (reference, data, data_name_list, color_list, bias_list, discrete_option, red_blue_option, \\
        skill_score_option, quadrans, max_std_dev, ref_color, component_name)

    Input:
        reference,:        the reference Std Deviation
        data:              Model standard deviations and correlations
        data_name_list:    Model letter-code list
        color_list:        Model color list
        bias_list:         (Optional) list of bias for each variable, otherwise None
        discrete_option:   1 if bias color bar is made of discrete colors, otherwise 0
        red_blue_option:   1 if color palette is standard blue-to-red centered on 0.
        skill_score_option 1 if to draw skill score isolines
        max_std_dev:       Maximum on standard deviation axis (position of outter circle)                                
        quadrans:          Number of quadrans (1 or 2)
        ref_color:         Color of reference circle
        component_name:    Name of statistical component (title of Taylor diagram)
    """

    # Cast data to double: to avoid a bug in Taylor diagram
    data = data.astype (Numeric.Float64)

    nb_data = len (data)

    # Create if necessary a new template
    new_template_name = 'new_'+ data_name_list[0]
    if new_template_name in canvas.listelements ("template"):
        template = canvas.gettemplate (new_template_name)
    else:
        template = canvas.createtemplate (new_template_name)
    template.data.x1 = 0.01
    template.data.x2 = 0.99
    template.data.y1 = 0.01
    template.data.y2 = 0.99

    # Create if necessary a Taylordiagram graphic method
    new_td_name = 'new_'+ string.replace (component_name, " ", "_")
    try:
        td = canvas.createtaylordiagram (new_td_name)
    except:
        td = canvas.gettaylordiagram (new_td_name)

    # Define the reference Std Deviation
    td.referencevalue = reference
    # and the color of circle associated to it
    td.referencecolor = ref_color
    # and the number of quadrans
    td.quadrans = quadrans

    # Specify the outter circle radius
    if max_std_dev:
        td.max = max_std_dev
    else:
        # Determine the maximum of both reference and data standard deviations
        td.max = max (max (data[:,0]), reference) * 1.2

    td.Marker.id = data_name_list

    #
    # Define the markers and ids for each point
    #

    if bias_list != []:
        # Each model/variable is represented on the diagrams as
        # a coloured dot with a one-or-two letter code. The dot color is
        # taken from a color scale (continuous or discrete)

        # Specify marker type and size
        td.Marker.symbol = ["dot"]*nb_data
        size = int (10 * symbol_size / 50.)
        td.Marker.size = [size]*nb_data

        #
        # Specify marker colors:
        # ---------------------

        # Determine range of bias values

        bias_max = max (bias_list)
        bias_min = min (bias_list)
        # actual_range = bias_max - bias_min
        # bias_max += actual_range * 0.1      # add 10% to max
        # bias_min -= actual_range * 0.1      # subtract 10% to min

        if red_blue_option:
            # In this case, the palette is centered on value 0.
            # compute maximum absolute value
            abs_bias_max = max (abs (bias_max), abs (bias_min))
            bias_max = abs_bias_max
            bias_min = - abs_bias_max

        if discrete_option:
            # Determine bias levels and associated colors
            levels = vcs.mkscale (bias_min, bias_max)
            nb_levels = len (levels)

            # Determine colors for each interval between levels
            color_step = (239 - 6) / (nb_levels - 2)
            level_color_list = range (6, 240, color_step)

            # For each point, convert from bias to color (index in colormap)
            level_max = levels[-1]
            level_min = levels[0]
            total_range = level_max - level_min

            color_list = []
            for bias in bias_list:
                # Convert from bias to index in level_color_list
                index =  int ( (bias - level_min) * (nb_levels - 1) / total_range)
                if index == nb_levels - 1:
                    index = nb_levels - 2
                # Convert to color index (into colormap)
                color_index = level_color_list [index]
                # Add to color list for all points
                color_list.append (color_index)

        else:
            # convert from bias to color index in range [6 - 239]

            # Determine color scale
            total_range = bias_max - bias_min
            # Compute bias to color index ratio
            # ratio = (239 - 6) / total_range

            # Determine marker color list: change Id colors by the same token
            color_list = [ int ( round ((bias - bias_min) * (239 - 6) / total_range + 6.))  for bias in bias_list]

        td.Marker.color = color_list

        # Marker ids color: all black
        td.Marker.id_color = [241]*nb_data

        # Offset names from markers symbols
        offset_length = 3. * symbol_size / 50.
        xoffsets = [offset_length]*nb_data
        yoffsets = [-offset_length]*nb_data

        # if two markers are close, change the offset of the second one
        # for i_index in range (1, nb_data):
        #   for j_index in range (i_index):
        #        distance_x = 1. - data[j_index, 0] / data[i_index, 0]
        #        distance_y = 1. - data[j_index, 1] / data[i_index, 1]
        #        if abs (distance_x) < 0.10 and abs (distance_y) < 0.10:
        #            yoffsets [i_index] = - yoffsets [i_index]

        td.Marker.xoffset = xoffsets
        td.Marker.yoffset = yoffsets

    else:
        # Each model/variable is represented on the diagrams as
        # a one-or-two letter code, optionally with a specific color.

        # Impossible to get no marker (bug)
        # instead make all markers white
        td.Marker.symbol = ["dot"]*nb_data
        td.Marker.size = [1]*nb_data
        td.Marker.color = [240]*nb_data

        # Marker ids: color
        td.Marker.id_color = color_list

    # Marker ids: size and font type
    # Adjust the size to diagram scale: typically fct of maximum StdDev
    scale = td.max
    td.Marker.id_size = [int (symbol_size * scale),]*nb_data
    td.Marker.id_font = [1,]*nb_data

    #
    # Plot Taylor diagram
    #

    if skill_score_option:
        # Determine Centered-pattern RMS as skill score function
        # Make a specific data-type off reference std dev value
        reference = reference_std_dev(reference)

        # Compute max_RMS: RMS for point at top left diagram corner
        # where std_dev = scale and correlation = 0.
        max_RMS = reference.compute_RMS_function (scale, 0.)

        # Determine the displayed levels for the skill score
        RMS_levels = vcs.mkscale (0., max_RMS)
        td.skillValues = RMS_levels[1:] # do not include the leading value wich is 0.

        # Plot the diagram with Centered-pattern RMS as skill score
        canvas.plot (data, td, template, skill = reference.compute_RMS_function, bg=background_mode)
    else:
        # Plot the diagram
        canvas.plot (data, td, template, bg=background_mode)

    if bias_list:
        # Create a 2D variable with the bias list (1D list)

        # Why a 2D variable: because, in order to plot the color bar,
        # we use a boxfill method that expects 2D data
        bias_data = cdms.createVariable ([bias_list, bias_list])

        if discrete_option:
            # Plot a discrete color bar
            pass
        else:
            # Plot a continuous color bar
            level_color_list = None
            levels = None

        plot_color_bar (canvas, bias_data, 0.67, 0.705, 0.1, 0.9, level_color_list, levels)


#=============================================================================
#=============================================================================

<Function plot_color_bar>= (<-U)

def plot_color_bar (canvas, data, x_left, x_right, y_top, y_bottom, color_list=None, level_list=None):
    """
    Function plot_color_bar
    -----------------------

    Plot a color bar at the specified position on the page,
    The color bar looks 
      - either as the legend of a "boxfill" plot:
        it is a continuous color bar ranging from blue to green to red to purple,
        along with appropriate scale marks.
      - or as the legend of an "isofill" plot: it is a set of color boxes

    Input:
        canvas:     canvas
        data:       dummy data. used to determine range of scale marks
        x_left:     x coordinate of color_bar left side
        x_right:    x coordinate of color_bar right side
        y_top:      y coordinate of color_bar top side
        y_bottom:   y coordinate of color_bar bottom side
        color_list: list of colors for a discrete ("isofill") color bar
                    if set to None, color bar is continuous 
        level_list: list of levels defining the range for each color

    Global var Output:
        color_bar_templ: the template for possible reuse

    Return nothing
    """

    # We use a trick: a template that allows only the legend to be displayed
    global color_bar

    if color_bar == None:
        # Create a new template based on standard template
        color_bar_templ = canvas.createtemplate('color_bar_templ')

        # Allows for legend display
        color_bar_templ.legend.priority = 1
        color_bar_templ.legend.line = "default"
        color_bar_templ.legend.texttable = "default"
        color_bar_templ.legend.textorientation = "defcenter"

        # Don't draw any of all other features
        color_bar_templ.logicalmask.priority = 0
        color_bar_templ.transformation.priority = 0
        color_bar_templ.file.priority =0
        color_bar_templ.function.priority =0
        color_bar_templ.source.priority = 0
        color_bar_templ.dataname.priority =0
        color_bar_templ.title.priority =0
        color_bar_templ.units.priority =0
        color_bar_templ.crdate.priority=0
        color_bar_templ.crtime.priority=0
        color_bar_templ.comment1.priority=0
        color_bar_templ.comment2.priority=0
        color_bar_templ.comment3.priority=0
        color_bar_templ.comment4.priority=0
        color_bar_templ.zname.priority=0
        color_bar_templ.tname.priority=0
        color_bar_templ.xunits.priority=0
        color_bar_templ.yunits.priority=0
        color_bar_templ.zunits.priority=0
        color_bar_templ.tunits.priority=0
        color_bar_templ.xvalue.priority=0
        color_bar_templ.yvalue.priority=0
        color_bar_templ.zvalue.priority=0
        color_bar_templ.tvalue.priority=0
        color_bar_templ.mean.priority=0
        color_bar_templ.min.priority=0
        color_bar_templ.max.priority=0
        color_bar_templ.xtic1.priority=0
        color_bar_templ.xtic2.priority=0
        color_bar_templ.xmintic1.priority=0
        color_bar_templ.xmintic2.priority=0
        color_bar_templ.ytic1.priority=0
        color_bar_templ.ytic2.priority=0
        color_bar_templ.ymintic1.priority=0
        color_bar_templ.ymintic2.priority=0
        color_bar_templ.xlabel1.priority=0
        color_bar_templ.xlabel2.priority=0
        color_bar_templ.ylabel1.priority=0
        color_bar_templ.ylabel2.priority=0
        color_bar_templ.box1.priority=0
        color_bar_templ.box2.priority=0
        color_bar_templ.box3.priority=0
        color_bar_templ.box4.priority=0
        color_bar_templ.xtic2.priority=0
        color_bar_templ.ytic2.priority=0
        color_bar_templ.data.priority=0

        # Create the color bar legend (namely 'Bias')
        bias_text = canvas.createtext ('bias_text', 'default', 'bias_text', 'default')
        bias_text.viewport = [0, 1, 0, 1]
        bias_text.worldcoordinate = [0, 1, 0, 1]
        bias_text.string =['Bias']
        bias_text.angle = 270
        bias_text.halign = 'center'
        bias_text.valign = 'bottom'
        bias_text.height = 29

        # Memorise the pair of graphic elements
        color_bar = [color_bar_templ, bias_text]
    else:
        color_bar_templ = color_bar[0]
        bias_text       = color_bar[1]

    # Adjust color-bar position and size
    color_bar_templ.legend.x1 = x_left
    color_bar_templ.legend.x2 = x_right
    color_bar_templ.legend.y1 = y_top
    color_bar_templ.legend.y2 = y_bottom
    # Adjust color-bar legend position
    bias_text.x = [x_left - 0.01]
    bias_text.y = [(y_top + y_bottom) / 2]

    # plot the color bar legend
    canvas.plot (bias_text, bg=background_mode)

    if color_list:
        # Create an "isofill" method if not yet existing
        if len(color_bar) > 2:
            isofill_method = color_bar[2]
        else:
            isofill_method = canvas.createisofill ('bias_isofill')
            color_bar.append (isofill_method)
        isofill_method.levels = level_list
        isofill_method.fillareacolors = color_list

        # plot the color bar
        canvas.plot (data, color_bar_templ, isofill_method, bg=background_mode)
    else:
        # plot the color bar
        canvas.plot (data, color_bar_templ, "boxfill", bg=background_mode)


#=============================================================================
#=============================================================================

<Function determine_model_list_from_file>= (<-U)

def determine_model_list_from_file (filename, model_list):
    """
    Function determine_model_list_from_file
    ---------------------------------------

    Determine the list of models to process from input file
    and from optional initial model list passed as a parameter.

    Call:
        model_list = determine_model_list_from_file (filename, model_list)

    Input:
        filename:           Input NetCDf file name
        model_list:         Initial model list (or empty list if none)
    """

    # Open first Comparison Statistics input file
    stats_file = cdms.open (filename)

    # list all variables in file
    variable_list = stats_file.listvariables()

    # if model-list has been given as a parameter
    if model_list:

        # Filter models with statistics available in input file
        processed_model_list = []
        for model_name in model_list:
            # Determine if statistics are available for this model
            Std_var_name = model_name + '_StandardDeviation'
            Corr_var_name =  model_name + '_Correlation'
            if (Std_var_name in variable_list) and (Corr_var_name in variable_list):
                # Stats available --> append to processed model list
                processed_model_list.append (model_name)
    else:
        <Extract model names from input file>

    stats_file.close()
    return processed_model_list

#-----------------------------------------------------------------------------

<Extract model names from input file>= (<-U)

# List all variables with '_StandardDeviation' as a suffix (Keep only their prefix)
Std_dev_list = [var[:-18] for var in variable_list if var.endswith ('_StandardDeviation')]
# List all variables with '_Correlation' as a suffix (Keep only their prefix)
Corr_dev_list = [var[:-12] for var in variable_list if var.endswith ('_Correlation')]

# Processed vars are those present in both lists
processed_model_list = [model for model in Std_dev_list if model in Corr_dev_list]

# Sort the list
processed_model_list.sort()


#=============================================================================
#=============================================================================

<Function read_needed_statistics_from_file>= (<-U)

def read_needed_statistics_from_file (filename, model_list, comp_index_list=None, bias_option=0):
    """
    Function read_needed_statistics_from_file
    -----------------------------------------

    Read and return the needed comparison statistics from input file
    for a model list passed as a parameter.
    Read and return also the variable name and units.

    Call:
        variable_name, units, ref_stats, model_stats, bias = \\
            read_needed_statistics_from_file (filename, model_list, comp_index_list, bias_option)
        or
        variable_name, units, ref_stats, model_stats = \\
            read_needed_statistics_from_file (filename, model_list, comp_index_list, 0)

    Input:
        filename:           Input NetCDf file name
        model_list:         Model list
        comp_index_list:    Index list of components statistics to be read in
                            if None , all components are to be read in
        bias_option:        True if bias to be read in and computed

    Output:
        variable_name:      Variable name (e.g. 'Sea Surface Temperature')
        units:              Variable units (e.g 'degree C')
        ref_stats:          Comparison statistics (standard deviations) for reference data
                            Array or scalar depending of number of components
        model_stats:        Array of comparison statistics (std-dev and correlations) for models
                            Rank of array depends of whether 1 or several components
        bias:               List of bias for all models
    """

    # Flag if one component only
    if comp_index_list and len (comp_index_list) == 1:
        only_one_component = 1
        component_index = comp_index_list[0]
    else:
        only_one_component = 0

    # Open first Comparison Statistics input file
    stats_file = cdms.open (filename)

    # Read in variable name
    if hasattr (stats_file, "obs_data_variable_long_name"):
        variable_name = stats_file.obs_data_variable_long_name
    elif hasattr (stats_file, "model_output_variable_long_name"):
        variable_name = stats_file.model_output_variable_long_name
    elif hasattr (stats_file, "obs_data_variable_name"):
        variable_name = stats_file.obs_data_variable_name
    elif hasattr (stats_file, "model_output_variable_name"):
        variable_name = stats_file.model_output_variable_name
    else:
        # By default, take the base name of the file
        base_filename = os.path.basename (filename)
        base_filename_wo_extension = os.path.splitext (base_filename) [0]
        variable_name = base_filename_wo_extension

    # Read observation data (reference) std deviation
    ReferenceStandardDeviation = stats_file('ReferenceStandardDeviation')
    
    # Read units
    if hasattr (ReferenceStandardDeviation, 'units'):
        units = ReferenceStandardDeviation.units
    else:
        units = ""

    if only_one_component:
        # Take only component std-dev
        ReferenceStandardDeviation = float (ReferenceStandardDeviation[component_index])

    if bias_option:
        # Read in reference global annual average
        referenceAverage = stats_file ('ReferenceAverage')
        # Convert it to a scalar
        referenceAverage = float (referenceAverage)

    # Loop over models: extract stats
    data_list = []
    for model_name in model_list:
        # Read model output correlation and std deviation
        standardDeviation = stats_file (model_name + '_StandardDeviation')
        correlation = stats_file (model_name + '_Correlation')

        # Read units if necessary
        if not units:
            if hasattr (standardDeviation, 'units'):
                units = standardDeviation.units

        if only_one_component:
            # Take only one component
            standardDeviation = float (standardDeviation [component_index])
            correlation = float (correlation [component_index])

            # Add the tupple (std-dev, correlation)
            data_list.append ((standardDeviation, correlation))
        else:
            # Convert from masked variables to arrays
            standardDeviation = standardDeviation.getValue()
            correlation = correlation.getValue()

            # Reshape with a second trailing dimension
            standardDeviation.shape = (len(standardDeviation), 1)
            correlation.shape = (len(correlation), 1)

            # Concatenate along second dimension
            data = Numeric.concatenate ((standardDeviation, correlation), 1)

            data_list.append (data)

    if bias_option:
        bias_list = []
        for model_name in model_list:
            # Read in model output global annual average
            average = stats_file (model_name + '_Average')
            # Convert it to a scalar
            average = float (average)

            # Compute bias
            bias = average - referenceAverage

            bias_list.append (bias)
    stats_file.close()

    # Gather all data in one array
    stats_data = Numeric.array (data_list)

    # Select statistics subset if required
    if comp_index_list and len (comp_index_list) > 1:
        stats_data = Numeric.take (stats_data, component_index_list, axis=1)
        ReferenceStandardDeviation = cdms.MV.take (ReferenceStandardDeviation, component_index_list)

    # Make resulting array a Masked Variable
    stats_data = cdms.asVariable (stats_data)

    if bias_option:
        return variable_name, units, ReferenceStandardDeviation, stats_data, bias_list
    else:
        return variable_name, units, ReferenceStandardDeviation, stats_data


#=============================================================================
#=============================================================================

<Class reference_std_dev>= (<-U)

class reference_std_dev (float):
    """
    This class implements the type: standard deviation of a reference variable.
    It is just a float value with one method that will be used in the computation of
    a test variable RMS.
    """
    < Function compute_RMS_function>

< Function compute_RMS_function>= (<-U)

def compute_RMS_function (self, s, R):
    """
    Function compute_RMS_function
    -----------------------------

    Compute and return the centered-pattern RMS of a test variable
    from its standard-deviation and correlation.

    Input:
        self:   reference-variable standard deviation
        s:      test variable standard deviation
        R:      test-variable correlation with reference variable
    """

    RMS_2 = self*self + s*s - 2*self*s*R
    # Compute square root
    return RMS_2** 0.5


#=============================================================================
#=============================================================================

<Main routine>= (<-U)

#=============================================================================
#=============================================================================

# Initialise a global variable
color_bar = None

<Read command-line parameters>
<Read model parameters>
<Initialise VCS graphic state>
<Determine model list>
<Read needed statistics from file>
<Create output directory if needed>
<Plot a Taylor diagram for total space-time field and each statistical subcomponent>
<Restore previous color palette if necessary>

#-----------------------------------------------------------------------------

<Read command-line parameters>= (<-U)

# Set default parameter values
stats_filename = 'model_vs_data_comp_stats.nc'
stats_filename_list = [stats_filename]
parameter_file = None
model_list = None
output_path = None
output_format = "ps"
output_orientation = "landscape"
add_title = ""
add_units = ""
bias_option = 0
color_option = 0
discrete_option = 0
red_blue_option = 0
symbol_size = 50
multi_components_option = 0
multi_variables_option = 0
# Statistical component selection subset
component_index_list = None
only_one_component = 0
# Default background mode
background_mode = 0


# Parse command line
try:
    options, arguments = getopt.getopt(sys.argv[1:], "hf:o:gpt:u:bcdri:svm:C:k", \
    ("help", "param-file=", "output=", "gif", "portrait", "title=", "units=", "bias", "color", "discrete", "red-blue", \
    "symbol-size=", "superimpose", "mix-variables", "models=", "components=", "background"))
except getopt.GetoptError, err:
    # print help information and exit:
    print err.msg
    usage()
    sys.exit(2)

# Analyse command line options
for option, arg in options:
    if option in ("-h", "--help"):
        usage()
        sys.exit(0)

    elif option in ("-f", "--param-file"):
        # if model-list specified from command line option
        if model_list:
            # Warn the user
            print "Option %s is incompatible with option '--models' or -m." % option
            print "%s option ignored." % option
        else:
            parameter_file = arg

    elif option in ("-o", "--output"):
        output_path = arg

    elif option in ("-g", "--gif"):
        output_format = "gif"

    elif option in ("-p", "--portrait"):
        output_orientation = "portrait"

    elif option in ("-t", "--title"):
        add_title = arg

    elif option in ("-u", "--units"):
        add_units = arg

    elif option in ("-b", "--bias"):
        bias_option = 1

    elif option in ("-c", "--color"):
        color_option = 1

    elif option in ("-d", "--discrete"):
        discrete_option = 1

    elif option in ("-r", "--red-blue"):
        red_blue_option = 1

    elif option in ("-i", "--symbol-size"):
        try:
            symbol_size = int(arg)
            if symbol_size > 100: symbol_size = 100
        except:
            print "Error in parameter (%s) with option %s!" % (arg, option)
            usage()
            sys.exit (3)

    elif option in ("-s", "--superimpose"):
        if multi_variables_option:
            print "Option %s is incompatible with option '--mix-variables' or -v." % option
            print "%s option ignored." % option
        else:
            multi_components_option = 1

    elif option in ("-v", "--mix-variables"):
        if multi_components_option:
            print "Option %s is incompatible with option '--superimpose' or -s." % option
            print "%s option ignored." % option
        else:
            multi_variables_option = 1

    elif option in ("-m", "--models"):
        # if model-list specified from parameter file
        if parameter_file:
            # Warn the user
            print "Option %s is incompatible with option '--param_file' or -f." % option
            print "%s option ignored." % option
        else:
            model_list = arg

    elif option in ("-C", "--components"):
        component_rank_list = arg.split (",")
        try:
            component_index_list = [int(rank_string) - 1 for rank_string in component_rank_list]

            # Redefine the components
            component_short_names = [component_short_names [index] for index in component_index_list]
            component_names = [component_names [index] for index in component_index_list]
            nb_components = len (component_names)

        except:
            print "Error in parameter (%s) with option %s!" % (arg, option)
            usage()
            sys.exit (3)
        else:
            if len (component_index_list) == 1:
                only_one_component = 1

    elif option in ("-k", "--background"):
        background_mode = 1


# Read command line arguments
if len(arguments) == 0:
    pass
elif len(arguments) == 1:
    # Input file name
    stats_filename = arguments[0]
    stats_filename_list = arguments
else:
    if multi_variables_option:
        stats_filename_list = arguments
    else:
        print "Unexpected argument: %s" % arguments[1]
        usage()
        sys.exit(1)

#-----------------------------------------------------------------------------

<Read model parameters>= (<-U)

# if no model list neither parameter file specified on the command line
if not model_list and not parameter_file:
    # if one specific parameter file in current directory
    if 'model_dictionary.py' in os.listdir("."):
        # Select that one
        parameter_file = 'model_dictionary.py'

# Parse parameter file
if parameter_file:
    # Read and parse this parameter file
    # Save module doc-string in case 'model_dictionary.py' has one too
    saved_doc = __doc__
    try:
        execfile (parameter_file)
    except IOError, err:
        # print help information and exit:
        print err.msg
        usage()
        sys.exit(4)
    except Exception, err:
        # print help information and exit:
        print "Error parsing file %s" % parameter_file
        print err.msg
        usage()
        sys.exit(4)

    __doc__ = saved_doc

    # if this file defined a variable named 'model_dictionary'
    if 'model_dictionary' in dir():
        # all_dictionary is a sum of dictionaries
        all_dictionary = model_dictionary.copy()

        # Prepare a list of models and other variables in appropriate order
        model_list = model_dictionary.keys()
        model_list.sort()

        # if this file defined a variable named 'other_dictionary'
        if 'other_dictionary' in dir():
            all_dictionary.update (other_dictionary)

            # Update the list of models and other variables in appropriate order
            other_list = other_dictionary.keys()
            other_list.sort()
            model_list = model_list + other_list

        nb_models = len (model_list)

elif model_list:
    model_list = model_list.split (",")
    # Reassign the total number of models
    nb_models = len (model_list)
    # There is no dictionary of model attributes for this model-list
    all_dictionary = {}

    # if one specific parameter file in current directory
    if 'model_dictionary.py' in os.listdir("."):
        print "Warning: Model list defined on command line."
        print "Parameter file model_dictionary.py is ignored."

#-----------------------------------------------------------------------------

<Initialise VCS graphic state>= (<-U)

canvas = vcs.init()

#if output_orientation == "portrait":
#    DOES NOT WORK: changing canvas size does not change PostScript output size
#    canvas.open()
#    canvas_width = canvas.canvasinfo()['width']
#    canvas_height = canvas.canvasinfo()['height']
#    # Use only a part of the page
#    # Define the plot area: the same height/width ratio as landscape
#    # but as wide as page in portrait orientation is
#    canvas.geometry (canvas_height, int (canvas_height * canvas_height / float (canvas_width) ))
#    # Force change by reopening
#    canvas.close()
#    canvas.open()

# Load red-blue palette if necessary
if bias_option and red_blue_option:

    previous_palette = canvas.getcolormap().name
    # Create new color palette
    canvas.createcolormap('blue_to_red_bias')
    # Set it the default one
    canvas.setcolormap('blue_to_red_bias')

    # Define the new palette color cells
    for i in range (6, 240):
        (r, g, b) = red_blue_map_dictionary[i]
        canvas.setcolorcell (i, r, g, b)


#-----------------------------------------------------------------------------

<Restore previous color palette if necessary>= (<-U)

# Load red-blue palette if necessary
if bias_option and red_blue_option:

    canvas.setcolormap (previous_palette)


#-----------------------------------------------------------------------------

<Determine model list>= (<-U)

if multi_variables_option:
    # Determine model list common in all input files
    processed_model_list = model_list

    for filename in stats_filename_list:
        processed_model_list = determine_model_list_from_file (filename, processed_model_list)
else:
    processed_model_list = determine_model_list_from_file (stats_filename, model_list)

#-----------------------------------------------------------------------------

<Read needed statistics from file>= (<-U)

if multi_variables_option:
    # Component to read in is only 1st component (Total space-time field)
    component_index_list = [0]

    ReferenceStandardDeviation = []
    Taylor_data = []
    variable_name_list = []
    for filename in stats_filename_list:
        # Read statistics for one variable, total space-time component
        variable_name, units, ref_std_dev, model_stats = \
            read_needed_statistics_from_file (filename, processed_model_list, component_index_list, 0)

        # Normalize standard deviations for all models
        model_stats [:, 0] =  model_stats [:, 0] / ref_std_dev
        # Normalize accordingly reference std-dev
        ref_std_dev = 1.

        # Build lists of statistics for all variables:
        ReferenceStandardDeviation.append (ref_std_dev)
        Taylor_data.append (model_stats)

        # Build list of variable names
        variable_name_list.append (variable_name)

    # At this point, "Taylor_data" is a list of Masked Variables
    # convert it to a list of Numeric array
    Taylor_data = [data.getValue() for data in Taylor_data]
    # convert it to a single array
    Taylor_data = Numeric.array (Taylor_data)

    # "ReferenceStandardDeviation" is a list of scalars
    nb_variables = len (variable_name_list)

else:
    if bias_option:
        variable_name, units, ReferenceStandardDeviation, Taylor_data, bias_list = \
            read_needed_statistics_from_file (stats_filename, processed_model_list, component_index_list, 1)
    else:
        variable_name, units, ReferenceStandardDeviation, Taylor_data = \
            read_needed_statistics_from_file (stats_filename, processed_model_list, component_index_list, 0)
        bias_list = []

    # At this point, "Taylor_data" is a Masked Variable
    # "ReferenceStandardDeviation" is a Masked Variable

#-----------------------------------------------------------------------------

<Create output directory if needed>= (<-U)

# if running in mode 1 (default) and more than one component selected
if not multi_variables_option and not multi_components_option and not only_one_component:
    # "output_path" supposed to be a directory path
    if output_path:
        # Add a trailing separator
        if output_path[-1] != "/":
            output_path += "/"
        # Test if output directory exist
        try:
            os.listdir (output_path)
        except OSError:
            # Given directory path does not exist or is not a valid directory
            # Try to create it
            try:
                os.makedirs (output_path)
            except OSError, error:
                print "Could not create directory path: %s" % output_path
                print "[Error no %d] %s" % (error.errno, error.strerror)
                sys.exit (3)


#-----------------------------------------------------------------------------

<Plot a Taylor diagram for total space-time field and each statistical subcomponent>= (<-U)
<Determine colors and symbols>
<Prepare a text legend for all taylor diagrams>
<Plot Taylor diagrams>

#-----------------------------------------------------------------------------

<Determine colors and symbols>= (<-U)
<Colors and symbols of models>
<Colors and symbols of other items>


<Colors and symbols of models>= (<-U)

# Position and style of title lines
title_text = canvas.createtextcombined ('td_title', 'std', 'td_title', 'ASD1')
# if additionnal user-defined title
if add_title:
    # There are two title lines
    title_text.x = [0.5, 0.5]
    title_text.y = [0.95, 0.05]
    title_string  = [add_title, ""] # to be completed later on
else:
    # There is only one automatic title line
    title_text.x = [0.5]
    title_text.y = [0.95]
    title_string  = [""]   # to be defined later on

# if model-attributes have been given as parameters
if all_dictionary:

    # Name of points: one-letter code taken from dictionary
    # List all letter-codes of processed models in same order as processing
    code_letter_list = [all_dictionary [model_name]['code_letter'] for model_name in processed_model_list]

    # Color of points and legend associated to each model: depend on command-line option
    # if multiple components or multiple variables on one diagram
    if multi_components_option or multi_variables_option:
        # Legend is black
        color_list = [241] * len (processed_model_list)
        # Points color will be determined later on
    else:
        if color_option:
            # Convert names into color code taken from a dictionary
            color_list = [all_dictionary [model_name]['color'] for model_name in processed_model_list]
        else:
            # All black
            color_list = [241] * len (processed_model_list)

else:
    # Determine code-letters after model/Variable names
    code_letter_list = []
    for model in processed_model_list:
        if model.upper() == 'MEAN':
            code_letter = 'X'
        elif model.upper() == 'MEDIAN':
            code_letter = 'Y'
        else:
            # Take first letter of name
            code_letter = model[0:1]
            # If letter already used
            if code_letter in code_letter_list:
                # Take first and second letter
                code_letter = model[0:2]
        # Append to code-letter list
        code_letter_list.append (code_letter)

    # Color of points and legend associated to each model: depend on command-line option
    # if multiple components or multiple variables on one diagram
    if multi_components_option or multi_variables_option:
        # Legend is black
        color_list = [241] * len (processed_model_list)
        # Points color will be determined later on
    else:
        if color_option:
            # Choose predefined colors
            # Do not use yellow color: skip color code 245
            standard_color_list = [241, 242, 243, 244, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255]

            color_list = []
            for model_index in range (len (processed_model_list)):
                # Cycle thru 14 color index in standard color list
                color = standard_color_list [model_index % 14]
                color_list.append (color)
        else:
            # All black
            color_list = [241] * len (processed_model_list)

<Colors and symbols of other items>= (<-U)

# Do not use yellow color: skip color code 245
standard_color_list = [241, 242, 243, 244, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255]

# if multiple components on one diagram
if multi_components_option:

    # Determine color for each component
    # Cycle thru 14 color index in standard color list
    component_colors = [standard_color_list [component_index % 14] for component_index in range (nb_components) ]

# else if multiple variables on one diagram
elif multi_variables_option:

    # Determine color for each variable
    # Cycle thru 14 color index in standard color list
    variable_colors = [standard_color_list [variable_index % 14] for variable_index in range (nb_variables) ]

#-----------------------------------------------------------------------------

<Prepare a text legend for all taylor diagrams>= (<-U)

legend_list = []
top_legend_y_offset = 0.9

# if multiple components or multiple variables on one diagram
if multi_components_option or multi_variables_option:

    # Start the legend list with component names or variable names
    if multi_components_option:
        name_list   = component_names
        name_colors = component_colors
        nb_names    = nb_components
    else:
        name_list   = variable_name_list
        name_colors = variable_colors
        nb_names    = nb_variables

    for name_index in range (nb_names):

        name = name_list [name_index]
        unique_name = str (hash (name))
        legend = canvas.createtextcombined (unique_name, 'default', unique_name, 'default')
        legend.color = name_colors [name_index]
        legend.string = name
        legend.x = [0.5]
        legend.y = [top_legend_y_offset - 0.03 * name_index]
        legend.height = 20

        legend_list.append (legend)

    next_legend_y_offset = legend.y[0] - 2 * 0.03
else:
    next_legend_y_offset = top_legend_y_offset

# Complete the legend with model names
for model_index in range (len (processed_model_list)):

    model_name = processed_model_list [model_index]

    model_legend = canvas.createtextcombined ('txt_' + str (hash (model_name)), 'default', 'txt_' + str (hash (model_name)), 'default')
    if bias_option:
        model_legend.color = 241 # black
    else:
        model_legend.color = color_list [model_index]
    model_legend.string = code_letter_list [model_index] + ": " + model_name
    model_legend.x = [0.8]
    model_legend.y = [next_legend_y_offset - 0.03 * model_index]
    model_legend.height = 20

    legend_list.append (model_legend)


# Temporary: Prepare a text element for units on the Standard Deviation axis
# Force 'units' given on the command line
if add_units: units = add_units
if units:
    units_legend = canvas.createtextcombined ('txt_' + str (hash (units)), 'default', 'txt_' + str (hash (units)),'default')
    units_legend.string = ['(' + units + ')']
    units_legend.angle = 270
    units_legend.halign = 'left'
    units_legend.valign = 'top'

    units_legend.x = [0.0324]
    units_legend.y = [0.68]
    units_legend.height = 24

#-----------------------------------------------------------------------------

<Plot Taylor diagrams>= (<-U)

# if multiple components on one diagram
if multi_components_option:

    # Determine maximum standard deviation in order to determine diagram scale:
    # "total space-time field" has max standard deviation over all components
    max_std_dev = max (ReferenceStandardDeviation[0], max (Taylor_data [:, 0, 0]))
    scale = max_std_dev * 1.2

    # Add 'units' information
    Taylor_data.units = units

    # Loop over all components
    for component_index in range (nb_components):
        component_name = component_names [component_index]

        reference = ReferenceStandardDeviation [component_index]
        # Convert to a scalar
        reference = float (reference)
        data = Taylor_data [:, component_index, :]

        # Determine color of data points and reference circle
        color = component_colors [component_index]
        color_list = [color] * len (processed_model_list)

        # Skill score isolines option: set to off
        skill_score_option = 0

        plot_taylor_diagram (reference, data, code_letter_list, color_list, [], 0, 0, skill_score_option, scale, 1, color, component_name)

    # Add a title
    title_string [-1] = variable_name + " Taylor diagram"
    title_text.string = title_string
    canvas.plot (title_text, bg=background_mode)

    # Add the legend
    for legend in legend_list:
        canvas.plot (legend, bg=background_mode)

    if not background_mode:
        print "This program will save the image in a file."
        print "Press <Enter> to continue..."
        sys.stdin.readline()

    # Output an image file: by default, named after variable name
    if not output_path:
        if output_format == "ps":
            output_path = variable_name + ".ps"
        else:
            output_path = variable_name + "gif"

    if output_orientation == "portrait":
        if output_format == "ps":
            canvas.postscript (output_path, 'p')
        else:
            canvas.gif (output_path, 'r', 'p')
    else:
        if output_format == "ps":
            canvas.postscript (output_path, 'l')
        else:
            canvas.gif (output_path, 'r', 'l')

# else if multiple variables on one diagram
elif multi_variables_option:

    # Determine maximum standard deviation for all models and for all variables
    # in order to determine diagram scale:
    max_std_dev = max (ReferenceStandardDeviation[0], max (Numeric.ravel (Taylor_data [:, :, 0])))
    scale = max_std_dev * 1.2

    # Loop over all variables
    for variable_index in range (nb_variables):
        variable_name = variable_name_list [variable_index]

        reference = ReferenceStandardDeviation [variable_index]     # Namely 1.0 since normalized
        data = Taylor_data [variable_index, :, :]

        # Determine color of data points and reference circle
        color = variable_colors [variable_index]
        color_list = [color] * len (processed_model_list)
        # Color of reference circle = black
        ref_color = 241

        # Skill score isolines option: set to on
        skill_score_option = 1

        plot_taylor_diagram (reference, data, code_letter_list, color_list, [], 0, 0, skill_score_option, scale, 1, ref_color, variable_name)

    # Add a title
    title_string [-1] = "Multi-variable Taylor diagram"
    title_text.string = title_string
    canvas.plot (title_text, bg=background_mode)

    # Add the legend
    for legend in legend_list:
        canvas.plot (legend, bg=background_mode)

    if not background_mode:
        print "This program will save the image in a file."
        print "Press <Enter> to continue..."
        sys.stdin.readline()

    # Output an image file: by default, named after first variable name
    if not output_path:
        if output_format == "ps":
            output_path = variable_name_list[0] + ".ps"
        else:
            output_path = variable_name_list[0] + "gif"

    if output_orientation == "portrait":
        if output_format == "ps":
            canvas.postscript (output_path, 'p')
        else:
            canvas.gif (output_path, 'r', 'p')
    else:
        if output_format == "ps":
            canvas.postscript (output_path, 'l')
        else:
            canvas.gif (output_path, 'r', 'l')
else:

    # Add 'units' information
    Taylor_data.units = units

    # Loop over all components

    for component_index in range (nb_components):
        component_name = component_names [component_index]

        if only_one_component:
            # Data read in from file already with 'component' dimension removed
            reference = ReferenceStandardDeviation
            data = Taylor_data
        else:
            reference = ReferenceStandardDeviation [component_index]
            # Convert to a scalar
            reference = float (reference)
            data = Taylor_data [:, component_index, :]

        # Determine if correlations have negative values
        min_correl = min (data [:, 1])
        # Determine if Taylor diagram with two quadrans is necessary
        if min_correl < -2e-2:
            quadrans = 2
        else:
            quadrans = 1
        
        # Two quadrans: not supported yet
        if quadrans == 2:
            print "Warning: because of negative correlation,"
            print "points may be out of drawing area for component '%s'..." % component_name
            quadrans = 1
            
            # Discard points out of draw area
            keep_condition = (data.getValue() [:, 1] >= -2e-2)
            # cdms.MV.compress does not work
            # data = cdms.MV.compress (keep_condition, data, dimension=1)
            indices = Numeric.arange (len (data))
            kept_indices = Numeric.compress (keep_condition, indices)
            data = cdms.MV.take (data, kept_indices)


        # Reference circle color: black
        reference_color = 241

        # Skill score isolines option: set to on
        skill_score_option = 1

        # clear canvas
        canvas.clear()

        if component_index == 0:
            plot_taylor_diagram (reference, data, code_letter_list, color_list, bias_list, discrete_option, red_blue_option, \
                skill_score_option, None, quadrans, reference_color, component_name)
        else:
            plot_taylor_diagram (reference, data, code_letter_list, color_list, [], 0, 0, \
                skill_score_option, None, quadrans, reference_color, component_name)

        # Add a title
        title_string [-1] = component_name
        title_text.string = title_string
        canvas.plot (title_text, bg=background_mode)

        # Add a legend
        for legend in legend_list:
            canvas.plot (legend, bg=background_mode)

        if not background_mode:
            print "This program will save the image in a file."
            print "Press <Enter> to continue..."
            sys.stdin.readline()

        # Output an image file
        component_short_name = component_short_names [component_index]
        image_file_name = string.replace (component_short_name, " ", "_")
        postscript_file_name = image_file_name  + ".ps"
        gif_file_name = image_file_name + ".gif"

        # if output specified on command line
        if output_path:
            if only_one_component:
                # Assume 'output_path' is a file name
                postscript_file_name = output_path
                gif_file_name = output_path
            else:
                postscript_file_name = os.path.join (output_path, postscript_file_name)
                gif_file_name = os.path.join (output_path, gif_file_name)

        if output_orientation == "portrait":
            if output_format == "ps":
                canvas.postscript (postscript_file_name, 'p')
            else:
                canvas.gif (gif_file_name, 'r', 'p')
        else:
            if output_format == "ps":
                canvas.postscript (postscript_file_name, 'l')
            else:
                canvas.gif (gif_file_name, 'r', 'l')

# print
# print "End...Press <Enter> to continue..."
# sys.stdin.readline()