Commit 25dfc48e authored by Jean Matthieu Haussaire's avatar Jean Matthieu Haussaire
Browse files

Adapt the obsoper to my case. It is hardcoded, with pieces of code delete, paths hardcoded ...

parent e2539dfb
......@@ -3,11 +3,14 @@ import datetime
import numpy as np
import os
#HJM
import time
from pycif.utils import path
from pycif.utils.check import verbose
from pycif.utils.datastores.dump import dump_datastore, read_datastore
from .check import check_inputs
from netCDF4 import Dataset
def obsoper(self, inputs, mode,
run_id=0, datei=datetime.datetime(1979, 1, 1),
......@@ -77,32 +80,41 @@ def obsoper(self, inputs, mode,
# Various initializations
# Get flexpart headers
subdir = subsimu_dates[0].strftime("%Y%m")
# subdir = subsimu_dates[0].strftime("%Y%m")
fp_header_glob = model.utils.flexpart_header.Flexpartheader()
try:
fp_header_glob.read_header(
os.path.join(model.run_dir_glob, obsvect.datastore.head(1)['station'][0].decode().
upper(),subdir, 'header'))
except AttributeError as e:
fp_header_glob.read_header(
os.path.join(model.run_dir_glob, obsvect.datastore.head(1)['station'][0].
upper(),subdir, 'header'))
# fp_header_glob = model.utils.flexpart_header.Flexpartheader()
# try:
# fp_header_glob.read_header(
# os.path.join(model.run_dir_glob, obsvect.datastore.head(1)['station'][0].decode().
# upper(),subdir, 'header'))
# except AttributeError as e:
# fp_header_glob.read_header(
# os.path.join(model.run_dir_glob, obsvect.datastore.head(1)['station'][0].
# upper(),subdir, 'header'))
if model.plugin.nested:
try:
fp_header_nest = model.utils.flexpart_header.Flexpartheader()
fp_header_nest.read_header(
os.path.join(model.run_dir_nest, obsvect.datastore.head(1)['station'][0].
decode().upper(), subdir, 'header_nest'))
except AttributeError as e:
fp_header_nest.read_header(
os.path.join(model.run_dir_nest, obsvect.datastore.head(1)['station'][0].
upper(), subdir, 'header_nest'))
# if model.plugin.nested:
# try:
# fp_header_nest = model.utils.flexpart_header.Flexpartheader()
# fp_header_nest.read_header(
# os.path.join(model.run_dir_nest, obsvect.datastore.head(1)['station'][0].
# decode().upper(), subdir, 'header_nest'))
# except AttributeError as e:
# fp_header_nest.read_header(
# os.path.join(model.run_dir_nest, obsvect.datastore.head(1)['station'][0].
# upper(), subdir, 'header_nest'))
# Trajectory life time
trajdays = fp_header_glob.trajdays
#trajdays = fp_header_nest.trajdays
fp_header_nest = model.utils.flexpart_header.Flexpartheader()
subdir = subsimu_dates[0].strftime("%Y%m%d")
hjm_path = os.path.join(model.run_dir_nest, obsvect.datastore.head(1)['station'][0].upper(),subdir, 'header')
fp_header_nest.read_header(hjm_path)
# Trajectory life time
trajdays = fp_header_nest.trajdays
# Get the domain definition
species = getattr(controlvect.components, 'fluxes').parameters.attributes[0]
tracer = getattr(getattr(controlvect.components, 'fluxes').parameters, species)
......@@ -149,83 +161,70 @@ def obsoper(self, inputs, mode,
print("di, df", di, df)
for obs_i, row in enumerate(obsvect.datastore.itertuples()):
#HJM TODO add station pick
# Or split my grid file into one folder per station
subdir = row.Index.strftime("%Y%m%d")
# Subdirectory for FLEXPART footprints
subdir = row.Index.strftime("%Y%m")
# For debugging
obs_check[obs_i] = obs_i
station = row.station
try:
runsubdir_nest = os.path.join(
model.run_dir_nest, station.decode().upper(), subdir)
runsubdir_glob = os.path.join(
model.run_dir_glob, station.decode().upper(), subdir)
except AttributeError as e:
runsubdir_nest = os.path.join(
model.run_dir_nest, station.upper(), subdir)
runsubdir_glob = os.path.join(
model.run_dir_glob, station.upper(), subdir)
#d1 = row.Index+datetime.timedelta(days=1)
#file_name = d1.strftime('grid_time_%Y%m%d%H%M%S.nc')
file_name = row.Index.strftime('grid_time_%Y%m%d%H%M%S.nc')
file_date = row.Index.strftime('%Y%m%d%H%M%S')
# Read nested grids
file_name = 'grid_time_nest_' + file_date + '_001'
runsubdir_nest = os.path.join(model.run_dir_nest, station.upper(),subdir)
# If file is not found, continue
if not os.path.isfile(os.path.join(runsubdir_nest, file_name)):
verbose("WARNING: file not found: " + os.path.join(runsubdir_nest, file_name))
verbose("HJM WARNING: file not found: " + os.path.join(runsubdir_nest, file_name))
continue
grid_nest, gtime, ngrid = model.utils.read.read_flexpart_grid(
grid_nest, gtime_nest, ngrid = model.utils.read.read_flexpart_grid(
runsubdir_nest, file_name, fp_header_nest, numscale=model.numscale)
grid_nest *= model.coeff*model.mmair/model.molarmass
# Read global footprints
# TODO read correction factor dry air
file_name = 'grid_time_' + file_date + '_001'
if not os.path.isfile(os.path.join(runsubdir_glob, file_name)):
verbose("Warning: file not found, ", os.path.join(runsubdir_glob, file_name))
continue
grid_glob, gtime_glob, ngrid_glob = model.utils.read.read_flexpart_grid(
runsubdir_glob, file_name, fp_header_glob, numscale=model.numscale)
grid_glob *= model.coeff*model.mmair/model.molarmass
# Background contribution from fluxes outside domain
hbkg = np.sum(grid_glob[:, :, 0:ngrid_glob], axis=2)
hbkg[ix1:ix2, iy1:iy2] = 0.0
# Index to state vector
ind = np.argmin(np.abs(tracer.dates[0:-1] - gtime_glob[0]))
grid_nest *= model.coeff*model.mmair/model.molarmass
# TODO: calculate on first iteration ,then store?
obs_bkg[obs_i] = np.sum(hbkg[:,:].T*controlvect.flx_bkg[ind,:,:])
obs_bkgerr[obs_i] = np.sum(hbkg[:,:].T*np.abs(controlvect.flx_bkg[ind,:,:]))*tracer.err
# HJM Get the bkg contribution using the Rödenbeck approach
file_baseline = subsimu_dates[0].strftime('/store/empa/em05/projects/VERIFY/TM5-4DVAR/BASELINES_ECMWF/data/TM5-to-FLE_v1_CH4_%Y_stations.nc')
with Dataset(file_baseline) as bkg:
#find the station index in the baseline file
ids = [''.join([c.decode().lower() for c in name]) for name in bkg['identification'][:]]
station_index = np.argmax([station in name for name in ids])
# find the time index in the baseline file
ts = np.array([pd.Timestamp(t[0],t[1],t[2],t[3],t[4]) for t in bkg['time1'][:]])
try:
time_index = np.where(ts == row.Index)[0][0]
obs_bkg[obs_i] = bkg['CH4_c_trans'][time_index,station_index]
except IndexError:
print('HJM ERROR: time not found', row.Index)
#HJM TODO obs_bkgerr ?
obs_bkgerr[obs_i] = 0
# Transport for boxes in regions
print('hjm ngrid',ngrid,fp_header_nest.numy,fp_header_nest.numx)
hbox = np.zeros((ngrid,nbox), np.float64)
for jy in range(fp_header_nest.numy):
for ix in range(fp_header_nest.numx):
if tracer.regions[jy, ix] > 0:
hbox[:,tracer.regions[jy, ix] - 1] += \
grid_nest[ix, jy, :ngrid]
hbox = hbox.flatten()
reg = np.unique(tracer.regions)
for r in reg:
if r>0:
wh = np.where(tracer.regions == r)
hbox[:,r-1] = np.sum( grid_nest[wh[1],wh[0],:ngrid], axis=0)
hbox = hbox.flatten()
print('hjm hbox',np.any(hbox))
# TODO: account for possible difference nested grid/inversion domain
hnest = grid_nest
istate = np.zeros(ngrid, dtype=int) - 1
# Calculate indices to state vector
for i, j in enumerate(gtime):
for i, j in enumerate(gtime_nest):
#if j > df:
# Avoid including the midnight files at period end
if j >= df:
......@@ -242,6 +241,7 @@ def obsoper(self, inputs, mode,
# ntstep: number of inversion intervals covered by the footprints
ntstep = int(np.max(istate[istate > -1]) - np.min(istate[istate > -1]) + 1)
print('hjm nstep',ntstep,nbox)
hx = np.zeros(ntstep*nbox)
px = np.zeros(ntstep*nbox)
......@@ -254,53 +254,31 @@ def obsoper(self, inputs, mode,
px[(ns-1)*nbox:ns*nbox] = controlvect.x[istate[i]*nbox:(istate[i]+1)*nbox]
hx[(ns-1)*nbox:ns*nbox] += hbox[i*nbox:(i+1)*nbox]
print('hjm',np.any(hx),np.any(px))
obs_model[obs_i] = np.dot(hx, px)
# Change in mixing ratio from best guess estimate
# TODO: This can possibly be calculated at first iteration only
obs_ghg[obs_i] = 0.
for i, itim in enumerate(gtime):
ind = np.argmin(np.abs(tracer.dates[0:-1] - itim))
obs_ghg[obs_i] += np.sum(hnest.T[i,:,:] * controlvect.flxall[ind,:,:])
# Get background contribution
if hasattr(obsvect, 'background'):
cini = np.zeros(obsvect.background.ncini)
for i, name in enumerate(obsvect.background.cini_names):
cini[i] = getattr(row, name)
if obsvect.background.optimize_cini is False:
obs_cinipos[obs_i] = np.sum(cini[:])
else:
# number of initial mixing ratio variables to optimize
ncini = obsvect.background.ncini
# Time step of initial concentration optimization; ni
mask = row.Index - datetime.timedelta(days=trajdays) - controlvect.background.cinitime[:] - datetime.timedelta(days=controlvect.background.cini_res) < datetime.timedelta(0)
ni = int(np.argmax(mask))
# Force last time step if last cini period (mask is all False)
if all(~mask):
ni=controlvect.background.ntcini - 1
obs_cinipos[obs_i] = np.dot(cini[:],
controlvect.x[npvar+ni*ncini:npvar+(ni+1)*ncini])
else:
obs_cinipos[:]=0
if getattr(tracer, 'offsets', False):
# Optimize offsets
obs_sim[obs_i] = obs_model[obs_i] + obs_ghg[obs_i] + obs_bkg[obs_i] + obs_cinipos[obs_i]
else:
# Optimize fluxes
obs_sim[obs_i] = obs_model[obs_i] + obs_bkg[obs_i] + obs_cinipos[obs_i]
obs_sim[obs_i] = obs_model[obs_i] + obs_bkg[obs_i]
print('hjm sim:', obs_model[obs_i], obs_bkg[obs_i])
# calculate gradient in observation space
# Jo'(p) = H^TR^-1(H(p) - y)
# calculate as: Jo'(p) = sum( H_i^T*ydel_i*R_i )
# Contribution to observation gradient from obs_i
print(row)
print(row.obs)
departure = obs_sim[obs_i] - row.obs
# if background missing then observation should have no influence
......@@ -330,12 +308,12 @@ def obsoper(self, inputs, mode,
obsvect.datastore['sim'] = obs_sim
obsvect.datastore['obs_bkgerr'] = obs_bkgerr
obsvect.datastore['obs_err'] = obs_err
obsvect.datastore['obs_model'] = obs_model
if dump_debug:
obsvect.datastore['obs_ghg'] = obs_ghg
obsvect.datastore['obs_bkg'] = obs_bkg
obsvect.datastore['obs_model'] = obs_model
obsvect.datastore['obs_sim'] = obs_sim
obsvect.datastore['obs_check'] = obs_check
......
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