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

Fixed some issues related to dates

parent e545544d
...@@ -69,6 +69,7 @@ model : ...@@ -69,6 +69,7 @@ model :
# Also projects information from the observation to the model space # Also projects information from the observation to the model space
# - fic_obsvect: observation vector from previous simulations # - fic_obsvect: observation vector from previous simulations
# - dump_debug: write out extra information (for debugging) # - dump_debug: write out extra information (for debugging)
# - sort_order: (optional) datastore sort order, ['index', 'station'] is default
# For FLEXPART the background plugin takes the following parameters # For FLEXPART the background plugin takes the following parameters
# - dir_initconc: directory containing concentration files # - dir_initconc: directory containing concentration files
# - file_initconc: name of concentration files # - file_initconc: name of concentration files
...@@ -84,6 +85,7 @@ obsvect: ...@@ -84,6 +85,7 @@ obsvect:
dump: True dump: True
dump_debug: True dump_debug: True
dump_type: nc dump_type: nc
sort_order: ['station', 'index']
background: background:
plugin: plugin:
...@@ -163,7 +165,7 @@ domain : ...@@ -163,7 +165,7 @@ domain :
nlat : 90 nlat : 90
nlev : 1 nlev : 1
type : deg type : deg
xmin_glob : -179. xmin_glob : -180.
ymin_glob : -90. ymin_glob : -90.
nlon_glob : 360 nlon_glob : 360
nlat_glob : 180 nlat_glob : 180
...@@ -189,9 +191,9 @@ mode: ...@@ -189,9 +191,9 @@ mode:
version: flexpart version: flexpart
reload_from_previous: False reload_from_previous: False
maxiter: 10 maxiter: 10
# zreduc: 0.001 df1: 1
# epsg: 0.01 nsim: 100
df1: 0.01 m: 1
save_out_netcdf: True save_out_netcdf: True
#################################################################### ####################################################################
...@@ -231,17 +233,17 @@ controlvect: ...@@ -231,17 +233,17 @@ controlvect:
hresol : regions hresol : regions
inc_ocean : true inc_ocean : true
fileregions : /home/eso/FLEX_INV/CH4/TEST_OUTPUT/regions_ghg.nc fileregions : /home/eso/FLEX_INV/CH4/TEST_OUTPUT/regions_ghg.nc
errtype : max # errtype : max
# errtype : avg # errtype : avg
err : 0.5 err : 0.5
# Lower limit flux error: unit (kg/m2/h) # Lower limit flux error: unit (kg/m2/h)
flxerr_ll: 1.e-8 flxerr_ll: 1.e-8
# Total error inversion domain: unit (Tg/y) # Total error inversion domain: unit (Tg/y)
# glob_err: 10. glob_err: 10.
numscale : 1.E12 numscale : 1.E12
xb_scale : 1. xb_scale : 1.
# periodflux : 5D
period : 1D period : 1D
# period : 10D
dir : /home/eso/FLEX_INV/CH4/TEST_OUTPUT/FLUXES/GHG/ dir : /home/eso/FLEX_INV/CH4/TEST_OUTPUT/FLUXES/GHG/
file : CH4_TOTAL_%Y_05x05.nc file : CH4_TOTAL_%Y_05x05.nc
file_glob : CH4_TOTAL_%Y_10x10.nc file_glob : CH4_TOTAL_%Y_10x10.nc
...@@ -260,15 +262,15 @@ controlvect: ...@@ -260,15 +262,15 @@ controlvect:
filelsm : /home/eso/FLEX_INV/CH4/TEST_INPUT/lsm_0.5x0.5_VERIFY.nc filelsm : /home/eso/FLEX_INV/CH4/TEST_INPUT/lsm_0.5x0.5_VERIFY.nc
dircorrel : /home/eso/repos/CIF/flextest/ dircorrel : /home/eso/repos/CIF/flextest/
dump_hcorr : True dump_hcorr : True
# sigma_land: 250. sigma_land: 250.
# sigma_sea: 1000. sigma_sea: 1000.
sigma_land: 1. # sigma_land: 1.
sigma_sea: 1. # sigma_sea: 1.
evalmin : 1.e-6 evalmin : 1.e-6
# evalmin : 1.e-15 # evalmin : 1.e-15
tcorrelations : tcorrelations :
# sigma_t : 30 # sigma_t : 30
sigma_t : 1 sigma_t : 10
dump_tcorr : True dump_tcorr : True
dircorrel : /home/eso/repos/CIF/flextest/ dircorrel : /home/eso/repos/CIF/flextest/
# evalmin : 1.e-6 # evalmin : 1.e-6
......
...@@ -80,7 +80,6 @@ def build_hcorrelations(zlat, zlon, lsm, ...@@ -80,7 +80,6 @@ def build_hcorrelations(zlat, zlon, lsm,
if regions: if regions:
landseamask2d = readnc(file_lsm, ['lsm']) landseamask2d = readnc(file_lsm, ['lsm'])
landseamask = map2scale(landseamask2d[np.newaxis, np.newaxis, :, :], tracer, tracer.domain, region_scale_area=False, region_max_val=True).flatten() landseamask = map2scale(landseamask2d[np.newaxis, np.newaxis, :, :], tracer, tracer.domain, region_scale_area=False, region_max_val=True).flatten()
import pdb; pdb.set_trace()
else: else:
landseamask = readnc(file_lsm, ['lsm']).flatten() landseamask = readnc(file_lsm, ['lsm']).flatten()
...@@ -105,8 +104,6 @@ def build_hcorrelations(zlat, zlon, lsm, ...@@ -105,8 +104,6 @@ def build_hcorrelations(zlat, zlon, lsm,
# Compute matrix of distance # Compute matrix of distance
dx = dist_matrix(zlat, zlon, projection) dx = dist_matrix(zlat, zlon, projection)
import pdb; pdb.set_trace()
# Compute the correlation matrix itself # Compute the correlation matrix itself
corr = np.exp(old_div(-dx, sigma)) corr = np.exp(old_div(-dx, sigma))
corr[sigma <= 0] = 0 corr[sigma <= 0] = 0
...@@ -133,7 +130,6 @@ def build_hcorrelations(zlat, zlon, lsm, ...@@ -133,7 +130,6 @@ def build_hcorrelations(zlat, zlon, lsm,
# mask = evalues >= evalmin # mask = evalues >= evalmin
# ESO: This is how it is done in flexinvert # ESO: This is how it is done in flexinvert
import pdb; pdb.set_trace()
mask = evalues >= evalmin*evalues.max() mask = evalues >= evalmin*evalues.max()
check.verbose("Truncating eigenvalues at " + str(evalmin*evalues.max())) check.verbose("Truncating eigenvalues at " + str(evalmin*evalues.max()))
......
...@@ -50,8 +50,8 @@ def build_tcorrelations(period, dates, sigma_t, ...@@ -50,8 +50,8 @@ def build_tcorrelations(period, dates, sigma_t,
- pd.DatetimeIndex(dates)[np.newaxis, :]), np.timedelta64(sigma_t, 'D')) - pd.DatetimeIndex(dates)[np.newaxis, :]), np.timedelta64(sigma_t, 'D'))
# Compute the correlation matrix itself # Compute the correlation matrix itself
# corr = np.exp(-dt ** 2) # corr = np.exp(-dt ** 2)
# ESO: # ESO:
corr = np.exp(-np.abs(dt)) corr = np.exp(-np.abs(dt))
# Component analysis # Component analysis
...@@ -75,8 +75,6 @@ def build_tcorrelations(period, dates, sigma_t, ...@@ -75,8 +75,6 @@ def build_tcorrelations(period, dates, sigma_t,
# Truncating values < evalmin # Truncating values < evalmin
mask = evalues >= evalmin mask = evalues >= evalmin
import pdb; pdb.set_trace()
return evalues[mask] ** 0.5, evectors[:, mask] return evalues[mask] ** 0.5, evectors[:, mask]
......
...@@ -8,7 +8,8 @@ from .utils.scalemaps import scale2map ...@@ -8,7 +8,8 @@ from .utils.scalemaps import scale2map
from pycif.utils.path import init_dir from pycif.utils.path import init_dir
def dump(self, cntrl_file, to_netcdf=False, dir_netcdf=None, run_id=None, **kwargs): def dump(self, cntrl_file, to_netcdf=False, dir_netcdf=None, run_id=None,
smallnum=1.e-27, **kwargs):
"""Dumps a control vector into a pickle file. """Dumps a control vector into a pickle file.
Does not save large correlations. Does not save large correlations.
...@@ -77,36 +78,26 @@ def dump(self, cntrl_file, to_netcdf=False, dir_netcdf=None, run_id=None, **kwar ...@@ -77,36 +78,26 @@ def dump(self, cntrl_file, to_netcdf=False, dir_netcdf=None, run_id=None, **kwar
# Make a copy to write out the original array # Make a copy to write out the original array
xa=copy.deepcopy(x) xa = copy.deepcopy(x)
xa/=np.float(self.model.numscale) xa /= np.float(self.model.numscale)
# If offsets, add prior fluxes # If offsets, add prior fluxes
if getattr(tracer, 'offsets', False): if getattr(tracer, 'offsets', False):
offsets = True offsets = True
print("Offsets=True")
x[:,0,:,:] = (x[:,0,:,:] + self.flxall)/np.float(self.model.numscale) x[:,0,:,:] = (x[:,0,:,:] + self.flxall)/np.float(self.model.numscale)
else: else:
print("Offsets=False") print("Offsets=False")
offsets = False offsets = False
# TODO: # Factor for regridding
# transform the scaling factor (0 - 1) to grid and apply to x. Check
# that sum of scaling factor is 1 and also check what to do with negative values?
# Transform scaling factor
scaling_factor = xb_grid / xb scaling_factor = xb_grid / xb
# Rescale using prior flux ratio grid to box # Rescale using prior flux ratio grid to box
x.values[xb_grid.values > 1.e-15] = x.values[xb_grid.values > 1.e-15]*scaling_factor.values[xb_grid.values > 1.e-15] x.values[xb_grid.values > smallnum] = x.values[xb_grid.values > smallnum]*scaling_factor.values[xb_grid.values > smallnum]
#x[xb_grid > 1.e-15] = x*scaling_factor
x[:,:,:,:] /= np.float(self.model.numscale) x[:,:,:,:] /= np.float(self.model.numscale)
std = np.reshape( std = np.reshape(
self.std[tracer.xpointer:tracer.xpointer + tracer.dim], self.std[tracer.xpointer:tracer.xpointer + tracer.dim],
(tracer.ndates, -1)) (tracer.ndates, -1))
......
...@@ -166,8 +166,7 @@ def init_bprod(cntrlv, options={}, **kwargs): ...@@ -166,8 +166,7 @@ def init_bprod(cntrlv, options={}, **kwargs):
tracer.domain, tracer.domain,
region_scale_area=False, region_max_val=False)) region_scale_area=False, region_max_val=False))
inv_intervals = tracer.dates[0::-1].shape[0] errsum = np.dot(cntrlv.std, np.tile(area_reg, tracer.ndates)/float(tracer.ndates))
errsum = np.dot(cntrlv.std, area_reg)/float(inv_intervals)
toterr = errsum*3600.*24.*365./1.e9/float(tracer.numscale) toterr = errsum*3600.*24.*365./1.e9/float(tracer.numscale)
cntrlv.errscalar = glob_err/toterr cntrlv.errscalar = glob_err/toterr
cntrlv.std = cntrlv.std*cntrlv.errscalar**2 cntrlv.std = cntrlv.std*cntrlv.errscalar**2
...@@ -175,7 +174,7 @@ def init_bprod(cntrlv, options={}, **kwargs): ...@@ -175,7 +174,7 @@ def init_bprod(cntrlv, options={}, **kwargs):
# TODO: scale cross correlations as well # TODO: scale cross correlations as well
# probably apply scaling to hcorr? # probably apply scaling to hcorr?
import pdb; pdb.set_trace() # import pdb; pdb.set_trace()
corr.sqrt_evalues = sqrt_evalues corr.sqrt_evalues = sqrt_evalues
corr.evectors = evectors corr.evectors = evectors
......
...@@ -67,7 +67,7 @@ def sqrtbprod(controlvect, chi, **kwargs): ...@@ -67,7 +67,7 @@ def sqrtbprod(controlvect, chi, **kwargs):
# Filling corresponding part in the control vector # Filling corresponding part in the control vector
xout[x_pointer:x_pointer + x_dim] = chi_tmp xout[x_pointer:x_pointer + x_dim] = chi_tmp
return xout * controlvect.std + controlvect.xb return xout * controlvect.std + controlvect.xb
......
...@@ -12,7 +12,7 @@ def create_domain(domain, ...@@ -12,7 +12,7 @@ def create_domain(domain,
Notes: Notes:
We assume that center coordinates are used in the domain definition We assume that outer edge coordinates are used in the domain definition
Todo: Todo:
...@@ -72,12 +72,15 @@ def create_domain(domain, ...@@ -72,12 +72,15 @@ def create_domain(domain,
domain.zlat_glob = zlat_glob domain.zlat_glob = zlat_glob
# Indices for inversion domain into global domain (relative to lower left corner) # Indices for inversion domain into global domain (relative to lower left corner)
ix1 = np.argmin(np.abs(lonc_glob - domain.lon[0] - 1)) ix1 = np.argmin(np.abs(lonc_glob - domain.lon[0] ))
iy1 = np.argmin(np.abs(latc_glob - domain.lat[0] - 1)) iy1 = np.argmin(np.abs(latc_glob - domain.lat[0] ))
ix2 = np.argmin(np.abs(lonc_glob - domain.lon[-1] - 1)) ix2 = np.argmin(np.abs(lonc_glob - domain.lon[-1] ))
iy2 = np.argmin(np.abs(latc_glob - domain.lat[-1] - 1)) iy2 = np.argmin(np.abs(latc_glob - domain.lat[-1] ))
domain.ix1 = ix1 domain.ix1 = ix1
domain.ix2 = ix2 domain.ix2 = ix2
domain.iy1 = iy1 domain.iy1 = iy1
domain.iy2 = iy2 domain.iy2 = iy2
# import pdb; pdb.set_trace()
...@@ -81,5 +81,5 @@ def read_glob(self, name, tracdir, tracfic, dates, ...@@ -81,5 +81,5 @@ def read_glob(self, name, tracdir, tracfic, dates,
flx[ddt, :, :] = xmod.interp(time=self.dates[ddt]) flx[ddt, :, :] = xmod.interp(time=self.dates[ddt])
else: else:
flx[ddt, :, :] = xmod.sel(time=self.dates[ddt], method='nearest') flx[ddt, :, :] = xmod.sel(time=self.dates[ddt], method='nearest')
return flx return flx
...@@ -37,14 +37,11 @@ def check_options(self, chi, finit, ...@@ -37,14 +37,11 @@ def check_options(self, chi, finit,
# Number of simulations # Number of simulations
# Can be explicitly given by the user, or will be determined from 'maxiter' # Can be explicitly given by the user, or will be determined from 'maxiter'
try:
niter = self.niter maxiter = self.maxiter
nsim = self.nsim niter = getattr(self, 'niter', maxiter)
nsim = getattr(self, 'nsim', 2*maxiter)
except AttributeError:
maxiter = self.maxiter
niter = maxiter
nsim = 2 * maxiter
# Default parameters # Default parameters
m = getattr(self, 'm', 5) m = getattr(self, 'm', 5)
......
...@@ -26,7 +26,7 @@ def minimize(self, finit, gradinit, chi0, **kwargs): ...@@ -26,7 +26,7 @@ def minimize(self, finit, gradinit, chi0, **kwargs):
# Initializing options (and filling missing values with default) # Initializing options (and filling missing values with default)
self = self.check_options(chi0, finit, **kwargs) self = self.check_options(chi0, finit, **kwargs)
# Running M1QN3 # Running M1QN3
xopt, fopt, gradopt, niter, nsim, epsg, mode = \ xopt, fopt, gradopt, niter, nsim, epsg, mode = \
self.m1qn3(finit, gradinit, chi0, **kwargs) self.m1qn3(finit, gradinit, chi0, **kwargs)
......
...@@ -47,10 +47,10 @@ def m1qn3(self, f, g, chi, **kwargs): ...@@ -47,10 +47,10 @@ def m1qn3(self, f, g, chi, **kwargs):
rm1 = getattr(self, 'rm1', 1e-4) rm1 = getattr(self, 'rm1', 1e-4)
rm2 = getattr(self, 'rm2', 0.99) rm2 = getattr(self, 'rm2', 0.99)
rmin = getattr(self, 'rmin', 1e-20) rmin = getattr(self, 'rmin', 1e-20)
nverbose = kwargs.get('verbose', 2) nverbose = kwargs.get('verbose', 2)
logfile = kwargs.get('logfile', None) logfile = kwargs.get('logfile', None)
# M1QN3 main inputs: # M1QN3 main inputs:
x = chi x = chi
......
...@@ -6,6 +6,8 @@ from __future__ import absolute_import ...@@ -6,6 +6,8 @@ from __future__ import absolute_import
import os import os
import numpy as np import numpy as np
import glob import glob
import re
from datetime import datetime, timedelta
from . import flexpart_header from . import flexpart_header
from . import mod_flexpart from . import mod_flexpart
...@@ -49,7 +51,6 @@ def read_flexpart_dir(subdir, nested=True, **kwargs): ...@@ -49,7 +51,6 @@ def read_flexpart_dir(subdir, nested=True, **kwargs):
grids[:,:,:,ii] = grid_fp grids[:,:,:,ii] = grid_fp
gtimes[:, ii] = gtime gtimes[:, ii] = gtime
return grids, gtimes return grids, gtimes
...@@ -86,13 +87,25 @@ def read_flexpart_grid(subdir, file_name, fp_header, **kwargs): ...@@ -86,13 +87,25 @@ def read_flexpart_grid(subdir, file_name, fp_header, **kwargs):
# Convert grid times to datetime format # Convert grid times to datetime format
gtime_dt = [] gtime_dt = []
# Convert FLEXPART julian days to datetime
for i in range(len(gtime)): for i in range(len(gtime)):
if gtime[i] == 0.: if gtime[i] == 0.:
break break
gtime_dt.append(j2d(gtime[i])) gtime_dt.append(j2d(gtime[i]))
# Generate the footprint dates as a series (for debugging), adding one hour
# gtime_dt.append(datetime.strptime(re.findall(r'\d+', file_name)[0], '%Y%m%d%H%M%S'))
# gtime_dt[0] += timedelta(hours=1)
# gtime_dt[0] = gtime_dt[0].replace(hour=0)
# for i in range(1, len(gtime)):
# if gtime[i] == 0.:
# break
# gtime_dt.append(gtime_dt[0] - timedelta(days=i))
# gtime_dt.reverse()
return grid_fp, gtime_dt, ngrid return grid_fp, gtime_dt, ngrid
......
...@@ -110,12 +110,12 @@ def obsoper(self, inputs, mode, ...@@ -110,12 +110,12 @@ def obsoper(self, inputs, mode,
iy2 = model.domain.iy2 iy2 = model.domain.iy2
else: else:
raise Exception("For FLEXPART, only hresol:regions are implemented in controlvect") raise Exception("For FLEXPART, only hresol:regions are implemented in controlvect")
# Loop through model periods and read model output # Loop through model periods and read model output
self.missingperiod = False self.missingperiod = False
# import pdb; pdb.set_trace()
for di, df in zip(subsimu_dates[:-1], subsimu_dates[1:]): for di, df in zip(subsimu_dates[:-1], subsimu_dates[1:]):
# Save to datastore for debugging purposes # Save to datastore for debugging purposes
...@@ -129,27 +129,23 @@ def obsoper(self, inputs, mode, ...@@ -129,27 +129,23 @@ def obsoper(self, inputs, mode,
# Loop over observation dates # Loop over observation dates
print("di, df", di, df) print("di, df", di, df)
for obs_i, row in enumerate(obsvect.datastore.itertuples()): for obs_i, row in enumerate(obsvect.datastore.itertuples()):
# for obs_i, row in enumerate(obsvect.datastore.iterrows()):
# Subdirectory for FLEXPART footprints # Subdirectory for FLEXPART footprints
subdir = row.Index.strftime("%Y%m") # itertuples subdir = row.Index.strftime("%Y%m")
# subdir = row[0].strftime("%Y%m")
# For debugging # For debugging
obs_check[obs_i] = obs_i obs_check[obs_i] = obs_i
station = row.station # itertuples station = row.station
# station = row[1].station
runsubdir_nest = os.path.join( runsubdir_nest = os.path.join(
model.run_dir_nest, station.upper(), subdir) model.run_dir_nest, station.upper(), subdir)
runsubdir_glob = os.path.join( runsubdir_glob = os.path.join(
model.run_dir_glob, station.upper(), subdir) model.run_dir_glob, station.upper(), subdir)
file_date = row.Index.strftime('%Y%m%d%H%M%S') # itertuples file_date = row.Index.strftime('%Y%m%d%H%M%S')
# file_date = row[0].strftime('%Y%m%d%H%M%S')
# Read nested grids # Read nested grids
file_name = 'grid_time_nest_' + file_date + '_001' file_name = 'grid_time_nest_' + file_date + '_001'
...@@ -176,15 +172,15 @@ def obsoper(self, inputs, mode, ...@@ -176,15 +172,15 @@ def obsoper(self, inputs, mode,
grid_glob *= model.coeff*model.mmair/model.molarmass grid_glob *= model.coeff*model.mmair/model.molarmass
# Background contribution from fluxes outside domain # Background contribution from fluxes outside domain
hbkg = np.sum(grid_glob[:, :, 0:ngrid-1], axis=2) hbkg = np.sum(grid_glob[:, :, 0:ngrid_glob], axis=2)
hbkg[ix1:ix2, iy1:iy2] = 0.0 hbkg[ix1:ix2, iy1:iy2] = 0.0
# Index to state vector # Index to state vector
ind = np.argmin(np.abs(tracer.dates[0::-1] - gtime_glob[0])) ind = np.argmin(np.abs(tracer.dates[0:-1] - gtime_glob[0]))
# TODO: calculate on first iteration ,then store? # TODO: calculate on first iteration ,then store?
obs_bkg[obs_i] = np.sum(hbkg[:,:].T*controlvect.flx_bkg[ind,:,:]) obs_bkg[obs_i] = np.sum(hbkg[:,:].T*controlvect.flx_bkg[ind,:,:])
obs_bkgerr[obs_i] = obs_bkg[obs_i]*tracer.err obs_bkgerr[obs_i] = np.sum(hbkg[:,:].T*np.abs(controlvect.flx_bkg[ind,:,:]))*tracer.err
# Transport for boxes in regions # Transport for boxes in regions
...@@ -204,11 +200,14 @@ def obsoper(self, inputs, mode, ...@@ -204,11 +200,14 @@ def obsoper(self, inputs, mode,
# Calculate indices to state vector # Calculate indices to state vector
for i,j in enumerate(gtime): for i,j in enumerate(gtime):
if j > df: #if j > df:
# Avoid including the midnight files at period end
if j >= df:
istate[i] = -1 istate[i] = -1
else: else:
# Discard 1st tracer date in the comparison # Discard 1st tracer date in the comparison
mask = j - tracer.dates[1::] <= datetime.timedelta(0) # (similar to writing dates + stateres)
mask = j - tracer.dates[1::] < datetime.timedelta(0)
istate[i] = int(np.argmax(mask)) istate[i] = int(np.argmax(mask))
if np.max(istate) < 0: if np.max(istate) < 0:
...@@ -235,14 +234,13 @@ def obsoper(self, inputs, mode, ...@@ -235,14 +234,13 @@ def obsoper(self, inputs, mode,
# TODO: This can possibly be calculated at first iteration only # TODO: This can possibly be calculated at first iteration only
obs_ghg[obs_i] = 0.