Commit faf58bac authored by Espen Sollum's avatar Espen Sollum
Browse files

Added missing file

parent 72f67d46
from builtins import str
import pandas as pd
import sys
import os
import shutil
import numpy as np
from datetime import timedelta
from pycif.utils.datastores.dump import dump_datastore, read_datastore
from pycif.plugins.obsvects.standard.utils.tstep import tstep
from pycif.plugins.obsvects.standard.utils.crop_monitor import crop_monitor
from pycif.utils.check import verbose
from pycif.utils import path
from pycif.utils.check.errclass import PluginError
def init_background(obsvect, **kwargs):
"""Initializes the background
Reads and regrids concentration fields
Reads sensitivity to initial conditions fields
Calculates contribution to observations from initial concentrations
obsvect (Plugin): the observation vector with all its attributes
obsvect, datastore updated with contribution to observations from initial
background = obsvect.background
model = obsvect.model
# Read global concentration fields. Regrid to model grid and
# filter by start and end date of simulation
conc_ini = background.read_conc(model, background.dir_initconc,
background.file_initconc, **kwargs)
n_x = model.domain.nlon_glob
n_y = model.domain.nlat_glob
n_z = len(model.domain.heights)
# cini is dimension number of observations x zones to optimize
# Stored in monitor as cini_1, cini_2 etc.
obsvect.n_cini = len(background.cini_lat)
obsvect.cini_names = ["cini_{}".format(i) for i in range(obsvect.n_cini)]
cini = np.empty((len(obsvect.datastore), obsvect.n_cini), order='C'); cini[:,:]=np.nan
# Read FLEXPART header to access trajdays
fp_header = None
# Loop over all observations and calculate background contribution
for obs_i, row in enumerate(obsvect.datastore.itertuples()):
# Subdirectory for FLEXPART grid_initial files
subdir = row.Index.strftime("%Y%m")
station = row.station
runsubdir_glob = os.path.join(
model.run_dir_glob, station.upper(), subdir)
file_date = row.Index.strftime('%Y%m%d%H%M%S')
file_name = 'grid_initial_' + file_date + '_001'
if not os.path.isfile(os.path.join(runsubdir_glob, file_name)):
# Read grid_initial file
gridini =, file_name,
n_x, n_y, n_z)
# Assume all trajdays are the same so read header once only
if fp_header is None:
fp_header = model.utils.flexpart_header.Flexpartheader()
fp_header.read_header(os.path.join(runsubdir_glob, 'header'))
trajdays = fp_header.trajdays
# Interpolate background concentrations to observation date minus trajdays
time_interpolate = row.Index - timedelta(days=trajdays)
conc = np.array(conc_ini.interp(coords={'time' : time_interpolate}))
for iz in range(len(model.domain.heights)):
print("background, iz:", iz)
for iy in range(model.domain.nlat_glob):
n = np.min(np.where(model.domain.lat_glob[iy] < background.cini_lat))
for ix in range(model.domain.nlon_glob):
cini[obs_i, n] += gridini[ix, iy, iz] * conc[ix, iy, iz]
# Add initial concentrations to datastore
for i, name in enumerate(obsvect.cini_names):
obsvect.datastore[name] = cini[:, i]
return obsvect
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