Commit cdba3f1d authored by Florentie, Liesbeth's avatar Florentie, Liesbeth
Browse files

update to included latitude and longitude in diagnostics outputfile optimizer.nc

parent d8ba3e4c
......@@ -351,7 +351,7 @@ class TotalColumnObservations(Observations):
# write data required by observation operator for sampling to file
f = io.CT_CDF(obsinputfile, method='create')
logging.debug('Creating new observations file for ObservationOperator (%s)' % obsinputfile)
logging.debug('Creating new observations file for ObservationOperator (%s) containing %d observations' % (obsinputfile,len(self.datalist)))
dimsoundings = f.add_dim('soundings', len(self.datalist))
dimlevels = f.add_dim('levels', self.getvalues('pressure').shape[1])
......
......@@ -210,7 +210,7 @@ class ObsPackObservations(Observations):
logging.debug("No observations found for this time period, nothing written to obs file")
else:
f = io.CT_CDF(obsinputfile, method='create')
logging.debug('Creating new observations file for ObservationOperator (%s)' % obsinputfile)
logging.debug('Creating new observations file for ObservationOperator (%s) containing %d observations' % (obsinputfile,len(self.datalist)))
dimid = f.add_dim('obs', len(self.datalist))
dim200char = f.add_dim('string_of200chars', 200)
......
"""CarbonTracker Data Assimilation Shell (CTDAS) Copyright (C) 2017 Wouter Peters.
"""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.
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.
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
You should have received a copy of the GNU General Public License along with this
program. If not, see <http://www.gnu.org/licenses/>."""
#!/usr/bin/env python
# optimizer.py
"""
Author : peters
Author : peters
Revision History:
File created on 28 Jul 2010.
......@@ -24,6 +24,8 @@ File created on 28 Jul 2010.
import os
import sys
import logging
import numpy as np
import da.tools.io4 as io
sys.path.append(os.getcwd())
from da.baseclasses.optimizer import Optimizer
......@@ -41,6 +43,123 @@ class CO2Optimizer(Optimizer):
All other methods are inherited from the base class Optimizer.
"""
def create_matrices(self):
""" Create Matrix space needed in optimization routine """
# mean state [X]
self.x = np.zeros((self.nlag * self.nparams,), float)
# deviations from mean state [X']
self.X_prime = np.zeros((self.nlag * self.nparams, self.nmembers,), float)
# mean state, transported to observation space [ H(X) ]
self.Hx = np.zeros((self.nobs,), float)
# deviations from mean state, transported to observation space [ H(X') ]
self.HX_prime = np.zeros((self.nobs, self.nmembers), float)
# observations
self.obs = np.zeros((self.nobs,), float)
# observation ids
self.obs_ids = np.zeros((self.nobs,), float)
# covariance of observations
# Total covariance of fluxes and obs in units of obs [H P H^t + R]
if self.algorithm == 'Serial':
self.R = np.zeros((self.nobs,), float)
self.HPHR = np.zeros((self.nobs,), float)
else:
self.R = np.zeros((self.nobs, self.nobs,), float)
self.HPHR = np.zeros((self.nobs, self.nobs,), float)
# localization of obs
self.may_localize = np.zeros(self.nobs, bool)
# rejection of obs
self.may_reject = np.zeros(self.nobs, bool)
# flags of obs
self.flags = np.zeros(self.nobs, int)
# species type
self.species = np.zeros(self.nobs, str)
# species type
self.sitecode = np.zeros(self.nobs, str)
# rejection_threshold
self.rejection_threshold = np.zeros(self.nobs, float)
# lat and lon
self.latitude = np.zeros(self.nobs, float)
self.longitude = np.zeros(self.nobs, float)
# species mask
self.speciesmask = {}
# Kalman Gain matrix
#self.KG = np.zeros((self.nlag * self.nparams, self.nobs,), float)
self.KG = np.zeros((self.nlag * self.nparams,), float)
def state_to_matrix(self, statevector):
allsites = [] # collect all obs for n=1,..,nlag
alllats = [] # collect all latitudes for n=1,..,nlag
alllons = [] # collect all longitudes for n=1,..,nlag
allobs = [] # collect all obs for n=1,..,nlag
allmdm = [] # collect all mdm for n=1,..,nlag
allids = [] # collect all model samples for n=1,..,nlag
allreject = [] # collect all model samples for n=1,..,nlag
alllocalize = [] # collect all model samples for n=1,..,nlag
allflags = [] # collect all model samples for n=1,..,nlag
allspecies = [] # collect all model samples for n=1,..,nlag
allsimulated = [] # collect all members model samples for n=1,..,nlag
allrej_thres = [] # collect all rejection_thresholds, will be the same for all samples of same source
for n in range(self.nlag):
samples = statevector.obs_to_assimilate[n]
members = statevector.ensemble_members[n]
self.x[n * self.nparams:(n + 1) * self.nparams] = members[0].param_values
self.X_prime[n * self.nparams:(n + 1) * self.nparams, :] = np.transpose(np.array([m.param_values for m in members]))
# Add observation data for all sample objects
if samples != None:
if type(samples) != list: samples = [samples]
for m in range(len(samples)):
sample = samples[m]
logging.debug('Lag %i, sample %i: rejection_threshold = %i, nobs = %i' %(n, m, sample.rejection_threshold, sample.getlength()))
allrej_thres.extend([sample.rejection_threshold] * sample.getlength())
allreject.extend(sample.getvalues('may_reject'))
alllocalize.extend(sample.getvalues('may_localize'))
allflags.extend(sample.getvalues('flag'))
allspecies.extend(sample.getvalues('species'))
allobs.extend(sample.getvalues('obs'))
allsites.extend(sample.getvalues('code'))
alllats.extend(sample.getvalues('lat'))
alllons.extend(sample.getvalues('lon'))
allmdm.extend(sample.getvalues('mdm'))
allids.extend(sample.getvalues('id'))
simulatedensemble = sample.getvalues('simulated')
for s in range(simulatedensemble.shape[0]):
allsimulated.append(simulatedensemble[s])
self.rejection_threshold[:] = np.array(allrej_thres)
self.obs[:] = np.array(allobs)
self.obs_ids[:] = np.array(allids)
self.HX_prime[:, :] = np.array(allsimulated)
self.Hx[:] = self.HX_prime[:, 0]
self.may_reject[:] = np.array(allreject)
self.may_localize[:] = np.array(alllocalize)
self.flags[:] = np.array(allflags)
self.species[:] = np.array(allspecies)
self.sitecode = allsites
self.latitude[:] = np.array(alllats)
self.longitude[:] = np.array(alllons)
self.X_prime = self.X_prime - self.x[:, np.newaxis] # make into a deviation matrix
self.HX_prime = self.HX_prime - self.Hx[:, np.newaxis] # make a deviation matrix
if self.algorithm == 'Serial':
for i, mdm in enumerate(allmdm):
self.R[i] = mdm ** 2
else:
for i, mdm in enumerate(allmdm):
self.R[i, i] = mdm ** 2
def set_localization(self, loctype='None'):
""" determine which localization to use """
......@@ -55,12 +174,12 @@ class CO2Optimizer(Optimizer):
elif self.nmembers == 150:
self.tvalue = 1.97591
elif self.nmembers == 200:
self.tvalue = 1.9719
else: self.tvalue = 0
self.tvalue = 1.9719
else: self.tvalue = 0
else:
self.localization = False
self.localizetype = 'None'
logging.info("Current localization option is set to %s" % self.localizetype)
if self.localization == True:
if self.tvalue == 0:
......@@ -72,9 +191,9 @@ class CO2Optimizer(Optimizer):
""" localize the Kalman Gain matrix """
import numpy as np
if not self.localization:
if not self.localization:
logging.debug('Not localized observation %i' % self.obs_ids[n])
return
return
if self.localizetype == 'CT2007':
count_localized = 0
for r in range(self.nlag * self.nparams):
......@@ -92,9 +211,174 @@ class CO2Optimizer(Optimizer):
self.algorithm = 'Serial'
else:
self.algorithm = 'Bulk'
logging.info("Current minimum least squares algorithm is set to %s" % self.algorithm)
def write_diagnostics(self, filename, type):
"""
Open a NetCDF file and write diagnostic output from optimization process:
- calculated residuals
- model-data mismatches
- HPH^T
- prior ensemble of samples
- posterior ensemble of samples
- prior ensemble of fluxes
- posterior ensemble of fluxes
The type designation refers to the writing of prior or posterior data and is used in naming the variables"
"""
# Open or create file
if type == 'prior':
f = io.CT_CDF(filename, method='create')
logging.debug('Creating new diagnostics file for optimizer (%s)' % filename)
elif type == 'optimized':
f = io.CT_CDF(filename, method='write')
logging.debug('Opening existing diagnostics file for optimizer (%s)' % filename)
# Add dimensions
dimparams = f.add_params_dim(self.nparams)
dimmembers = f.add_members_dim(self.nmembers)
dimlag = f.add_lag_dim(self.nlag, unlimited=False)
dimobs = f.add_obs_dim(self.nobs)
dimstate = f.add_dim('nstate', self.nparams * self.nlag)
dim200char = f.add_dim('string_of200chars', 200)
# Add data, first the ones that are written both before and after the optimization
savedict = io.std_savedict.copy()
savedict['name'] = "statevectormean_%s" % type
savedict['long_name'] = "full_statevector_mean_%s" % type
savedict['units'] = "unitless"
savedict['dims'] = dimstate
savedict['values'] = self.x.tolist()
savedict['comment'] = 'Full %s state vector mean ' % type
f.add_data(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "statevectordeviations_%s" % type
savedict['long_name'] = "full_statevector_deviations_%s" % type
savedict['units'] = "unitless"
savedict['dims'] = dimstate + dimmembers
savedict['values'] = self.X_prime.tolist()
savedict['comment'] = 'Full state vector %s deviations as resulting from the optimizer' % type
f.add_data(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "modelsamplesmean_%s" % type
savedict['long_name'] = "modelsamplesforecastmean_%s" % type
savedict['units'] = "mol mol-1"
savedict['dims'] = dimobs
savedict['values'] = self.Hx.tolist()
savedict['comment'] = '%s mean mole fractions based on %s state vector' % (type, type)
f.add_data(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "modelsamplesdeviations_%s" % type
savedict['long_name'] = "modelsamplesforecastdeviations_%s" % type
savedict['units'] = "mol mol-1"
savedict['dims'] = dimobs + dimmembers
savedict['values'] = self.HX_prime.tolist()
savedict['comment'] = '%s mole fraction deviations based on %s state vector' % (type, type)
f.add_data(savedict)
# Continue with prior only data
if type == 'prior':
savedict = io.std_savedict.copy()
savedict['name'] = "sitecode"
savedict['long_name'] = "site code propagated from observation file"
savedict['dtype'] = "char"
savedict['dims'] = dimobs + dim200char
savedict['values'] = self.sitecode
savedict['missing_value'] = '!'
f.add_data(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = 'latitude'
savedict['dims'] = dimobs
savedict['units'] = "degrees_north"
savedict['missing_value'] = -999.9
savedict['values'] = self.latitude
f.add_data(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = 'longitude'
savedict['units'] = "degrees_east"
savedict['missing_value'] = -999.9
savedict['dims'] = dimobs
savedict['values'] = self.longitude
f.add_data(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "observed"
savedict['long_name'] = "observedvalues"
savedict['units'] = "mol mol-1"
savedict['dims'] = dimobs
savedict['values'] = self.obs.tolist()
savedict['comment'] = 'Observations used in optimization'
f.add_data(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "obs_num"
savedict['dtype'] = "int64"
savedict['long_name'] = "Unique_observation_number"
savedict['units'] = ""
savedict['dims'] = dimobs
savedict['values'] = self.obs_ids.tolist()
savedict['comment'] = 'Unique observation number across the entire ObsPack or satellite retrieval distribution'
f.add_data(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "modeldatamismatchvariance"
savedict['long_name'] = "modeldatamismatch variance"
savedict['units'] = "[mol mol-1]^2"
if self.algorithm == 'Serial':
savedict['dims'] = dimobs
else: savedict['dims'] = dimobs + dimobs
savedict['values'] = self.R.tolist()
savedict['comment'] = 'Variance of mole fractions resulting from model-data mismatch'
f.add_data(savedict)
# Continue with posterior only data
elif type == 'optimized':
savedict = io.std_savedict.copy()
savedict['name'] = "totalmolefractionvariance"
savedict['long_name'] = "totalmolefractionvariance"
savedict['units'] = "[mol mol-1]^2"
if self.algorithm == 'Serial':
savedict['dims'] = dimobs
else: savedict['dims'] = dimobs + dimobs
savedict['values'] = self.HPHR.tolist()
savedict['comment'] = 'Variance of mole fractions resulting from prior state and model-data mismatch'
f.add_data(savedict)
savedict = io.std_savedict.copy()
savedict['name'] = "flag"
savedict['long_name'] = "flag_for_obs_model"
savedict['units'] = "None"
savedict['dims'] = dimobs
savedict['values'] = self.flags.tolist()
savedict['comment'] = 'Flag (0/1/2/99) for observation value, 0 means okay, 1 means QC error, 2 means rejected, 99 means not sampled'
f.add_data(savedict)
#savedict = io.std_savedict.copy()
#savedict['name'] = "kalmangainmatrix"
#savedict['long_name'] = "kalmangainmatrix"
#savedict['units'] = "unitless molefraction-1"
#savedict['dims'] = dimstate + dimobs
#savedict['values'] = self.KG.tolist()
#savedict['comment'] = 'Kalman gain matrix of all obs and state vector elements'
#dummy = f.add_data(savedict)
f.close()
logging.debug('Diagnostics file closed')
################### End Class CO2Optimizer ###################
if __name__ == "__main__":
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment