Commit 7de798f1 authored by brunner's avatar brunner
Browse files

No commit message

No commit message
parent fc970343
...@@ -91,7 +91,7 @@ class Obs(Observations): ...@@ -91,7 +91,7 @@ class Obs(Observations):
for ncfile in ncfilelist: for ncfile in ncfilelist:
infile = os.path.join(self.obspack_dir, 'data', ncfile + '.nc') infile = os.path.join(self.obspack_dir, 'data', ncfile + '.nc')
print(infile) # print(infile)
ncf = io.ct_read(infile, 'read') ncf = io.ct_read(infile, 'read')
idates = ncf.get_variable('time_components') idates = ncf.get_variable('time_components')
# print(idates) # print(idates)
...@@ -128,6 +128,9 @@ class Obs(Observations): ...@@ -128,6 +128,9 @@ class Obs(Observations):
lats = empty(shape=(len(obspacknum))) lats = empty(shape=(len(obspacknum)))
lons = empty(shape=(len(obspacknum))) lons = empty(shape=(len(obspacknum)))
alts = empty(shape=(len(obspacknum))) alts = empty(shape=(len(obspacknum)))
print(lats_1d.shape)
print(lats.shape)
# sys.exit()
lats[:] = lats_1d lats[:] = lats_1d
lons[:] = lons_1d lons[:] = lons_1d
alts[:] = alts_1d alts[:] = alts_1d
......
...@@ -117,22 +117,23 @@ class ObservationOperator(object): ...@@ -117,22 +117,23 @@ class ObservationOperator(object):
# for ncfile in ncfilelist: # for ncfile in ncfilelist:
# infile = os.path.join(ncfile + '.nc') # infile = os.path.join(ncfile + '.nc')
#for dt in rrule.rrule(rrule.HOURLY, dtstart=dacycle['time.start'], until=dacycle['time.start']+timedelta(hours=168)): for dt in rrule.rrule(rrule.HOURLY, dtstart=dacycle['time.start'], until=dacycle['time.start']+timedelta(hours=504)):
# for ens in range(0,self.forecast_nmembers): for ens in range(0,self.forecast_nmembers):
# ens = str(ens).zfill(3) ens = str(ens).zfill(3)
# cdo.setunit("'kg m-2 s-1' -expr,GPP_"+ens+"_F=CO2_GPP_F*parametermap -merge "+os.path.join(dacycle['da.bio.input'], 'gpp_%s.nc' % dt.strftime('%Y%m%d%H')), input = "/scratch/snx3000/parsenov/projname/input/parameters."+ens+".nc", output = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_"+ens+"_%s.nc" % dt.strftime('%Y%m%d%H'))) cdo.setunit("'kg m-2 s-1' -expr,GPP_"+ens+"_F=CO2_GPP_F*parametermap -merge "+os.path.join(dacycle['da.bio.input'], 'gpp_%s.nc' % dt.strftime('%Y%m%d%H')), input = os.path.join(dacycle['dir.da_run'],"input/parameters."+ens+".nc"), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_"+ens+"_%s.nc" % dt.strftime('%Y%m%d%H')))
# cdo.setunit("'kg m-2 s-1' -expr,RESP_"+ens+"_F=CO2_RESP_F*parametermap -merge "+os.path.join(dacycle['da.bio.input'], 'ra_%s.nc' % dt.strftime('%Y%m%d%H')), input = "/scratch/snx3000/parsenov/projname/input/parameters."+ens+".nc", output = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_"+ens+"_%s.nc" % dt.strftime('%Y%m%d%H'))) cdo.setunit("'kg m-2 s-1' -expr,RESP_"+ens+"_F=CO2_RESP_F*parametermap -merge "+os.path.join(dacycle['da.bio.input'], 'ra_%s.nc' % dt.strftime('%Y%m%d%H')), input = os.path.join(dacycle['dir.da_run'],"input/parameters."+ens+".nc"), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_"+ens+"_%s.nc" % dt.strftime('%Y%m%d%H')))
# cdo.merge(input = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_???_%s.nc" % dt.strftime('%Y%m%d%H')), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_%s.nc" % dt.strftime('%Y%m%d%H'))) cdo.merge(input = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_???_%s.nc" % dt.strftime('%Y%m%d%H')), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_%s.nc" % dt.strftime('%Y%m%d%H')))
# cdo.merge(input = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_???_%s.nc" % dt.strftime('%Y%m%d%H')), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_%s.nc" % dt.strftime('%Y%m%d%H'))) cdo.merge(input = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_???_%s.nc" % dt.strftime('%Y%m%d%H')), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_%s.nc" % dt.strftime('%Y%m%d%H')))
# os.chdir(dacycle['da.obsoperator.home']) os.chdir(dacycle['da.obsoperator.home'])
# os.system('python run_chain.py ctdas '+dacycle['time.start'].strftime('%Y-%m-%d')+' 0 '168' -j meteo icbc emissions biofluxes int2lm post_int2lm cosmo') # os.system('python run_chain.py ctdas '+dacycle['time.start'].strftime('%Y-%m-%d')+' 0 504 -j meteo icbc emissions biofluxes int2lm post_int2lm cosmo')
os.chdir(dacycle['dir.da_run'])
self.extract_model_data(dacycle,self.forecast_nmembers) self.extract_model_data(dacycle,self.forecast_nmembers)
for i in range(0,self.forecast_nmembers): for i in range(0,self.forecast_nmembers):
idx = str(i+1).zfill(3) idx = str(i+1).zfill(3)
cosmo_file = os.path.join('/store/empa/em05/parsenov/cosmo_data/model_'+idx+'_%s.nc' % self['time.sample.stamp']) cosmo_file = os.path.join('/store/empa/em05/parsenov/cosmo_data/model_'+idx+'_%s.nc' % dacycle['time.sample.stamp'])
ifile = Dataset(cosmo_file, mode='r') ifile = Dataset(cosmo_file, mode='r')
model_data[i,:] = np.squeeze(ifile.variables['CO2'][:])#*29./44.)#*1000000. # in ppm model_data[i,:] = np.squeeze(ifile.variables['CO2'][:])*29./44.)*1E6 # in ppm
ifile.close() ifile.close()
for j,data in enumerate(zip(ids,obs,mdm)): for j,data in enumerate(zip(ids,obs,mdm)):
...@@ -164,12 +165,11 @@ class ObservationOperator(object): ...@@ -164,12 +165,11 @@ class ObservationOperator(object):
for ens in range(0,ensnum): for ens in range(0,ensnum):
ens = str(ens).zfill(3) ens = str(ens).zfill(3)
files2cat=[] files2cat=[]
for dt in rrule.rrule(rrule.HOURLY, dtstart=dacycle['time.start'], until=dacycle['time.start']+timedelta(hours=24)): for dt in rrule.rrule(rrule.HOURLY, dtstart=dacycle['time.start'], until=dacycle['time.start']+timedelta(hours=168)):
#for dt in rrule.rrule(rrule.HOURLY, dtstart=dacycle['time.start'], until=dacycle['time.start']+timedelta(hours=168)):
dt=dt.strftime('%Y%m%d%H') dt=dt.strftime('%Y%m%d%H')
co2_in_fn = cosmo_out+'lffd'+dt+'.nc' co2_in_fn = cosmo_out+'lffd'+dt+'.nc'
co2_out_fn = cosmo_out+'CO2_'+ens+'_'+dt+'.nc' co2_out_fn = cosmo_out+'CO2_'+ens+'_'+dt+'.nc'
cdo.expr("CO2=CO2_BG-GPP_"+ens+"+RESP_"+ens+"+CO2_A_CH+CO2_A", input = "-selname,CO2_BG,GPP_"+ens+",RESP_"+ens+",CO2_A_CH,CO2_A "+co2_in_fn, output = co2_out_fn) cdo.expr("CO2=(CO2_BG-GPP_"+ens+"+RESP_"+ens+"+CO2_A_CH+CO2_A)/(1.-QV)", input = "-selname,QV,CO2_BG,GPP_"+ens+",RESP_"+ens+",CO2_A_CH,CO2_A "+co2_in_fn, output = co2_out_fn)
files2cat.append(co2_out_fn) files2cat.append(co2_out_fn)
cdo.cat(input = files2cat, output = cosmo_out+"CO2_"+ens+"_"+time_stamp+".nc") cdo.cat(input = files2cat, output = cosmo_out+"CO2_"+ens+"_"+time_stamp+".nc")
......
...@@ -26,6 +26,7 @@ import logging ...@@ -26,6 +26,7 @@ import logging
import datetime as dtm import datetime as dtm
#from string import strip #from string import strip
import numpy as np
from numpy import array, logical_and, sqrt from numpy import array, logical_and, sqrt
sys.path.append(os.getcwd()) sys.path.append(os.getcwd())
sys.path.append('../../') sys.path.append('../../')
...@@ -69,7 +70,7 @@ class ObsPackObservations(Observations): ...@@ -69,7 +70,7 @@ class ObsPackObservations(Observations):
# Step 1: Read list of available site files in package # Step 1: Read list of available site files in package
infile = os.path.join(self.obspack_dir, 'summary', '%s_dataset_summary.txt' % (self.obspack_id,)) infile = os.path.join(self.obspack_dir, '%s_dataset_summary.txt' % (self.obspack_id,))
f = open(infile, 'r') f = open(infile, 'r')
lines = f.readlines() lines = f.readlines()
f.close() f.close()
...@@ -89,7 +90,7 @@ class ObsPackObservations(Observations): ...@@ -89,7 +90,7 @@ class ObsPackObservations(Observations):
for ncfile in ncfilelist: for ncfile in ncfilelist:
infile = os.path.join(self.obspack_dir, 'data', 'nc', ncfile + '.nc') infile = os.path.join(self.obspack_dir, 'data', ncfile + '.nc')
ncf = io.ct_read(infile, 'read') ncf = io.ct_read(infile, 'read')
idates = ncf.get_variable('time_components') idates = ncf.get_variable('time_components')
dates = array([dtm.datetime(*d) for d in idates]) dates = array([dtm.datetime(*d) for d in idates])
...@@ -105,16 +106,37 @@ class ObsPackObservations(Observations): ...@@ -105,16 +106,37 @@ class ObsPackObservations(Observations):
if 'ccggAllData' in ncfile: if 'ccggAllData' in ncfile:
obspackid = ncf.get_variable('id').take(subselect, axis=0) obspackid = ncf.get_variable('id').take(subselect, axis=0)
else: else:
obspackid = ncf.get_variable('obspack_id').take(subselect, axis=0) # obspackid = ncf.get_variable('obspack_id').take(subselect, axis=0)
obspackid = [s.tostring().lower() for s in obspackid] obspackid = ['obs_' + str(int(s)) for s in obspacknum]
obspackid = list(map(str.strip,str(obspackid))) # obspackid = [s.tostring().lower() for s in obspackid]
# obspackid = list(map(str.strip,str(obspackid)))
datasetname = ncfile # use full name of dataset to propagate for clarity datasetname = ncfile # use full name of dataset to propagate for clarity
lats = ncf.get_variable('latitude').take(subselect, axis=0) if (ncf.get_variable('latitude')).shape==(1,):
lons = ncf.get_variable('longitude').take(subselect, axis=0) lats = np.empty(np.size(obspackid))
alts = ncf.get_variable('altitude').take(subselect, axis=0) lats[:] = ncf.get_variable('latitude')
else:
lats = ncf.get_variable('latitude').take(subselect, axis=0)
if (ncf.get_variable('longitude')).shape==(1,):
lons = np.empty(np.size(obspackid))
lons[:] = ncf.get_variable('longitude')
else:
lons = ncf.get_variable('longitude').take(subselect, axis=0)
if (ncf.get_variable('altitude')).shape==(1,):
alts = np.empty(np.size(obspackid))
alts[:] = ncf.get_variable('altitude')
else:
alts = ncf.get_variable('altitude').take(subselect, axis=0)
obs = ncf.get_variable('value').take(subselect, axis=0) obs = ncf.get_variable('value').take(subselect, axis=0)
species = ncf.get_attribute('dataset_parameter') if (obs.max()<1.): # pavle fix this: if units are mol mol-1
flags = ncf.get_variable('obs_flag').take(subselect, axis=0) obs = obs*1E6
species = 'co2'
# species = ncf.get_attribute('dataset_parameter')
# flags = ncf.get_variable('obs_flag').take(subselect, axis=0)
flags = np.empty(np.size(obspackid))
flags[:] = 1
ncf.close() ncf.close()
for n in range(len(dates)): for n in range(len(dates)):
...@@ -325,10 +347,10 @@ class ObsPackObservations(Observations): ...@@ -325,10 +347,10 @@ class ObsPackObservations(Observations):
for obs in self.datalist: # first loop over all available data points to set flags correctly for obs in self.datalist: # first loop over all available data points to set flags correctly
obs.mdm = 1E-6 # default is very high model-data-mismatch, until explicitly set by script obs.mdm = 0.1 # default is very high model-data-mismatch, until explicitly set by script
obs.flag = 1 obs.flag = 1
identifier = obs.code identifier = obs.code
species, site, method, lab, datasetnr = identifier.split('_') # species, site, method, lab, datasetnr = identifier.split('_')
# nr_obs_per_day = len([c.code for c in self.datalist if c.code == obs.code and c.xdate.day == obs.xdate.day and c.flag == 0]) # nr_obs_per_day = len([c.code for c in self.datalist if c.code == obs.code and c.xdate.day == obs.xdate.day and c.flag == 0])
# obs.mdm = site_info[identifier]['error'] * sqrt(nr_obs_per_day) * self.global_R_scaling # obs.mdm = site_info[identifier]['error'] * sqrt(nr_obs_per_day) * self.global_R_scaling
# obs.may_localize = site_info[identifier]['may_localize'] # obs.may_localize = site_info[identifier]['may_localize']
......
...@@ -293,6 +293,7 @@ class StateVector(object): ...@@ -293,6 +293,7 @@ class StateVector(object):
newmember = EnsembleMember(member) newmember = EnsembleMember(member)
newmember.param_values = np.dot(C, rands) + newmean newmember.param_values = np.dot(C, rands) + newmean
newmember.param_values[-1] = 0
self.ensemble_members[lag].append(newmember) self.ensemble_members[lag].append(newmember)
logging.debug('%d new ensemble members were added to the state vector # %d' % (self.nmembers, (lag + 1))) logging.debug('%d new ensemble members were added to the state vector # %d' % (self.nmembers, (lag + 1)))
......
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