Commit 945546cf authored by brunner's avatar brunner
Browse files

No commit message

No commit message
parent 352ba28b
This diff is collapsed.
#!/usr/bin/env python
# model.py
import logging
identifier = 'RandomizerObservationOperator'
version = '1.0'
################### 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.
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...
"""
def __init__(self, dacycle=None):
""" The instance of an ObservationOperator is application dependent """
self.ID = identifier
self.version = version
self.restart_filelist = []
self.output_filelist = []
self.outputdir = None # Needed for opening the samples.nc files created
logging.info('Observation Operator object initialized: %s' % self.ID)
# The following code allows the object to be initialized with a dacycle object already present. Otherwise, it can
# be added at a later moment.
if dacycle != None:
self.dacycle = dacycle
else:
self.dacycle = {}
def get_initial_data(self):
""" This method places all initial data needed by an ObservationOperator in the proper folder for the model """
def setup(self,dacycle):
""" Perform all steps necessary to start the observation operator through a simple Run() call """
self.dacycle = dacycle
self.outputdir = dacycle['dir.output']
def prepare_run(self):
""" Prepare the running of the actual forecast model, for example compile code """
import os
# Define the name of the file that will contain the modeled output of each observation
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'])
def run(self):
import da.tools.io4 as io
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import os
from netCDF4 import Dataset
import sys
from datetime import datetime, timedelta
import subprocess
# from cdo import *
from dateutil import rrule
from datetime import datetime, timedelta
f = io.CT_CDF(self.simulated_file, method='create')
logging.debug('Creating new simulated observation file in ObservationOperator (%s)' % self.simulated_file)
dimid = f.createDimension('obs_num', size=None)
dimid = ('obs_num',)
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)
dimmember = ('nmembers',)
savedict = io.std_savedict.copy()
savedict['name'] = "flask"
savedict['dtype'] = "float"
savedict['long_name'] = "mole_fraction_of_trace_gas_in_air"
savedict['units'] = "mol mol-1"
savedict['dims'] = dimid + dimmember
savedict['comment'] = "Simulated model value created by COSMO"
f.add_data(savedict,nsets=0)
# Open file with x,y,z,t of model samples that need to be sampled
f_in = io.ct_read(self.dacycle['ObsOperator.inputfile'],method='read')
# Get simulated values and ID
ids = f_in.get_variable('obs_num')
obs = f_in.get_variable('observed')
mdm = f_in.get_variable('modeldatamismatch')
f_in.close()
# Loop over observations, add random white noise, and write to file
shape = (self.forecast_nmembers,mdm.size)
model_data=np.empty(shape=shape) # 3x7
# sys.exit()
# self.obspack_dir = dacycle.dasystem['obspack.input.dir']
# infile = os.path.join(self.obspack_dir, 'summary', '%s_dataset_summary.txt' % (self.obspack_id,))
# infile = "/store/empa/em05/parsenov/obspack/summary/obspack_co2_1_GLOBALVIEWplus_v3.2_2017-11-02_dataset_summary.txt"
# f = open(infile, 'r')
# lines = f.readlines()
# f.close()
# ncfilelist = []
# for line in lines:
# if not line.startswith('# dataset:'): continue
# items = line.split(':')
# ncfile = items[1].strip()
# ncfilelist += [ncfile]
# for ncfile in ncfilelist:
# infile = os.path.join(ncfile + '.nc')
# here comes int2lm
# cosmo = os.path.join(dacycle.dasystem['cosmo_path'],"run_chain.py")
# os.system('python '+cosmo+'ctdas'+dacycle['time.start'].strftime('%Y-%m-%d')+'0 672 -j meteo icbc emissions biofluxes int2lm post_int2lm')
# HERE MULTIPLY COSMO FLUXES INT2LM WITH CTDAS PARAMS
for ens in range(1,self.forecast_nmembers+1):
ens = str(ens).zfill(3)
os.system('ln *nc /scratch/snx3000/parsenov/ctdas/2013040100_0_168/int2lm/output /scratch/snx3000/parsenov/ctdas/cosmo_ensemble/'+ens)
os.system('rm /scratch/snx3000/parsenov/ctdas/cosmo_ensemble/*f.nc')
for dt in rrule.rrule(rrule.HOURLY, dtstart=self['time.start'], until=self['time.start']+timedelta(hours=168)):
print(dt)
cdo.mul("parameters."+ens+".nc", input = "/scratch/snx3000/parsenov/ctdas/2013040100_0_168/int2lm/output/laf"+dt.strftime('%Y%m%d%H')+"f.nc", output = "/scratch/snx3000/parsenov/ctdas/cosmo_ensemble/"+ens+"laf"+dt.strftime('%Y%m%d%H')+"f.nc")
# here comes COSMO
sys.exit()
os.system('python '+cosmo+' ctdas '+self['time.start'].strftime('%Y-%m-%d')+' 0 672 -j cosmo')
extract_model_data(self.forecast_nmembers)
for i in range(0,self.forecast_nmembers):
idx=str(i+1)
cosmo_file = os.path.join('/scratch/snx3000/parsenov/processing_chain/ctdas_'+idx+'/%s/cosmo/output/model_'+idx+'_%s.nc' % self['time.sample.stamp'])
ifile = Dataset(cosmo_file, mode='r')
model_data[i,:] = np.squeeze(ifile.variables['CO2'][:])#*29./44.)#*1000000. # in ppm
ifile.close()
for j,data in enumerate(zip(ids,obs,mdm)):
f.variables['obs_num'][j] = data[0]
f.variables['flask'][j,:] = model_data[:,j]
# print model_data[:,j]
f.close()
# f.variables['flask'][i,:] = data[1]+np.random.randn(self.forecast_nmembers)
#### WARNING ACHTUNG PAZNJA POZOR VNEMANIE data[2] is model data mismatch (=1000) by default in tools/io4.py!!! pavle
# print i # i is num of hour
# f.variables['flask'][i,:] = data[1]+np.random.randn(self.forecast_nmembers)*data[2]
# print 'pavle',data[1]
# print 'pavle',np.random.randn(self.forecast_nmembers)
# print 'pavle',data[2]
# print 'i',i
# print 'data',data
# print 'model data',model_data
# print 'ids',ids # number of obs [11841779 11841780 and such] dim =7
# print 'obs',obs # observation data [0.00040033 0.0004001 and suchh] dim =7
# print 'mdm',mdm # model data mismatch (=1000.) dim = 7
# print self.forecast_nmembers
# sys.exit()
# Report success and exit
logging.info('ObservationOperator finished successfully, output file written (%s)' % self.simulated_file)
def run_forecast_model(self):
self.prepare_run()
self.run()
def extract_model_data(ens):
import os
import sys
from datetime import datetime, timedelta
import subprocess
from . import site_height
#from cdo import *
cdo = Cdo()
files2cat=[]
time_stamp = self['time.sample.stamp']
cosmo_out = self['da.cosmo_path']+'/'+ens+'/'+time_stamp+'cosmo/output'
os.chdir(cosmo_out)
for time in tools.iter_hours(starttime, hstart, hstop):
co2_in_fn = time.strftime(cosmo_out+'/lffd%Y%m%d%H.nc')
co2_out_fn = time.strftime(cosmo_out+'/CO2_'+ens+'_%Y%m%d%H.nc')
hhl_fn = time.strftime(cosmo_out+'/lffd%Y%m%d%Hc.nc')
cdo.expr("CO2=CO2_BG-CO2_GPP+CO2_RESP+CO2_A_CH+CO2_A", input = "-selname,CO2_BG,CO2_GPP,CO2_RESP,CO2_A_CH,CO2_A "+co2_in_fn, output = co2_out_fn)
files2cat.append(co2_out_fn)
cdo.selname("HHL", input = hhl_fn, output = "hhl.nc")
cdo.cat(input = files2cat, output = "CO2_"+time_stamp+".nc")
sites = ("lhw","brm","jfj","ssl")
cdo.remapnn("lon=7.99_lat=46.54,", input = "CO2_"+time_stamp+".nc", output = "CO2_jfj_"+time_stamp+".nc")
cdo.remapnn("lon=7.99_lat=46.54,", input = "hhl.nc", output = "hhl_jfj.nc")
cdo.remapnn("lon=8.40_lat=47.48,", input = "CO2_"+time_stamp+".nc", output = "CO2_lhw_"+time_stamp+".nc")
cdo.remapnn("lon=8.40_lat=47.48,", input = "hhl.nc", output = "hhl_lhw.nc")
cdo.remapnn("lon=8.18_lat=47.19,", input = "CO2_"+time_stamp+".nc", output = "CO2_brm_"+time_stamp+".nc")
cdo.remapnn("lon=8.18_lat=47.19,", input = "hhl.nc", output = "hhl_brm.nc")
cdo.remapnn("lon=7.92_lat=47.92,", input = "CO2_"+time_stamp+".nc", output = "CO2_ssl_"+time_stamp+".nc")
cdo.remapnn("lon=7.92_lat=47.92,", input = "hhl.nc", output = "hhl_ssl.nc")
for s,ss in enumerate(sites):
site_height.main(ss, time_stamp)
cdo.intlevel("860", input = ",CO2_60lev_"+ens+"_lhw_"+time_stamp+".nc", output = "modelled_"+ens+"_lhw_"+time_stamp+".nc")
cdo.intlevel("797", input = ",CO2_60lev_"+ens+"_brm_"+time_stamp+".nc", output = "modelled_"+ens+"_brm_"+time_stamp+".nc")
cdo.intlevel("3580", input = "CO2_60lev_"+ens+"_jfj_"+time_stamp+".nc", output = "modelled_"+ens+"_jfj_"+time_stamp+".nc")
cdo.intlevel("1205", input = "CO2_60lev_"+ens+"_ssl_"+time_stamp+".nc", output = "modelled_"+ens+"_ssl_"+time_stamp+".nc")
cdo.cat(input = "modelled_"+ens+"_brm_"+time_stamp+".nc modelled_"+ens+"_jfj_"+time_stamp+".nc modelled_"+ens+"_lhw_"+time_stamp+".nc modelled_"+ens+"_ssl_"+time_stamp+".nc ", output = "model_"+ens+"_"+time_stamp+".nc")
################### End Class ObservationOperator ###################
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
"""
if __name__ == "__main__":
pass
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