Added compilation of flexpart utilities to setup.py. Various developments.

examples/config_inv_flexpart_ch4.yml

+####################
+# pyCIF config file #
+#####################
+
+# Define here all parameters for pyCIF following Yaml syntax
+# For details on Yaml syntax, please see:
+# http://docs.ansible.com/ansible/latest/YAMLSyntax.html
+
+# For non-specified parameters, one can either comment the corresponding line
+# or give the value "null"
+
+# Some parameters are switches that need to be activated or not.
+# pyCIF expects a boolean True/False
+# Yaml accepts the following syntax to be converted to booleans:
+# y|Y|yes|Yes|YES|n|N|no|No|NO|true|True|TRUE|false|False|FALSE|on|On|ON|off|Off|OFF
+# This offers flexibility, but one should stick to a style for consistency
+
+#####################################################################
+#####################################################################
+# pyCIF parameters
+
+# Verbose level
+# 1 = Basic information
+# 2 = Debugging
+verbose : 2
+
+# Log file (to be saved in $wordkir)
+logfile: pyCIF.logtest
+
+# Execution directory
+workdir : /home/eso/repos/CIF/flextest/
+
+# Initial date
+# Use the following compatible format for the date:
+# YYYY-mm-dd or YYYY-mm-dd HH:mm
+datei : 2012-01-01
+
+# End date (same as for the initial date)
+datef : 2012-01-31
+#####################################################################
+#####################################################################
+# Transport model
+# Accepts any model registered in pyvar.models
+# fp_run_dir - directory from which to read FLEXPART footprints
+# nested - yes if using nested output domain
+# coeff - coefficient to convert from grid units of ppt to ppbv/ppmv (1e-3/1e-6)
+# (flexinvert uses ppb, pycif use 1e-6)
+# mmair - standard molecular mass of air
+# molarmass - Molar mass (in flux files, e.g. CO2-C=12, CH4=16)
+# numscale - numeric scaling
+
+model :
+  plugin:
+    name    : FLEXPART
+    version : std
+    periods : 1MS
+    nested  : yes
+  run_dir_nest : /home/eso/FLEX_INV/TEST_OUTPUT/FLEXOUT/GHG/NEST/
+  run_dir_glob : /home/eso/FLEX_INV/TEST_OUTPUT/FLEXOUT/GHG/NEST/
+  coeff : 1.e-3
+  mmair : 28.97
+  molarmass : 16.
+  numscale : 1.E12
+
+
+#####################################################################
+#####################################################################
+# How to build your observation vector and observation uncertainties if needed
+# Also projects information from the observation to the model space
+# - fic_obsvect: observation vector from previous simulations
+obsvect:
+  plugin:
+    name: standard
+    version: std
+    # fic_obsvect : /home/eso/repos/CIF/flextest/monitor_obsvect.nc
+    # file_obsvect : /home/eso/repos/CIF/flextest/monitor_obsvect.nc
+  dump: True
+  dump_type: nc
+#####################################################################
+#####################################################################
+
+#####################################################################
+#####################################################################
+# Measurements to account for
+# Main keys are:
+# - infos: infos on tracers to include in the simulation
+# - fic_monitor: standard txt file including all observations
+measurements :
+  plugin:
+    name: standard
+    version: std
+
+ # File where to save data, if does not exist. Reads from there if exists
+  # fic_monitor :  /home/eso/repos/CIF/flextest/monitor_allobs.nc
+  # file_monit :  /home/eso/repos/CIF/flextest/monitor_allobs.nc
+  # file_monitor :  /home/eso/repos/CIF/flextest/monitor_allobs.nc
+  dump_type : nc
+  species :
+    CH4 :
+      dailyPM :
+        provider : WDCGG
+        format   : std
+        dir_obs  : /home/eso/FLEX_INV/TEST_INPUT/OBS/GHG/*.cn.*
+#        rescale  : false
+        rescale  : true
+        na_values : -999999.99
+        default_unit : ppb
+        default_duration: 1
+        dump     : True
+# Minimum measeurement error	
+        measerr  : 5.0
+
+####################################################################
+####################################################################
+
+# spaces to control/observation spaces
+# - read_only: do not run the transport model, read output from existing run
+obsoperator:
+  plugin:
+    name: fp
+    version: std
+  autorestart: false
+  read_only: true
+#####################################################################
+#####################################################################
+
+####################################################################
+#####################################################################
+# Domain definition
+# domain :
+#   plugin :
+#     name : FLEXPART
+#     version : std
+#   xmin : -180
+#   xmax : 180
+#   ymin : -90
+#   ymax : 90
+#   nlon : 360
+#   nlat : 180
+#   nlev : 1
+#   type : deg
+
+# Domain definition
+# coordinates refers to edge of inversion domain
+# For FLEXPART, include definition of global domain
+domain :
+  plugin :
+    name : FLEXPART
+    version : std
+  xmin : -15.
+  xmax :  35.
+  ymin : 30.
+  ymax : 75.
+  nlon : 100
+  nlat : 90
+  nlev : 1
+  type : deg
+  xmin_glob : -179.
+  ymin_glob : -90.
+  nlon_glob : 360
+  nlat_glob : 180
+  dx_glob : 1.
+  dy_glob : 1.
+  
+
+
+mode:
+  plugin:
+    name: 4dvar
+    version: std
+  minimizer:
+    plugin:
+#      name: congrad
+      name: M1QN3
+      version: std
+    simulator:
+      plugin:
+        name: gausscost
+        version: flexpart
+      reload_from_previous: False
+    maxiter: 10
+#    zreduc: 0.001
+#    epsg: 0.01
+    df1: 0.01
+  save_out_netcdf: True
+
+####################################################################
+#####################################################################
+# Arguments to define the state vector
+# These are specific to the state vector and inversion method you chose
+# Please refer to the documentation to know what to include here
+# For the standard FLEXPART, include the following:
+# - ficlsm: land-sea mask (must be consistent with LMDZ grid)
+# - correl: Use correlation or not in B
+# - dircorrel: path to pre-computed correlations
+# - sigma_land: spatial correlation length for prior errors over land (km)
+# - sigma_sea: spatial correlation length for prior errors over ocean (km)
+# - tracers: list of tracers to put in the state vector (with definition arguments):
+#     - calcstd: calculate global standard deviation
+#     - resol: resolution at which fields are scaled
+#            (choice = bands,regions,pixels;
+#             if regions, provide a netcdf file ficregion
+#             if bands, define a list of latitudes as band limits (n+1 for n bands)
+#     - inc_ocean: for 'regions', optimize ocean boxes (default true)
+#     - periodflux: period of variation for increments within a month (days)
+#     - glob_err (optional) = uncertainty on global budget
+#     - glob_err (optional) = uncertainty on global budget
+controlvect:
+  plugin:
+    name: standard
+    version: std
+  components:
+    fluxes:
+      parameters:
+        CH4 :
+          plugin:
+            name: 'FLEXPART'
+            version: 'nc'
+        # resol : hpixels
+          type : flexpart
+          to_netcdf : True
+          dir_netcdf : /home/eso/repos/CIF/flextest/
+          hresol : regions
+          inc_ocean : true
+          ficregions : /home/eso/FLEX_INV/TEST_OUTPUT/regions_ghg.nc
+          fileregions : /home/eso/FLEX_INV/TEST_OUTPUT/regions_ghg.nc
+#          errtype : max
+          errtype : avg
+          err : 0.5
+#	   Lower limit flux error: unit (kg/m2/h)
+          flxerr_ll:   1.e-8
+#	   Total error inversion domain: unit (Tg/y) -- todo: not implemented
+          glob_err: 10.
+          numscale : 1.E12
+          xb_scale : 1.
+#          periodflux : 5D
+          period : 10D
+          dir : /home/eso/FLEX_INV/TEST_OUTPUT/FLUXES/GHG/
+#          fic : CH4_TOTAL_%Y_05x05.nc
+          file : CH4_TOTAL_%Y_05x05.nc
+#          fic : CH4_TOTAL_%Y_05x05.nc
+          file_glob : CH4_TOTAL_%Y_10x10.nc
+          
+          varname_flx:  emisch4
+          lonname_flx:  longitude
+          latname_flx:  latitude
+          timename_flx: time
+          offsets : True
+#          offsets : False
+
+# Initial concentrations
+          dir_initconc : /home/eso/FLEX_INV/TEST_INPUT/INIT_CONC/GHG/
+          file_initconc:  ch4_noaa_%Y%M.nc
+          varname_init:   CH4
+
+
+
+          hcorrelations :
+#            landsea: True
+            landsea: False
+            ficlsm : /home/eso/FLEX_INV/TEST_INPUT/lsm_0.5x0.5_VERIFY.nc
+            filelsm : /home/eso/FLEX_INV/TEST_INPUT/lsm_0.5x0.5_VERIFY.nc
+            dircorrel : /home/eso/repos/CIF/flextest/
+            dump_hcorr : True
+            sigma_land: 250
+            sigma_sea: 1000
+          tcorrelations :
+            sigma_t : 30
+            dump_tcorr : True
+            dircorrel  : /home/eso/repos/CIF/flextest/
+            evalmin : 1.e-5
+

pycif/plugins/controlvects/standard/__init__.py

@@ -65,14 +65,14 @@ def ini_data(plugin, **kwargs):
 
                 # Gets read from model if not already defined
                 try:
-                    tracer.read(*range(4))
+                    tracer.read(*list(range(4)))
+#                    tracer.read(*range(4))
 #                except PluginError:
                 # except TypeError:
                 #     tracer.read(*list(range(4)))
                 except PluginError:
                     if comp in ['fluxes', 'inicond', 'latcond', 'topcond']:
                         tracer.read = getattr(plugin.model, comp).read
-                        print "ESO: TypeError"
 
                     else:
                         raise Exception('WARNING: {} was not defined in the '
@@ -80,7 +80,6 @@ def ini_data(plugin, **kwargs):
                                         'control vector'.format(comp))
                 
                 except Exception:
-                    print "ESO: exception"
                     pass
 
                 # Gets the domain and

pycif/plugins/controlvects/standard/build_hcorr.py

@@ -61,7 +61,7 @@ def build_hcorrelations(zlat, zlon, lsm,
     except IOError:
         check.verbose("Computing hcorr")
 
-        # 'regions' is WIP, for now use identity matrix
+        # TODO 'regions'
         if regions:
             if hresoldim is None:
                 raise Exception("hresoldim missing!")
@@ -124,7 +124,7 @@ def dump_hcorr(nlon, nlat, sigma_sea, sigma_land,
 
     ncell = evalues.size
 
-    # ESO TODO: change file name for spatial regions
+    # TODO change file name for spatial regions
     file_dump = '{}/horcor_{}x{}_cs{}_cl{}.bin'.format(
         dir_dump, nlon, nlat, sigma_sea, sigma_land)
     
@@ -150,7 +150,7 @@ def read_hcorr(nlon, nlat, sigma_sea, sigma_land, dir_dump, hresoldim=None):
     
     """
 
-# ESO TODO fix file name
+# TODO file name
     file_dump = '{}/horcor_{}x{}_cs{}_cl{}.bin'.format(
         dir_dump, nlon, nlat, sigma_sea, sigma_land)
     

pycif/plugins/controlvects/standard/build_tcorr.py

@@ -51,9 +51,6 @@ def build_tcorrelations(period, dates, sigma_t,
         
         # Compute the correlation matrix itself
         corr = np.exp(-dt ** 2)
-
-        # ESO: test to compare with FIP
-        # corr = np.exp(-np.abs(dt))
         
         # Component analysis
         evalues, evectors = np.linalg.eigh(corr)

pycif/plugins/controlvects/standard/dump.py

@@ -81,7 +81,7 @@ def dump(self, cntrl_file, to_netcdf=False, dir_netcdf=None, run_id=None, **kwar
             xb /= np.float(self.model.numscale)
                 
 
-            # ESO: unaggregated fluxes, for testing
+            # store the unaggregated fluxes, for testing
             arr = self.flxall[:, np.newaxis, :, :]/np.float(self.model.numscale)
             xb_grid = xr.DataArray(arr,
                             coords={'time': tracer.dates[:-1], 'lev': tracer.levels},

pycif/plugins/controlvects/standard/init_bprod.py

@@ -32,7 +32,6 @@ def init_bprod(cntrlv, options={}, **kwargs):
             tracer = getattr(component.parameters, trcr)
             
             # Deals with horizontal correlations if any
-            # ESO: check what to do with regions
             if hasattr(tracer, 'hcorrelations') and tracer.hresol == 'hpixels':
                 corr = tracer.hcorrelations
                 dump_hcorr = getattr(corr, 'dump_hcorr', False)
@@ -86,9 +85,8 @@ def init_bprod(cntrlv, options={}, **kwargs):
                 tracer.chi_hresoldim = sqrt_evalues.size
 
             elif hasattr(tracer, 'hcorrelations') and tracer.hresol == 'regions':
-                # ESO: work in progress on regions
-                # TODO possibly don't need a separate if block for regions, just
-                # pass a regions / tracer.hresol keyword to buil_hcorrelations
+                # TODO don't need a separate block for regions, instead
+                # pass tracer.hresol/regions keyword to build_hcorr
                 corr = tracer.hcorrelations
                 dump_hcorr = getattr(corr, 'dump_hcorr', False)
                 dircorrel = getattr(corr, 'dircorrel', '')

pycif/plugins/controlvects/standard/init_xb.py

@@ -43,7 +43,7 @@ def init_xb(cntrlv, **kwargs):
                     datei, datef, getattr(tracer, 'period', ''),
                     subperiod=getattr(tracer, 'subperiod', ''))
 
-            # ESO: n dates, but n-1 intervals            
+            # n dates, but n-1 intervals. TODO make consistent with rest of CIF
             tracer.ndates = len(tracer.dates) - 1
             
             # Keeping a pointer to the correct location in the whole control
@@ -64,9 +64,8 @@ def init_xb(cntrlv, **kwargs):
 
             # Filling Xb and uncertainties
             if getattr(tracer, 'type') == 'flexpart':
-                # ESO: Reading fluxes and nested fluxes for flexpart
+                # Reading fluxes and nested fluxes for flexpart
                 # For now, only aggregation to regions is supported
-                # TODO: correlation/error cov
                 flxall = tracer.read(
                     trcr, tracer.dir, tracer.file, tracer.dates,
                     comp_type=comp,
@@ -101,7 +100,6 @@ def init_xb(cntrlv, **kwargs):
                     **kwargs) * xb_scale \
                     + xb_value
 
-#                import pdb; pdb.set_trace()
                 cntrlv.flx_bkg = flx_bkg
                 
 

pycif/plugins/controlvects/standard/sqrtbprod.py

@@ -24,8 +24,7 @@ def sqrtbprod(controlvect, chi, **kwargs):
             ndates = tracer.ndates
             
             # Dealing with non-diagonal spatial B
-#            if tracer.hresol == 'hpixels' and hasattr(tracer, 'hcorrelations'):
-# ESO: allow use with regions
+            # Adding 'regions' eventually
             if (tracer.hresol == 'hpixels' or tracer.hresol == 'regions') and hasattr(tracer, 'hcorrelations'):
                 corr = tracer.hcorrelations
                 sqrt_evalues = corr.sqrt_evalues
@@ -91,8 +90,6 @@ def sqrtbprod_ad(controlvect, dx, **kwargs):
             chi_dim = tracer.chi_dim
             ndates = tracer.ndates
 
-            # ESO TODO Look at this...
-            # x * std
             xstd = dx[x_pointer: x_pointer + x_dim] \
                    * controlvect.std[x_pointer: x_pointer + x_dim]
             

pycif/plugins/controlvects/standard/utils/scalemaps.py

@@ -116,7 +116,6 @@ def map2scale(xmap, tracer, dom, scale_area=False):
                 iband += 1
 
     elif tracer.hresol == 'regions':
-        # TODO: add option to scale with area size
         x = np.zeros((ndates, tracer.vresoldim, tracer.nregions))
         if scale_area:
             for regID, reg in enumerate(np.unique(tracer.regions)):
@@ -128,8 +127,6 @@ def map2scale(xmap, tracer, dom, scale_area=False):
                 x[..., regID] = np.sum(xmap[..., tracer.regions == reg], axis=2)
 
 
-#        import pdb; pdb.set_trace()
-
     elif tracer.hresol == 'global':
         x = xmap.sum(axis=(2, 3))[..., np.newaxis]
 

pycif/plugins/domains/flexpart/__init__.py

+from create_domain import create_domain
+from read_domain import read_grid
+
+domain_flexpart = True

pycif/plugins/domains/flexpart/create_domain.py

+import numpy as np
+import pycif.utils.check as check
+
+def create_domain(domain,
+                  **kwargs):
+    """Creates a grid if needed
+
+    Args:
+        domain (dictionary): dictionary defining the domain.
+
+    Returns:
+         
+
+    Notes:
+        We assume that center coordinates are used in the domain definition 
+
+    Todo:
+        
+
+    """
+
+    nlon = domain.nlon
+    nlat = domain.nlat
+
+    dlon = (domain.xmax - domain.xmin)/nlon
+    dlat = (domain.ymax - domain.ymin)/nlat
+
+    # Corner coordinates
+    lonc = np.arange(domain.xmin, domain.xmin + nlon*dlon, dlon)
+    latc = np.arange(domain.ymin, domain.ymin + nlat*dlat, dlat)
+    
+    # Center coordinates
+    lon = lonc + dlon/2.
+    lat = latc + dlat/2.
+
+    # Meshgrids
+    zlon, zlat = np.meshgrid(lon, lat)
+
+    lonc = np.append(lonc, lonc[-1] + dlon)
+    latc = np.append(latc, latc[-1] + dlat)
+    
+    zlonc, zlatc = np.meshgrid(lonc, latc)
+
+    domain.zlon = zlon
+    domain.zlat = zlat
+    domain.zlonc = zlonc
+    domain.zlatc = zlatc
+    domain.lon = lon
+    domain.lat = lat
+
+    # Adding global grid definition
+    xmin_glob = domain.xmin_glob
+    ymin_glob = domain.ymin_glob
+    nlon_glob = domain.nlon_glob
+    nlat_glob = domain.nlat_glob
+    dlon_glob = domain.dx_glob
+    dlat_glob = domain.dy_glob
+
+    # Corner coordinates
+    lonc_glob = np.arange(xmin_glob, xmin_glob + nlon_glob*dlon_glob, dlon_glob)
+    latc_glob = np.arange(ymin_glob, ymin_glob + nlat_glob*dlat_glob, dlat_glob)
+    
+    # Center coordinates
+    lon_glob = lonc_glob + dlon_glob/2.
+    lat_glob = latc_glob + dlat_glob/2.
+
+    # Indices for inversion domain into global domain (relative to lower left corner)
+    ix1 = np.argmin(np.abs(lonc_glob - domain.lon[0] - 1))
+    iy1 = np.argmin(np.abs(latc_glob - domain.lat[0] - 1))
+    ix2 = np.argmin(np.abs(lonc_glob - domain.lon[-1] - 1))
+    iy2 = np.argmin(np.abs(latc_glob - domain.lat[-1] - 1))
+
+    domain.ix1 = ix1
+    domain.ix2 = ix2
+    domain.iy1 = iy1
+    domain.iy2 = iy2

pycif/plugins/domains/flexpart/read_domain.py

+import numpy as np
+import pycif.utils.check as check
+
+def read_grid(domain,
+              **kwargs):
+    """Reads a grid from an existing file
+
+    Args:
+        domain (Plugin): dictionary defining the domain. Should include
+        ficgrid to be able to read the grid from a file
+
+    Return:
+        Grid dictionary with meshgrids for center lon/lat and corner lon/lat
+
+    Notes:
+        For now, this function is not used, so create_domain() will be called.
+
+
+    """
+
+
+    try:
+        zlon = np.loadtxt(domain.ficlon)
+        zlat = np.loadtxt(domain.ficlat)
+    
+        nlon = zlon.size
+        nlat = zlat.size
+        
+        # Corner coordinates
+        dlon = np.ptp(zlon) / (nlon - 1) / 2.
+        zlonc = zlon - dlon
+        zlonc = np.append(zlonc, zlonc[-1] + 2 * dlon)
+        
+        dlat = np.ptp(zlat) / (nlat - 1) / 2.
+        zlatc = zlat - dlat
+        zlatc = np.append(zlatc, zlatc[-1] + 2 * dlat)
+        
+        # Meshgrids
+        zlon, zlat = np.meshgrid(zlon, zlat)
+        zlonc, zlatc = np.meshgrid(zlonc, zlatc)
+        
+        # Saving information to domain attributes
+        domain.nlon = nlon
+        domain.nlat = nlat
+        domain.zlon = zlon
+        domain.zlat = zlat
+        domain.zlonc = zlonc
+        domain.zlatc = zlatc
+    
+    except (IOError, AttributeError):
+        check.verbose("Couldn't read longitudes and latitudes.\n"
+                      "Making them from given coordinates in config file")
+        domain.create_domain()
+    

pycif/plugins/fluxes/flexpart/__init__.py

+import os
+from types import MethodType
+
+from pycif.utils import path
+from pycif.utils.check import verbose
+
+from read import read
+from read_glob import read_glob
+from write import write
+
+requirements = {'domain': {'name': 'FLEXPART', 'version': 'std', 'empty': False}}
+

pycif/plugins/fluxes/flexpart/read.py

+import os
+import pandas as pd
+import pycif.utils.check as check
+from pycif.utils import path
+import xarray as xr
+import datetime
+import numpy as np
+from pycif.utils.netcdf import readnc
+from pycif.utils.dates import j2d
+
+
+def read(self, name, tracdir, tracfic, dates,
+                interpol_flx=False, **kwargs):
+    """Get fluxes from pre-computed fluxes and load them into a pycif
+    variables
+
+    Args:
+        self: the model Plugin
+        name: the name of the component
+        tracdir, tracfic: flux directory and file format
+        dates: list of dates to extract
+        interpol_flx (bool): if True, interpolates fluxes at time t from
+        values of surrounding available files
+
+    """
+
+    # Available files in the directory