observationoperator.py 6.56 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
"""CarbonTracker Data Assimilation Shell (CTDAS) Copyright (C) 2017 Wouter Peters. 
Users are recommended to contact the developers (wouter.peters@wur.nl) to receive
updates of the code. See also: http://www.carbontracker.eu. 

This program is free software: you can redistribute it and/or modify it under the
terms of the GNU General Public License as published by the Free Software Foundation, 
version 3. This program is distributed in the hope that it will be useful, but 
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 
FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. 

You should have received a copy of the GNU General Public License along with this 
program. If not, see <http://www.gnu.org/licenses/>."""
13
14
15
16
#!/usr/bin/env python
# model.py

"""
Peters, Wouter's avatar
Peters, Wouter committed
17
18
.. module:: observationoperator
.. moduleauthor:: Wouter Peters 
19
20
21
22
23
24
25
26

Revision History:
File created on 30 Aug 2010.

"""

import logging

27
28
identifier = 'RandomizerObservationOperator'
version = '1.0'
29
30
31
32
33
34
35
36
37
38

################### Begin Class ObservationOperator ###################
class ObservationOperator(object):
    """
    Testing
    =======
    This is a class that defines an ObervationOperator. This object is used to control the sampling of
    a statevector in the ensemble Kalman filter framework. The methods of this class specify which (external) code
    is called to perform the sampling, and which files should be read for input and are written for output.

39
    The baseclasses consist mainly of empty methods that require an application specific application. The baseclass will take observed values, and perturb them with a random number chosen from the model-data mismatch distribution. This means no real operator will be at work, but random normally distributed residuals will come out of y-H(x) and thus the inverse model can proceed. This is mainly for testing the code...
40
41
42

    """

43
    def __init__(self, dacycle=None):
44
        """ The instance of an ObservationOperator is application dependent """
45
46
47
        self.ID = identifier
        self.version = version
        self.restart_filelist = []
48
        self.output_filelist = []
49
        self.outputdir = None # Needed for opening the samples.nc files created 
50

51
        logging.info('Observation Operator object initialized: %s' % self.ID)
52

53
        # The following code allows the object to be initialized with a dacycle object already present. Otherwise, it can
54
55
        # be added at a later moment.

56
57
        if dacycle != None:
            self.dacycle = dacycle
58
        else:
59
            self.dacycle = {}
60

61
    
62
    def get_initial_data(self):
63
        """ This method places all initial data needed by an ObservationOperator in the proper folder for the model """
64

65
    def setup(self,dacycle):
66
67
        """ Perform all steps necessary to start the observation operator through a simple Run() call """

68
        self.dacycle = dacycle
brunner's avatar
brunner committed
69
        self.outputdir = dacycle['dir.output']
70
71
72
73

    def prepare_run(self):
        """ Prepare the running of the actual forecast model, for example compile code """

brunner's avatar
brunner committed
74
        import os
75

brunner's avatar
brunner committed
76
        # Define the name of the file that will contain the modeled output of each observation
77

brunner's avatar
brunner committed
78
79
        self.simulated_file = os.path.join(self.outputdir, 'samples_simulated.%s.nc' % self.dacycle['time.sample.stamp'])
        self.forecast_nmembers = int(self.dacycle['da.optimizer.nmembers'])
80

81
    def validate_input(self):
82
83
84
        """ Make sure that data needed for the ObservationOperator (such as observation input lists, or parameter files)
            are present.
        """
85
    def save_data(self):
86
87
        """ Write the data that is needed for a restart or recovery of the Observation Operator to the save directory """

88
    def run(self):
brunner's avatar
brunner committed
89
90
91
92
        """
         This Randomizer will take the original observation data in the Obs object, and simply copy each mean value. Next, the mean 
         value will be perturbed by a random normal number drawn from a specified uncertainty of +/- 2 ppm
        """
93

brunner's avatar
brunner committed
94
95
        import da.tools.io4 as io
        import numpy as np
96

97
        # Create a flask output file to hold simulated values for later reading
98

brunner's avatar
brunner committed
99
100
101
        f = io.CT_CDF(self.simulated_file, method='create')
        logging.debug('Creating new simulated observation file in ObservationOperator (%s)' % self.simulated_file)
        
102
        dimid = f.createDimension('obs_num', size=None)
brunner's avatar
brunner committed
103
        dimid = ('obs_num',)
104
105
106
107
108
109
110
111
112
113
        savedict = io.std_savedict.copy() 
        savedict['name'] = "obs_num"
        savedict['dtype'] = "int"
        savedict['long_name'] = "Unique_Dataset_observation_index_number"
        savedict['units'] = ""
        savedict['dims'] = dimid
        savedict['comment'] = "Unique index number within this dataset ranging from 0 to UNLIMITED."
        f.add_data(savedict,nsets=0)

        dimmember = f.createDimension('nmembers', size=self.forecast_nmembers)
brunner's avatar
brunner committed
114
        dimmember = ('nmembers',)
115
116
117
118
119
120
121
122
123
        savedict = io.std_savedict.copy() 
        savedict['name'] = "flask"
        savedict['dtype'] = "float"
        savedict['long_name'] = "mole_fraction_of_trace_gas_in_air"
        savedict['units'] = "mol tracer (mol air)^-1"
        savedict['dims'] = dimid + dimmember
        savedict['comment'] = "Simulated model value created by RandomizerObservationOperator"
        f.add_data(savedict,nsets=0)

brunner's avatar
brunner committed
124
        # Open file with x,y,z,t of model samples that need to be sampled
125
126
127

        f_in = io.ct_read(self.dacycle['ObsOperator.inputfile'],method='read') 

brunner's avatar
brunner committed
128
        # Get simulated values and ID
129
130
131
132
133

        ids = f_in.get_variable('obs_num')
        obs = f_in.get_variable('observed')
        mdm = f_in.get_variable('modeldatamismatch')

brunner's avatar
brunner committed
134
        # Loop over observations, add random white noise, and write to file
135

brunner's avatar
brunner committed
136
137
138
        for i,data in enumerate(zip(ids,obs,mdm)):
            f.variables['obs_num'][i] = data[0]
            f.variables['flask'][i,:] = data[1]+np.random.randn(self.forecast_nmembers)*data[2]
139

brunner's avatar
brunner committed
140
141
        f.close()
        f_in.close()
142

brunner's avatar
brunner committed
143
        # Report success and exit
144

brunner's avatar
brunner committed
145
        logging.info('ObservationOperator finished successfully, output file written (%s)' % self.simulated_file)
146
147
148
149
150
151

    def run_forecast_model(self):
        self.prepare_run()
        self.validate_input()
        self.run()
        self.save_data()
152
153
154
155


################### End Class ObservationOperator ###################

156
157
158
159
160
161
class RandomizerObservationOperator(ObservationOperator):
    """ This class holds methods and variables that are needed to use a random number generated as substitute
        for a true observation operator. It takes observations and returns values for each obs, with a specified 
        amount of white noise added 
    """

162
163
164
165


if __name__ == "__main__":
    pass