Inversions in the project VERIFY

Methane net total emissions in Europe (WP4)

  • inversions with surface stations:

    1. constraints are put into a CIF monitor file for a set of core stations:

    Show/Hide Code

      1#####################################################
      2# This code is provided by: I. Pison
      3# contact: isabelle.pison@lsce.ipsl.fr
      4# date: June 2022
      5# project: H2020 VERIFY
      6# relation to CIF: uses plugins, not embedded in.
      7# It is an example, to be used as a basis
      8# for developing other work.
      9# No guarantee as to its compatibility or results
     10# is attached since it is not part of the CIF
     11# Please don't forget acknowledgements if you use it.
     12#####################################################
     13
     14import pandas as pd
     15import numpy as np
     16from pycif.utils.datastores.dump import dump_datastore, read_datastore
     17from pycif.utils.datastores.empty import init_empty
     18import glob
     19import datetime
     20import sys
     21
     22dir_data = '/home/chimereicos/VERIFYdataWP4/Obs2021/core/'
     23spec = 'ch4'
     24max_unc = 2000.
     25
     26# Initializes an empty datastore
     27ds = pd.DataFrame()
     28
     29# put errors from Szenasi et al., 2021: 
     30# A pragmatic protocol for characterising errors 
     31# in atmospheric inversions of methane emissions over Europe, 
     32# Tellus B: Chemical and Physical Meteorology, 
     33# 10.1080/16000889.2021.1914989 
     34TNO_obs_error_file = 'Obs_error_file.txt'
     35f = open(TNO_obs_error_file, 'r')
     36TNO_obs_error = f.readlines()[1:]
     37f.close()
     38TNO_error_station = []
     39TNO_station_category = []
     40TNO_LBC_error = np.empty([len(TNO_obs_error)])
     41TNO_T_error = np.empty([len(TNO_obs_error)])
     42TNO_repr_error = np.empty([len(TNO_obs_error)])
     43for ii in range(len(TNO_obs_error)):
     44    TNO_obs_error_list = TNO_obs_error[ii].split()
     45    TNO_error_station.append(TNO_obs_error_list[0])
     46    TNO_station_category.append(TNO_obs_error_list[1])
     47    if TNO_obs_error_list[2] != 'NaN' : 
     48        TNO_LBC_error[ii] = int(TNO_obs_error_list[2])
     49        TNO_T_error[ii] = int(TNO_obs_error_list[3])
     50        TNO_repr_error[ii] = int(TNO_obs_error_list[4])
     51    else:
     52        TNO_LBC_error[ii] = np.nan
     53        TNO_T_error[ii] = np.nan
     54        TNO_repr_error[ii] = np.nan
     55# squares of errors
     56montain_seor_mean = np.empty([len(TNO_obs_error)])
     57costal_seor_mean = np.empty([len(TNO_obs_error)])
     58other_seor_mean = np.empty([len(TNO_obs_error)])
     59for jj in range(len(TNO_obs_error)):
     60    if TNO_station_category[jj] == 'M':
     61        montain_seor_mean[jj] = np.sqrt((TNO_LBC_error[jj])**2 +\
     62                                    (TNO_T_error[jj])**2 +\
     63                                    (TNO_repr_error[jj])**2)
     64    else:
     65        montain_seor_mean[jj] = np.nan
     66    if TNO_station_category[jj] == 'C':
     67        costal_seor_mean[jj] = np.sqrt((TNO_LBC_error[jj])**2 +\
     68                                   (TNO_T_error[jj])**2 +\
     69                                   (TNO_repr_error[jj])**2)
     70    else:
     71        costal_seor_mean[jj] = np.nan
     72    if TNO_station_category[jj] == 'O':
     73        other_seor_mean[jj] = np.sqrt((TNO_LBC_error[jj])**2 +\
     74                                  (TNO_T_error[jj])**2 +\
     75                                  (TNO_repr_error[jj])**2)
     76    else:
     77        other_seor_mean[jj] = np.nan
     78# Mountain station or not, according to S. Houweling information
     79dt = pd.read_csv('/home/users/ipison/CHIMERE/VERIFY/intercomp_EU/CH4_observation_sitelist.csv', sep=',')
     80
     81# Read VERIFY provided files
     82list_files = glob.glob("{}/{}_*.csv".format(dir_data,spec))
     83nb_header_lines = 9
     84
     85
     86for fi in list_files:
     87
     88   fc=open(fi,'r')
     89   alllines = fc.readlines()
     90   data = alllines[nb_header_lines - 1]
     91   columns =  data.split(' ')
     92   fc.close()
     93   station = fi.split('_')[nbstat].upper()
     94   # time zone for selecting night/day hours
     95   for l in alllines[:nb_header_lines]:
     96      info = l.split(' ')[0]
     97      if info == 'timezone':
     98         utc = l.split('+')
     99         if len(utc) > 1:
    100            shift = int(utc[1])
    101         else:
    102            shift = 0
    103   df = pd.read_csv(fi, delimiter=',', header= nb_header_lines-1, names = columns) 
    104   # Information to put into the monitor file:
    105   # date
    106   df.loc[df['minute'] > 60, 'minute'] = 0
    107   df['Time']=df['year'].astype(str)+'/'+\
    108         df['month'].astype(str)+'/'+\
    109         df['day'].astype(str)+'/'+df['hour'].astype(str)+'/'+df['minute'].astype(str)
    110   df['date'] = pd.to_datetime(df['Time'],format='%Y/%m/%d/%H/%M')
    111   # duration (hours)
    112   df = df.assign(duration = 1.) 
    113   # station
    114   df.rename(columns={'siteid': 'station', 'longitude': 'lon', 'latitude': 'lat'}, inplace=True)
    115   # filter on time and type of station
    116   # 12:00 to 16:00 local time for not-mountain, 0:00 to 6:00 local time for mountain
    117   is_mountain = dt['Mountain site'].loc[dt['ID'] == station].values[0]
    118   df['local time'] =  df['hour'] + shift
    119   if is_mountain == 'No':
    120     df = df.loc[ (df['local time'] < 16) & (df['local time'] >= 12)]
    121   else:
    122     df = df.loc[ (df['local time'] < 6) & (df['local time'] >= 0)]
    123   # network
    124   df = df.assign(network = 'VERIFYcore')
    125   # parameter
    126   df = df.assign(parameter = spec.upper())
    127   # obs
    128   df.rename(columns={'value':'obs'}, inplace = True)
    129   # filter out negative concentrations!
    130   df = df.loc[df['obs'] >= 0.]
    131   # obserror
    132   # filter out negative uncertainties
    133   df = df.loc[df['value_unc'] >=0. ]
    134   if station in TNO_error_station:
    135       print(station + ' is in TNO_error_station')
    136       if (len(df) == (df['value_unc'] >= 0.).sum()) and \
    137          (len(df) == (df['value_unc'] - df['obs'] <0.).sum()) and \
    138          (len(df) == (df['value_unc'] < max_unc).sum()):
    139              print('OK error')
    140              if ~np.isnan(TNO_LBC_error[TNO_error_station.index(station)]):
    141                  df['obserror'] = np.sqrt( \
    142                        (TNO_LBC_error[TNO_error_station.index(station)])**2 + \
    143                        (TNO_T_error[TNO_error_station.index(station)])**2 + \
    144                        (TNO_repr_error[TNO_error_station.index(station)])**2 +
    145                        df['value_unc']**2 \
    146                        )
    147              elif TNO_station_category[TNO_error_station.index(station)] == 'C':
    148                  df['obserror'] = np.sqrt( \
    149                        (np.nanmean(costal_seor_mean))**2 +
    150                        df['value_unc']**2 \
    151                        )
    152              elif TNO_station_category[TNO_error_station.index(station)] == 'M':
    153                  df['obserror'] = np.sqrt( \
    154                        (np.nanmean(montain_seor_mean))**2 +
    155                        df['value_unc']**2 \
    156                        )          
    157              elif TNO_station_category[TNO_error_station.index(station)] == 'O':
    158                  df['obserror'] = np.sqrt( \
    159                        (np.nanmean(other_seor_mean))**2 +
    160                        df['value_unc']**2 \
    161                        )
    162       else:
    163         print('not OK for error')
    164         sys.exit()
    165         
    166                  
    167   else:
    168       print(station + ' is NOT in TNO_error_station')
    169       sys.exit()
    170   # alt (m a.s.l)
    171   df.rename(columns={'altitude': 'alt'}, inplace = True)
    172
    173   # extend the output datastore
    174   column2save = ['date', 'duration', 'station', 'network', 'parameter', 'lon', 'lat', 'obs', 'obserror', 'alt']
    175   ds = ds.append(df.loc[:, column2save])
    176
    177# Dump the datastore
    178dump_datastore(ds, "{}/monitor_VERIFYcore_{}.nc".format(dir_data, spec), dump_default=False, col2dump=ds.columns)
    179
    180# Checking the consistency of the datastore when reading
    181ds2 = read_datastore("monitor.nc")
    182
    183
    
    1. one job per year is run, based on the following configuration file yml.sed:

    Show/Hide Code

      1#####################
      2# PYCIF config file #
      3#####################
      4
      5#####################################################################
      6#####################################################################
      7# PYCIF parameters
      8pycif_root: &pycif_root /home/users/ipison/cif/
      9workdir: &workdir /home/chimereicos/VERIFYoutputs/surface-inv3-_YYYY_
     10outdir: &outdir !join [*workdir, /outdir]
     11
     12# Verbose level
     13verbose: 1
     14
     15# Log file (to be saved in $wordkir)
     16#logfile: tltestlog
     17logfile: pyCIF_inv.logtest
     18
     19# Initial date
     20datei: _YYYY_-01-01 00:00:00
     21datef: _YYNN_-01-01 00:00:00
     22
     23#####################################################################
     24#####################################################################
     25# Running mode for PYCIF
     26mode:
     27  minimizer:
     28    maxiter: 20
     29    epsg: 0.1
     30    df1: 0.01
     31    plugin:
     32      name: M1QN3
     33      version: std
     34    simulator:
     35      plugin:
     36        name: gausscost
     37        version: std
     38      reload_from_previous: true
     39  plugin:
     40    name: 4dvar
     41    version: std
     42  save_out_netcdf: true
     43
     44#####################################################################
     45#####################################################################
     46# Observation operator
     47obsoperator:
     48  plugin:
     49    name: standard
     50    version: std
     51    type: obsoperator
     52  autorestart: True
     53#####################################################################
     54#####################################################################
     55controlvect:
     56  plugin:
     57    name: standard
     58    version: std
     59    type: controlvect
     60#  save_out_netcdf: True
     61  
     62#####################################################################
     63#####################################################################
     64# Transport model
     65model :
     66  plugin:
     67    name: CHIMERE
     68    version : std
     69    type: model
     70
     71  direxec: !join [*pycif_root, /model_sources/chimereGES/]
     72  nphour_ref: 6
     73  ichemstep: 2
     74  ideepconv: 0
     75  nivout: 17
     76  nlevemis: 1
     77
     78  nzdoms: 6
     79  nmdoms: 1
     80 
     81  force_clean_run: True
     82  useRAMonly: True 
     83  ignore_input_dates: True
     84
     85####################################################################
     86#####################################################################
     87obsvect:
     88  plugin:
     89    name: standard
     90    version: std
     91    type: obsvect
     92  dir_obsvect: !join [*outdir, /ref_obsvect]
     93  dump_type: nc
     94
     95#####################################################################
     96#####################################################################
     97platform:
     98  plugin:
     99    name: LSCE
    100    version: obelix
    101
    102#####################################################################
    103#####################################################################
    104# Domain definition
    105domain :
    106  plugin:
    107   name    : CHIMERE
    108   version : std
    109  repgrid : /home/chimereges/PYCIF_TEST_DATA/CHIMERE/domains
    110  domid : EUROCOMex3
    111  nlev: 17
    112  p1: 997
    113  pmax: 200
    114
    115#####################################################################
    116#####################################################################
    117# Chemical scheme
    118chemistry:
    119  plugin:
    120    name: CHIMERE
    121    version: gasJtab
    122  schemeid: ges.espigrad.humid
    123  dir_precomp : /home/users/ipison/CHIMERE/VERIFY/WP4_surface/
    124
    125#####################################################################
    126#####################################################################
    127datavect:
    128  plugin:
    129    name: standard
    130    version: std
    131  components:
    132    meteo:
    133      plugin:
    134        name: CHIMERE
    135        version: std      
    136        type: meteo
    137      dir: /home/satellites15/afortems/data_EUROCOM/%Y/
    138      file: METEO.%Y%m%d%H.24.nc
    139      file_freq: 1D
    140
    141    concs:
    142      parameters:
    143        CH4: 
    144          dir: /home/chimereicos/VERIFYdataWP4/Obs2021/core/
    145          file: monitor_VERIFYcore_ch4.nc
    146          vertical_interpolation:
    147            file_statlev: /home/users/ipison/CHIMERE/VERIFY/WP4_surface/stations_levels_CHIMERE17lev.txt
    148 
    149    flux:
    150      parameters:
    151        CH4:
    152          plugin:
    153            name: CHIMERE
    154            version: AEMISSIONS
    155            type: flux
    156          dir: /home/chimereicos/VERIFYdataWP4/Fluxes_intercomp_EU/
    157          file: AEMISSIONS.%Y%m0100.24.nc
    158          file_freq: 1D
    159          hresol : hpixels
    160          vresol : vpixels
    161          tresol: 24H
    162          nlev : 1
    163          err : 1.0 
    164          hcorrelations:
    165             landsea: True 
    166             sigma_land: 200 
    167             sigma_sea: 1000
    168             filelsm: /home/users/ipison/CHIMERE/VERIFY/intercomp_EU/lsmEUROCOMex3.nc
    169             dump_hcorr: True
    170          tcorrelations:
    171             sigma_t: 2500H
    172             dump_tcorr: True
    173    inicond:
    174      parameters:
    175        CH4: 
    176          plugin:
    177            name: CAMS
    178            version: netcdf
    179            type: field        
    180          dir: /home/comdata2/Fields/CAMSv19r1/
    181          comp_type: inicond
    182          file_freq: 1MS
    183          file:  cams73_v19r1_ch4_conc_surface_inst_%Y01.nc
    184          hresol: hpixels
    185          vresol: vpixels
    186          nlev: 34
    187          err: 0.1
    188          hcorrelations:
    189             sigma: 1000 
    190    latcond:
    191      parameters:
    192        CH4: 
    193          plugin:
    194            name: CAMS
    195            version: netcdf
    196            type: field        
    197          dir: /home/comdata2/Fields/CAMSv19r1/
    198          file:  cams73_v19r1_ch4_conc_surface_inst_%Y%m.nc
    199          comp_type: latcond
    200          file_freq: 1MS
    201          hresol: hpixels
    202          vresol: vpixels
    203          tresol : 2D
    204          nlev: 34 
    205          err: 0.1
    206          hcorrelations: 
    207             sigma: 1000 
    208          tcorrelations: 
    209             sigma_t: 120H
    210
    211    topcond:
    212      parameters:
    213        CH4: 
    214          plugin:
    215            name: CAMS
    216            version: netcdf
    217            type: field     
    218          dir: /home/comdata2/Fields/CAMSv19r1/
    219          file: cams73_v19r1_ch4_conc_surface_inst_%Y%m.nc
    220          comp_type: topcond
    221          file_freq: 1MS
    222          hresol: global
    223          tresol : 2D
    224          err: 0.1
    225          tcorrelations:
    226            sigma_t: 120H
    227        
    
    1. post-processing of the results include:

      1. maps of fluxes and increments:

      Show/Hide Code

        1#####################################################
        2# This code is provided by: I. Pison
        3# contact: isabelle.pison@lsce.ipsl.fr
        4# date: June 2022
        5# project: H2020 VERIFY
        6# relation to CIF: uses plugins, not embedded in.
        7# It is an example, to be used as a basis
        8# for developing other work.
        9# No guarantee as to its compatibility or results
       10# is attached since it is not part of the CIF
       11# Please don't forget acknowledgements if you use it.
       12#####################################################
       13
       14
       15import numpy as np
       16import pandas as pd
       17from pycif.utils.netcdf import readnc
       18from pycif.utils.datastores import dump
       19import matplotlib.pyplot as plt
       20import matplotlib.colors as colors
       21from matplotlib import ticker
       22import datetime
       23import xarray as xr
       24import cartopy.crs as crs
       25from mpl_toolkits.axes_grid1.inset_locator import InsetPosition
       26import cartopy
       27from cartopy.feature import ShapelyFeature
       28from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
       29
       30dir_results = '/home/chimereicos/VERIFYoutputs/surface-inv2-' 
       31dir_out = '/controlvect/flux/'
       32
       33year_beg = 2005
       34year_end = 2017
       35Na = 6.022e23
       36fs = 15
       37list_pr = []
       38list_po =[]
       39list_diff= []
       40
       41for yy in range(year_beg, year_end + 1, 1):
       42  print(yy)
       43  fcv = '{}{}{}/controlvect_flux_CH4.nc'.format(dir_results, yy, dir_out)
       44  flpr = readnc(fcv,['xb_phys'])
       45  flpo = readnc(fcv,['x_phys'])
       46  list_pr.append(flpr.mean(axis=0)[0,:,:]/ Na * 16.e6) # molec -> microgCH4
       47  list_po.append(flpo.mean(axis=0)[0,:,:]/ Na * 16.e6) # molec -> microgCH4
       48  diff = (flpo - flpr)/flpr * 100.
       49  list_diff.append(diff.mean(axis=0)[0,:,:])
       50
       51clon = readnc(fcv,['longitudes_corner'])
       52clat = readnc(fcv,['latitudes_corner'])
       53
       54############
       55# To plot: #
       56############
       57rivers = cartopy.feature.NaturalEarthFeature(
       58    category='physical', name='rivers_lake_centerlines',
       59    scale='10m', facecolor='none', edgecolor='blue')
       60borders = cartopy.feature.NaturalEarthFeature(
       61    category='cultural', name='admin_0_boundary_lines_land',
       62    scale='10m', facecolor='none', edgecolor='black')
       63
       64#############
       65# plot maps #
       66#############
       67# color scale with Ncolor colors from cmap base colormap for fluxes
       68Ncolor = 15
       69cmap = plt.cm.jet
       70cmaplist = cmap(np.linspace(0, 1, Ncolor))
       71cmap_plot = cmap.from_list('custom', cmaplist, Ncolor)
       72
       73# same for increments
       74Ncolor2 = 11
       75cmap2 = plt.cm.bwr
       76cmaplist2 = cmap2(np.linspace(0, 1, Ncolor2))
       77cmap2 = cmap2.from_list('custom', cmaplist2, Ncolor2)
       78
       79lognorm = colors.LogNorm()
       80
       81##plot a given year
       82yy = 2010
       83nby = yy - year_beg
       84fig = plt.figure(figsize=(21, 10))
       85gs = fig.add_gridspec(nrows = 2, ncols = 4, width_ratios = [2,2,2,0.1], 
       86                 left=0., right=0.97, bottom=0, top=1,
       87                  height_ratios = [1, 1 ] ,hspace = 0.01, wspace = 0.25)
       88subprior = fig.add_subplot(gs[0,1], projection=crs.PlateCarree())
       89subprior.pcolor(clon, clat, list_pr[nby],
       90                         norm=lognorm,
       91                         cmap = cmap,
       92                         transform=crs.PlateCarree())
       93subprior.coastlines(resolution='50m')
       94subprior.add_feature(borders,linewidth=0.5,linestyle='--')
       95subprior.set_title("Prior", fontsize = fs)
       96
       97subpost = fig.add_subplot(gs[0,2], projection=crs.PlateCarree())
       98im = subpost.pcolor(clon, clat, list_po[nby],
       99                         norm=lognorm,
      100                         cmap = cmap,
      101                         transform=crs.PlateCarree())
      102subpost.coastlines(resolution='50m')
      103subpost.add_feature(borders,linewidth=0.5,linestyle='--')
      104subpost.set_title("Posterior", fontsize = f)
      105
      106echelle = fig.add_subplot(gs[0,3])
      107ip = InsetPosition(echelle, [0,0,1,1]) 
      108echelle.set_axes_locator(ip)
      109p00 = echelle.get_position()
      110p01 = echelle.get_position()
      111p02 = subpost.get_position()
      112p00_new = [p00.x0, p02.y0, p00.width, p02.height]
      113echelle.set_position(p00_new)
      114formatter = ticker.LogFormatter(10, labelOnlyBase=False) 
      115cbar = plt.colorbar(im, format=formatter, cax=echelle)
      116cbar.ax.set_title('molec.cm$^{-2}$.s$^{-1}$', fontsize = fs) 
      117
      118subd = fig.add_subplot(gs[1,2], projection=crs.PlateCarree())
      119im2 = subd.pcolor(clon, clat, list_diff[nby],
      120                           vmin=-100, vmax=100,
      121                           cmap = cmap2,
      122                           transform=crs.PlateCarree())
      123subd.coastlines(resolution='50m')
      124subd.add_feature(borders,linewidth=0.5,linestyle='--')
      125subd.set_title("Increments: posterior - prior (%)", fontsize = fs)
      126echelle2 = fig.add_subplot(gs[1,3])
      127ip = InsetPosition(echelle2, [0,0,1,1]) 
      128echelle2.set_axes_locator(ip)
      129p00 = echelle2.get_position()
      130p01 = echelle2.get_position()
      131p02 = subd.get_position()
      132p00_new = [p00.x0, p02.y0, p00.width, p02.height]
      133echelle2.set_position(p00_new)
      134cbar2 = plt.colorbar(im2, cax=echelle2, ax = subd)
      135cbar2.ax.set_title('%', fontsize = fs) 
      136
      137plt.savefig('flux{}.png'.format(yy))
      138
      139
      140
      
      1. time series of mixing ratios at the stations (obs, prior, post) and maps of statistical indicators at the stations (ratios posterior to prior):

      Show/Hide Code

        1#####################################################
        2# This code is provided by: I. Pison
        3# contact: isabelle.pison@lsce.ipsl.fr
        4# date: June 2022
        5# project: H2020 VERIFY
        6# relation to CIF: uses plugins, not embedded in.
        7# It is an example, to be used as a basis
        8# for developing other work.
        9# No guarantee as to its compatibility or results
       10# is attached since it is not part of the CIF
       11# Please don't forget acknowledgements if you use it.
       12#####################################################
       13
       14import numpy as np
       15import pandas as pd
       16from pycif.utils.netcdf import readnc
       17from pycif.utils.datastores import dump
       18import matplotlib.pyplot as plt
       19import datetime
       20import cartopy.crs
       21from cartopy.feature import ShapelyFeature
       22from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
       23from mpl_toolkits.axes_grid1.inset_locator import InsetPosition
       24
       25dir_results = '/home/chimereicos/VERIFYoutputs/surface-inv2-' 
       26dir_out = '/obsoperator/fwd_'
       27dir_monit = '/obsvect/concs/CH4/'
       28
       29year_beg = 2005
       30year_end = 2017
       31
       32fs = 30 # fontsize
       33
       34# read monitors of all years
       35ds_all_years = pd.DataFrame()
       36for yy in range(year_beg, year_end + 1, 1):
       37  print(yy)
       38  fmonitpr = '{}{}{}-001{}/monitor.nc'.format(dir_results, yy, dir_out, dir_monit)
       39  fmonitpo = '{}{}{}-002{}/monitor.nc'.format(dir_results, yy, dir_out, dir_monit) 
       40  dspr = dump.read_datastore(fmonitpr)
       41  dspo = dump.read_datastore(fmonitpo)
       42  dsprpo = pd.concat( [dspr, dspo['maindata']['sim']], axis = 1)
       43  ds_all_years = pd.concat( [dsprpo, ds_all_years], axis = 0)
       44
       45list_stations = list(ds_all_years[('metadata', 'station')].unique())
       46list_lat = []
       47list_lon = []
       48lat_etiq = []
       49lon_etiq = []
       50for stat in list_stations:
       51  lat = ds_all_years[('metadata', 'lat')].loc[ds_all_years[('metadata', 'station')] == stat].iloc[0]
       52  lon = ds_all_years[('metadata', 'lon')].loc[ds_all_years[('metadata', 'station')] == stat].iloc[0]
       53  list_lat.append(lat)
       54  list_lon.append(lon)
       55  if stat =='ZSF':
       56    lat_etiq.append(lat - 0.3)
       57    lon_etiq.append(lon + 0.2) 
       58  elif stat =='KIT':
       59    lat_etiq.append(lat - 0.5)
       60    lon_etiq.append(lon + 0.2) 
       61  else:
       62    lat_etiq.append(lat + 0.2)
       63    lon_etiq.append(lon + 0.2)
       64    
       65dstat = pd.DataFrame(columns = ['station','Prior bias', 'Posterior bias', 'Prior rms', 'Posterior rms', 'Prior corr', 'Posterior corr'])
       66dstat['station'] = list_stations
       67for stat in list_stations:
       68  print(stat)
       69  ds = ds_all_years.loc[ ds_all_years[('metadata', 'station')] == stat ]
       70  ds['var'] = ds[('maindata','obserror')].pow(2)
       71  # compute bias, standard deviation, time correlation
       72  ds['prior difference'] = ds[('maindata', 'obs')] - ds[('maindata', 'sim')]
       73  ds['post difference'] = ds[('maindata', 'obs')] - ds['sim']
       74  dstat['Prior bias'].loc[dstat['station'] == stat ] = ds['prior difference'].mean()
       75  dstat['Posterior bias'].loc[dstat['station'] == stat ] = ds['post difference'].mean()
       76  dstat['Prior rms'].loc[dstat['station'] == stat ] =  np.sqrt(ds['prior difference'].pow(2).mean())
       77  dstat['Posterior rms'].loc[dstat['station'] == stat ] =  np.sqrt(ds['post difference'].pow(2).mean())
       78  dstat['Prior corr'].loc[dstat['station'] == stat ] =  ds[('maindata', 'obs')].corr(ds[('maindata', 'sim')])
       79  dstat['Posterior corr'].loc[dstat['station'] == stat] = ds[('maindata', 'obs')].corr(ds['sim'])
       80  # monthly means
       81  ds2 = ds.resample('MS', on = ('metadata', 'date'), closed = 'left')
       82  dsm = ds2.mean()
       83  counts = ds2.size()
       84  dsm['meanobserror'] =  dsm['var'].pow(0.5)/counts.pow(0.5)
       85# plot times series at each station
       86  plt.figure(figsize=(21,11))
       87  ax = plt.gca()
       88  ax.tick_params(axis = 'both', which = 'major', labelsize = fs/2)
       89  ax.set_xticks([datetime.datetime(i,1,1) for i in range(year_beg, year_end + 1)])
       90  plt.title('Methane mixing ratios at {}'.format(stat), fontsize = fs)
       91  plt.xlabel('Date', size = fs)
       92  plt.ylabel('[CH4] (ppb)', size = fs)
       93  plt.errorbar(dsm.index, dsm[('maindata','obs')], \
       94               dsm['meanobserror'], label = 'measurements', \
       95               linestyle = '-', marker = 'o', ms = 1)
       96  plt.plot(dsm.index, dsm[('maindata','sim')], \
       97               label = 'prior', \
       98               linestyle = '-', marker = 'o', ms = 1)
       99  plt.plot(dsm.index, dsm['sim'], \
      100               label = 'analysis', \
      101               linestyle = '-', marker = 'o', ms = 1)
      102  plt.xlim(datetime.datetime(2005,1,1,0),datetime.datetime(2018,1,1,0))  
      103  plt.legend(fontsize = fs/2)
      104  plt.savefig('ts_{}.png'.format(stat))
      105  plt.close()  
      106
      107# maps for statistics
      108Ncolor = 10
      109cmap = plt.cm.bwr
      110cmaplist = cmap(np.linspace(0, 1, Ncolor))
      111cmap_plot = cmap.from_list('custom', cmaplist, Ncolor)
      112cmin, cmax = 0, 2
      113rivers = cartopy.feature.NaturalEarthFeature(
      114    category='physical', name='rivers_lake_centerlines',
      115    scale='10m', facecolor='none', edgecolor='blue')
      116borders = cartopy.feature.NaturalEarthFeature(
      117    category='cultural', name='admin_0_boundary_lines_land',
      118    scale='10m', facecolor='none', edgecolor='black')
      119latmin=33.
      120latmax=73.
      121longmin=-15.
      122longmax=35.
      123fig = plt.figure(figsize=(21, 10))
      124gs = fig.add_gridspec(nrows = 2, ncols = 3, height_ratios = [1, 0.02 ] )
      125s0 = fig.add_subplot(gs[0, 0], projection=cartopy.crs.PlateCarree())
      126s0.scatter(list_lon, list_lat, c=np.abs(dstat['Posterior bias']/dstat['Prior bias']), 
      127                    s=20, cmap = cmap, marker = 's', vmin = 0, vmax = 2, edgecolors = 'k',
      128                    linewidths = 0.4)
      129s0.set_title('Ratio of posterior to prior mean difference\n between simulation and measurements')
      130s0.set_extent([longmin, longmax, latmin, latmax])
      131s0.coastlines(linewidth=0.5, resolution='50m')
      132s0.add_feature(borders, linewidth=0.5, linestyle='--')
      133
      134s1 = fig.add_subplot(gs[0, 1], projection=cartopy.crs.PlateCarree())
      135im = s1.scatter(list_lon, list_lat, c=dstat['Posterior rms']/dstat['Prior rms'], 
      136                    s=20, cmap = cmap, marker = 's',vmin = 0, vmax = 2, edgecolors = 'k',
      137                    linewidths = 0.4)
      138s1.set_title('Ratio of posterior to prior quadratic mean\n of the differences\n between simulation and measurements')
      139s1.set_extent([longmin, longmax, latmin, latmax])
      140s1.coastlines(linewidth=0.5, resolution='50m')
      141s1.add_feature(borders, linewidth=0.5, linestyle='--')
      142
      143s2 = fig.add_subplot(gs[0, 2], projection=cartopy.crs.PlateCarree())
      144s2.scatter(list_lon, list_lat, c=np.abs(dstat['Posterior corr'])/np.abs(dstat['Prior corr']), 
      145                    s=20, cmap = cmap, marker = 's', vmin = 0, vmax = 2, edgecolors = 'k',
      146                    linewidths = 0.4)
      147s2.set_title('Ratio of posterior to prior absolute values\n  of time correlation\n between simulation and measurements')
      148s2.set_extent([longmin, longmax, latmin, latmax])
      149s2.coastlines(linewidth=0.5, resolution='50m')
      150s2.add_feature(borders, linewidth=0.5, linestyle='--')
      151
      152cb = fig.add_subplot(gs[1,1])
      153cbar = plt.colorbar(im, cax = cb, ax = s0 , orientation = 'horizontal')
      154cbar.ax.tick_params(labelsize=18)
      155plt.savefig('map_statistics.png')
      156
      157
      
      1. yearly emission budgets for various regions (EU27+UK, EU27 countries, groups of countries):

      Show/Hide Code

        1#####################################################
        2# This code is provided by: I. Pison
        3# contact: isabelle.pison@lsce.ipsl.fr
        4# date: June 2022
        5# project: H2020 VERIFY
        6# relation to CIF: uses plugins, not embedded in.
        7# It is an example, to be used as a basis
        8# for developing other work.
        9# No guarantee as to its compatibility or results
       10# is attached since it is not part of the CIF
       11# Please don't forget acknowledgements if you use it.
       12#####################################################
       13
       14import numpy as np
       15import pandas as pd
       16from pycif.utils.netcdf import readnc
       17from pycif.utils.datastores import dump
       18import matplotlib.pyplot as plt
       19import matplotlib.colors as colors
       20import datetime
       21import xarray as xr
       22import cartopy.crs as crs
       23import netCDF4 as nc
       24from mpl_toolkits.axes_grid1.inset_locator import InsetPosition
       25import calendar
       26
       27# read masks for various regions
       28fmasks = '/home/chimereicos/VERIFYoutputs/masks_for_EU27+UK.nc'
       29fm = nc.Dataset(fmasks)
       30nom_reg = 'EU_27+UK'
       31reg = np.array(fm[nom_reg])
       32mask = (reg > 0 )
       33
       34dir_results = '/home/chimereicos/VERIFYoutputs/surface-inv2-' 
       35dir_out = '/controlvect/flux/'
       36
       37year_beg = 2005
       38year_end = 2017
       39dy = 0.5
       40dx = 0.5
       41rearth=637103000. 
       42Na = 6.022e23
       43fcv = '{}{}{}/controlvect_flux_CH4.nc'.format(dir_results, year_beg, dir_out)
       44zlon = readnc(fcv,['longitudes'])[0,:]
       45zlat = readnc(fcv,['latitudes'])[:,0]
       46nlat = len(zlat)
       47nlon = len(zlon)
       48area = np.zeros((nlat,nlon))
       49for i in range(nlat):
       50  for j in range(nlon):
       51     lat1 = zlat[i]
       52     lat2 =  zlat[i] + dy
       53     # formula with lat between -90:90
       54     area[i,j] = rearth * rearth * dx * np.pi/180. * \
       55                 (np.cos(lat1 * np.pi/180. + np.pi/2.) - \
       56                  np.cos(lat2 * np.pi/180. + np.pi/2.))
       57                  
       58yrflregpr = []
       59yrflregpo = []
       60summ = np.zeros((2,12))
       61# read controlvect
       62for yy in range(year_beg, year_end + 1, 1):
       63  print(yy)
       64  fcv = '{}{}{}/controlvect_flux_CH4.nc'.format(dir_results, yy, dir_out)
       65  ds = xr.open_dataset(fcv)
       66  flpr = readnc(fcv,['xb_phys'])
       67  flpo = readnc(fcv,['x_phys'])
       68  flpran = np.sum(np.sum(flpr,axis=1),axis=0) * area * 3600. / Na * 16.e-3
       69  flpoan = np.sum(np.sum(flpo,axis=1),axis=0) * area * 3600. / Na * 16.e-3
       70  flreg = np.sum(flpran * mask) * 1e-3 * 1e-3 # kgCH4.yr -> t/an -> kt/an
       71  yrflregpr.append(flreg)
       72  flreg = np.sum(flpoan * mask) * 1e-3 * 1e-3 # kgCH4.yr -> t/an -> kt/an
       73  yrflregpo.append(flreg)
       74  beg = 0
       75  for m in range(12):
       76    mm = m + 1
       77    nbd = calendar.mdays[mm] + ( mm ==2 and calendar.isleap(yy))
       78    flprmm = np.sum(np.sum(flpr,axis=1)[beg:beg + nbd * 24],axis=0) * area * 3600. / Na * 16.e-3
       79    flpomm = np.sum(np.sum(flpo,axis=1)[beg:beg + nbd * 24],axis=0) * area * 3600. / Na * 16.e-3
       80    beg = beg + nbd * 24
       81    summ[0, m] = summ[0, m] + np.sum(flprmm * mask) * 1e-3 * 1e-3
       82    summ[1, m] = summ[1, m] + np.sum(flpomm * mask) * 1e-3 * 1e-3
       83
       84####################
       85# plot time series #
       86####################
       87fs = 25
       88plt.figure(figsize=(21,11))
       89plt.title('Methane emissions in '+nom_reg, fontsize=fs)
       90plt.xlabel('Year', fontsize=fs)
       91plt.ylabel('ktCH$_4$', fontsize=fs)
       92list_ticks = [ i  for i in range(year_end - year_beg + 1)]
       93list_names = [ year_beg + i for i in range(year_end - year_beg + 1)]
       94plt.xticks(list_ticks, list_names, fontsize=15)
       95plt.tick_params(axis='y', labelsize=30)
       96plt.yticks(fontsize=fs)
       97plt.bar(list_ticks, yrflregpr, label = 'prior', width=-0.4, align='edge')
       98plt.bar(list_ticks, yrflregpo, label = 'posterior', width=0.4, align='edge')
       99plt.legend(fontsize=fs)
      100plt.savefig('ts_fluxes_{}_yearly.png'.format(nom_reg))
      101plt.close()
      102
      103plt.figure(figsize=(21,11))
      104plt.title('Methane emissions: average seasonal cycle over '+str(year_beg)+'-'+str(year_end)+' in '+nom_reg, fontsize=fs)
      105plt.xlabel('Month', fontsize=fs)
      106plt.ylabel('ktCH$_4$', fontsize=fs)
      107ax = plt.gca()
      108ax.tick_params(axis = 'both', which = 'major', labelsize = fs/2)
      109plt.xticks([i+1 for i in range(12)] , [str(i+1) for i in range(12) ])
      110nbyears = year_end - year_beg + 1
      111plt.plot([i+1 for i in range(12)], summ[0,:]/nbyears, label = 'prior', \
      112               linestyle = '--', marker = 'o', ms = 8)
      113plt.plot([i+1 for i in range(12)], summ[1,:]/nbyears, label = 'posterior', \
      114              linestyle = '-', marker = 'o', ms = 8)
      115plt.legend(fontsize=fs-5)
      116plt.savefig('ts_seasonal_{}_mass.png'.format(nom_reg))
      117
      
      1. time series of emissions for various regions (EU27+UK, EU27 countries, groups of countries):

      Show/Hide Code

       1#####################################################
       2# This code is provided by: I. Pison
       3# contact: isabelle.pison@lsce.ipsl.fr
       4# date: June 2022
       5# project: H2020 VERIFY
       6# relation to CIF: uses plugins, not embedded in.
       7# It is an example, to be used as a basis
       8# for developing other work.
       9# No guarantee as to its compatibility or results
      10# is attached since it is not part of the CIF
      11# Please don't forget acknowledgements if you use it.
      12#####################################################
      13
      14import numpy as np
      15import pandas as pd
      16from pycif.utils.netcdf import readnc
      17from pycif.utils.datastores import dump
      18import matplotlib.pyplot as plt
      19import matplotlib.colors as colors
      20import datetime
      21import xarray as xr
      22import cartopy.crs as crs
      23import netCDF4 as nc
      24from mpl_toolkits.axes_grid1.inset_locator import InsetPosition
      25
      26# read masks for various regions
      27fmasks = '/home/chimereicos/VERIFYoutputs/masks_for_EU27.nc'
      28fm = nc.Dataset(fmasks)
      29reg_name = 'EU_27'
      30reg = np.array(fm[reg_name])
      31
      32dir_results = '/home/chimereicos/VERIFYoutputs/surface-inv-' 
      33dir_out = '/controlvect/flux/'
      34
      35year_beg = 2005
      36year_end = 2017
      37
      38# read controlvect
      39list_dates = []
      40tsprEU = []
      41tspoEU = []
      42for yy in range(year_beg, year_end + 1, 1):
      43  print(yy)
      44  fcv = '{}{}{}/controlvect_flux_CH4.nc'.format(dir_results, yy, dir_out)
      45  ds = xr.open_dataset(fcv)
      46  dates = ds['time_phys'].to_dataframe()
      47  list_dates.append(dates['time_phys'].loc[(dates.index.day == 1) & (dates.index.hour == 1)].tolist())
      48  flpr = readnc(fcv,['xb_phys'])
      49  flpo = readnc(fcv,['x_phys'])
      50  # over the selected region
      51  flprEU = flpr[:-1,:,:,:] * reg
      52  tsprEU.append(np.mean(flprEU))
      53  flpoEU = flpo[:-1,:,:,:] * reg
      54  tspoEU.append(np.mean(flpoEU))
      55   
      56####################
      57# plot time series #
      58####################
      59plt.figure(figsize=(21,11))
      60plt.title('Methane emissions: average in {}'.format(reg_name), fontsize=15)
      61plt.xlabel('Date', fontsize=15)
      62plt.ylabel('molec.cm$^{-2}$.s$^{-1}$', fontsize=15)
      63list_ticks = []
      64for ld in list_dates:
      65    list_ticks.append(ld[0])
      66list_names = []
      67for d in list_ticks:
      68 list_names.append(d.strftime('%Y'))
      69plt.xticks(list_ticks, list_names, fontsize=15)
      70plt.tick_params(axis='y', labelsize=30)
      71plt.yticks(fontsize=15)
      72plt.bar(list_ticks, tsprEU, label = 'prior', width=-100, align='edge')
      73plt.bar(list_tciks, tspoEU, label = 'posterior', width=100, align='edge')
      74plt.legend(fontsize=15)
      75plt.savefig('ts_fluxes_{}.png'.format(reg_name))
      76plt.close()
      77
      
  • inversions with TROPOMI data:

    1. plot raw TROPOMI data:

      Show/Hide Code

        1###################################################################
        2# This code is provided by: I. Pison and R. Plauchu
        3# contact: isabelle.pison@lsce.ipsl.fr, robin.plauchu@lsce.ipsl.fr
        4# date: June 2022
        5# project: H2020 VERIFY, ANR ARGONAUT, CNES Tosca ARGOS
        6# relation to CIF: none
        7# It is an example, to be used as a basis
        8# for developing other work.
        9# No guarantee as to its compatibility or results
       10# is attached since it is not part of the CIF
       11# Please don't forget acknowledgements if you use it.
       12####################################################################
       13
       14import calendar
       15import os
       16import datetime
       17import sys
       18import netCDF4 as nc
       19import numpy as np
       20import matplotlib.pyplot as plt
       21import matplotlib.axes as maxes
       22from mpl_toolkits.axes_grid1 import make_axes_locatable
       23import cartopy.crs as crs
       24import cartopy.feature as cfeature
       25from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
       26from matplotlib.collections import PolyCollection
       27
       28dir_sat = '/home/chimereges2/ipison/testdownload/'
       29spec = 'CH4'
       30domain = 'EUROCOM'
       31
       32# selec time slot when the satellite goes over the targeted area
       33hh_min = '14'
       34hh_max = '10'
       35
       36# Projection
       37projection = crs.PlateCarree()
       38
       39# Satellite data
       40mission_name       = 'S5P'
       41processing_stream  = '*'
       42product_identifier = 'L2__{}____'.format(spec)
       43
       44# Time period
       45year   = 2018
       46yearn = 2021
       47month  = 4
       48monthn = 2
       49day    = 30
       50dayn   = 15
       51
       52# from readme: Sentinel-5P-Methane-Product-Readme-File.pdf, dated 202-12-11
       53Qa = 0.5 # strictly superior to 0.5
       54
       55# For the figures:
       56# Layout parameters
       57subplot_title_pad = 15
       58subplot_title_fontsize  = 6
       59gridlines_fontsize      = 5
       60colorbar_title_fontsize = 5
       61colorbar_tick_fontsize  = 5
       62colorbar_shrink=0.8
       63colorbar_pad=0.1
       64# Font
       65nice_fonts = {
       66        #"text.usetex": True,
       67        "font.family": "sans-serif", # Font with or without "empattement"
       68        "font.serif" : "Helvetica",
       69        }
       70plt.rcParams.update(nice_fonts)
       71# color map
       72Ncolor = 6
       73cmap  = plt.cm.jet
       74cmaplist = cmap(np.linspace(0, 1, Ncolor))
       75cmap = cmap.from_list('custom', cmaplist, Ncolor)
       76#-----------------------------------------------------------------------
       77# Use CHIMERE's domain files to set the extent of the maps to plot
       78# Also possible: provide hard-coded lon_min,lon_max,lat_min,lat_max
       79print('---------------------------------------------------------------')
       80print("Start read COORD")
       81dirdom = '/home/users/afortems/PYTHON/PYVAR_CHIMERE_svn_Europe_ok/reposcode/reposcode/project_example/domains/'
       82coord_file        = '{}/HCOORD/COORD_{}'.format(dirdom,domain)
       83coord_corner_file = '{}/HCOORD//COORDcorner_{}'.format(dirdom,domain)
       84# Domain size
       85listdom = '{}/domainlist.nml'.format(dirdom)
       86f = open(listdom, 'r')
       87availabledom = f.readlines()
       88for line in availabledom:
       89    tmp = line.split() #str.split(line)
       90    if tmp[0] == domain:
       91        nlon_target = int(tmp[1])
       92        nlat_target = int(tmp[2])
       93print('nlon, nlat :',nlon_target,nlat_target)
       94# Reading coord files
       95f1 = open(coord_corner_file, 'r')
       96coord1 = f1.readlines()
       97f1.close()
       98# Putting coords into matrix of coord corners
       99zlonc_target = np.zeros((nlon_target + 1, nlat_target + 1), float)
      100zlatc_target = np.zeros((nlon_target + 1, nlat_target + 1), float)
      101k = 0
      102for i in range(nlat_target + 1):
      103    for j in range(nlon_target + 1):
      104        tmp = coord1[k].split() #str.split(coord1[k])
      105        k += 1
      106        zlonc_target[j, i] = float(tmp[0])
      107        zlatc_target[j, i] = float(tmp[1])
      108lon_min, lon_max = zlonc_target.min(), zlonc_target.max()
      109lat_min, lat_max = zlatc_target.min(), zlatc_target.max()
      110print('Mask :',lon_min,lon_max,lat_min,lat_max)
      111print("End read COORD")
      112
      113
      114print('---------------------------------------------------------------')
      115sat_data_id = mission_name+'_'+processing_stream+'_'+product_identifier
      116fic_out_part_name = 'TROPOMI_{}_{}'.format(spec,domain)
      117file_list_name    = 'filelist_{}_'.format(spec)
      118# Time - Period covered
      119for yy in range(year, yearn+1):
      120  print(yy)
      121  endm = 12
      122  begm = 1
      123  if (yy == yearn): endm = monthn
      124  if (yy == year) : begm = month
      125  for mm in range(begm, endm+1):
      126    print(mm)
      127    endday = ndays = calendar.mdays[mm] + ((mm == 2) and calendar.isleap(yy))
      128    begday = 1
      129    if (mm == monthn): endday = dayn
      130    if( mm == month ) : begday = day
      131    yyyymm = str(yy)+"%2.2i"%mm
      132    savefig_name_monthly = dir_sat + fic_out_part_name + '_pixel_Qa' + str(Qa) + '_' + yyyymm + '.png'
      133    list_vrts_m = []
      134    list_z_m = []
      135    nb_obs_trop_m = 0
      136    for dd in range(begday,endday+1):
      137      print( mm, dd)
      138      yyyymmdd = str(yy)+"%2.2i"%mm+"%2.2i"%dd
      139      savefig_name = dir_sat + fic_out_part_name + '_pixel_Qa' + str(Qa) + '_' + yyyymmdd + '.png'
      140   
      141      # List of netcdf files for one day included in [hh_min; hh_max] as
      142      # an aera may be covered on several scanlines while the satellite is orbiting
      143      datetime0 = datetime.datetime(year=yy, month=mm, day=dd, hour=int(hh_min), minute=0)
      144      os.system('ls ' + dir_sat + sat_data_id + datetime0.strftime("%Y%m%dT") + '{' + hh_min + '..' + hh_max + '}*' +\
      145         datetime0.strftime("%Y%m%dT") + '*.nc > ' + dir_sat + file_list_name + str(yy)+"%2.2i"%mm+"%2.2i"%dd+'.txt')
      146      ff = open(dir_sat + file_list_name + str(yy)+"%2.2i"%mm+"%2.2i"%dd+'.txt','r')
      147      sat_files = ff.readlines()
      148      ff.close()
      149      print('#-----------------------------------------------------------------------')
      150      for theline in sat_files:
      151        print(theline)
      152      print('#-----------------------------------------------------------------------')
      153      if len(sat_files) > 1:
      154        print('More than 1 orbit is selected')
      155      if len(sat_files) == 0: continue
      156      #-----------------------------------------------------------------------
      157      # Read each netcdf file and select variables of interest
      158      count = 0
      159      list_vrts = []
      160      list_z = []
      161      list_time_coverage_start = []
      162      list_orbit = []
      163      nb_obs_trop = 0
      164      for theline in sat_files:
      165        count += 1
      166        print('#-----------------------------------------------------------------------')
      167        print('File n°:',count,'/',len(sat_files))
      168        print(str.split(theline,'/')[-1])
      169
      170        s = str.split(theline)
      171
      172        thefile = nc.Dataset(str.strip(theline), 'r')
      173        # read groups and such in netcdf4 file
      174        product_grp = thefile.groups['PRODUCT'] 
      175        support_data_grp = product_grp.groups['SUPPORT_DATA']
      176        geoloc_grp  = support_data_grp.groups['GEOLOCATIONS']
      177        # selection on Qa XXet dans domaine EUROCOMXX
      178        mask = (product_grp.variables['qa_value'][:] > Qa) & \
      179              (product_grp.variables['longitude'][:] < lon_max) & \
      180              (product_grp.variables['longitude'][:] > lon_min) & \
      181              (product_grp.variables['latitude' ][:] < lat_max) & \
      182              (product_grp.variables['latitude' ][:] > lat_min)
      183        CH4  = product_grp.variables['methane_mixing_ratio'][:,:,:][mask]
      184        LAT  = product_grp.variables['latitude'][:,:,:][mask]
      185        LON  = product_grp.variables['longitude'][:,:,:][mask]
      186        LATC = geoloc_grp.variables['latitude_bounds'][:,:,:,:][mask]
      187        LONC = geoloc_grp.variables['longitude_bounds'][:,:,:,:][mask]
      188        nb_obs_loc = CH4.shape[0]
      189        print('Nb of obs :',CH4.shape)
      190        if(nb_obs_loc == 0): 
      191          continue
      192        nb_obs_trop = nb_obs_trop + nb_obs_loc
      193        nb_obs_trop_m = nb_obs_trop_m + nb_obs_loc
      194        # Make vertices for poly-collection
      195        list_vrts.append(np.array([LONC.T,LATC.T]).T)
      196        list_z.append(CH4)
      197        list_time_coverage_start.append(thefile.time_coverage_start)
      198        list_orbit.append(thefile.orbit)
      199        list_vrts_m.append(np.array([LONC.T,LATC.T]).T)
      200        list_z_m.append(CH4)
      201        thefile.close()
      202      if (nb_obs_trop == 0 ): continue
      203      ## figure "template"
      204      fig = plt.figure(figsize=(14, 5))
      205      gs = fig.add_gridspec(nrows = 1, ncols = 1, figure=fig,
      206                             width_ratios  = [1],
      207                             height_ratios = [1],
      208                             wspace = 0.2,
      209                             hspace = 0.2)
      210      ax0 = fig.add_subplot(gs[0,0],projection=projection)
      211      ax0.set_global()
      212      ax0.set_title('TROPOMI-S5P L2 {} Qa={} on {} - time: {}, orbit no'.format(spec,Qa,domain,list_time_coverage_start[:],list_orbit[:]),
      213                     pad      = subplot_title_pad,
      214                     fontsize = subplot_title_fontsize)
      215      ax0.text(0.95, 0.008, 'Nb of pix: {}'.format(nb_obs_trop),
      216                verticalalignment='bottom', horizontalalignment='right',
      217                transform=ax0.transAxes,
      218                color='red', fontsize=5)
      219      ax0.set_extent([lon_min-1, lon_max+1,
      220                       lat_min-1, lat_max+1], crs=projection)
      221      ax0.coastlines(color='black',    resolution='50m',  linewidth=0.3)
      222      ax0.add_feature(cfeature.BORDERS.with_scale('50m'), linewidth=0.2)
      223      gl_ax0 = ax0.gridlines(
      224                              crs=projection,
      225                              draw_labels=True,
      226                              linewidth=0.4, color='gray', alpha=0.5, linestyle='--')
      227      # Make label disappear top, bottom, left, right
      228      gl_ax0.xlabels_bottom = False
      229      gl_ax0.ylabels_right  = False
      230      # Change format type
      231      gl_ax0.xformatter = LONGITUDE_FORMATTER
      232      gl_ax0.yformatter = LATITUDE_FORMATTER
      233      # Monitor font options
      234      gl_ax0.xlabel_style = {'size': 5, 'color': 'gray'}
      235      gl_ax0.ylabel_style = {'size': 5, 'color': 'gray'}
      236
      237      ## color-in the figure
      238      # Display each satellite pixel as individual polygones
      239      for i, toplot in enumerate(list_vrts):
      240        collection = PolyCollection(toplot, array=list_z[i], cmap=cmap)
      241        ax0.add_collection(collection)
      242      # Manage colorbar
      243      divider0 = make_axes_locatable(ax0)
      244      cax0 = divider0.append_axes("right", size="3%", pad=colorbar_pad, axes_class=maxes.Axes)
      245      cb_ax0 = fig.colorbar(collection,
      246                             cax=cax0,
      247                             ax=ax0,
      248                             orientation='vertical',
      249                             extend='both',)
      250      cb_ax0.set_label(     "mixing ratio (ppb)",
      251                             fontsize=colorbar_title_fontsize)
      252      cb_ax0.ax.tick_params(direction='in',
      253                             labelsize=colorbar_tick_fontsize)
      254      cb_ax0.formatter.set_powerlimits((0, 0))
      255      cb_ax0.ax.yaxis.set_offset_position('left')
      256      cb_ax0.update_ticks()
      257
      258      plt.savefig(savefig_name, bbox_inches ="tight", orientation ='landscape', dpi=240)
      259      print('File created : ', savefig_name)
      260      print('#-----------------------------------------------------------------------')
      261   
      262
      263
      264
      265
      266    if (nb_obs_trop_m == 0 ): continue
      267    ## figure "template"
      268    figm = plt.figure(figsize=(14, 5))
      269    gsm = figm.add_gridspec(nrows = 1, ncols = 1, figure=figm,
      270                             width_ratios  = [1],
      271                             height_ratios = [1],
      272                             wspace = 0.2,
      273                             hspace = 0.2)
      274    ax0m = figm.add_subplot(gsm[0,0],projection=projection)
      275    ax0m.set_global()
      276    ax0m.set_title('TROPOMI-S5P L2 {} Qa={} on {} - year {} month {}'.format(spec,Qa,domain,yy,mm),
      277                     pad      = subplot_title_pad,
      278                     fontsize = subplot_title_fontsize)
      279    ax0m.text(0.95, 0.008, 'Nb of pix: {}'.format(nb_obs_trop_m),
      280                verticalalignment='bottom', horizontalalignment='right',
      281                transform=ax0m.transAxes,
      282                color='red', fontsize=5)
      283    ax0m.set_extent([lon_min-1, lon_max+1,
      284                     lat_min-1, lat_max+1], crs=projection)
      285    ax0m.coastlines(color='black',    resolution='50m',  linewidth=0.3)
      286    ax0m.add_feature(cfeature.BORDERS.with_scale('50m'), linewidth=0.2)
      287    gl_ax0m = ax0m.gridlines(
      288                              crs=projection,
      289                              draw_labels=True,
      290                              linewidth=0.4, color='gray', alpha=0.5, linestyle='--')
      291    # Make label disappear top, bottom, left, right
      292    gl_ax0m.xlabels_bottom = False
      293    gl_ax0m.ylabels_right  = False
      294    # Change format type
      295    gl_ax0m.xformatter = LONGITUDE_FORMATTER
      296    gl_ax0m.yformatter = LATITUDE_FORMATTER
      297    # Monitor font options
      298    gl_ax0m.xlabel_style = {'size': 5, 'color': 'gray'}
      299    gl_ax0m.ylabel_style = {'size': 5, 'color': 'gray'}
      300
      301    ## color-in the figure
      302    # Display each satellite pixel as individual polygones
      303    for i, toplot in enumerate(list_vrts_m):
      304        collectionm = PolyCollection(toplot, array=list_z_m[i], cmap=cmap)
      305        ax0m.add_collection(collectionm)
      306    # Manage colorbar
      307    divider0m = make_axes_locatable(ax0m)
      308    cax0m = divider0m.append_axes("right", size="3%", pad=colorbar_pad, axes_class=maxes.Axes)
      309    cb_ax0m = figm.colorbar(collectionm,
      310                             cax=cax0m,
      311                             ax=ax0m,
      312                             orientation='vertical',
      313                             extend='both',)
      314    cb_ax0m.set_label(     "mixing ratio (ppb)",
      315                             fontsize=colorbar_title_fontsize)
      316    cb_ax0m.ax.tick_params(direction='in',
      317                             labelsize=colorbar_tick_fontsize)
      318    cb_ax0m.formatter.set_powerlimits((0, 0))
      319    cb_ax0m.ax.yaxis.set_offset_position('left')
      320    cb_ax0m.update_ticks()
      321
      322    plt.savefig(savefig_name_monthly, bbox_inches ="tight", orientation ='landscape', dpi=240)
      323    print('File created : ', savefig_name_monthly)
      324    print('#-----------------------------------------------------------------------')
      325
      
    2. put TROPOMI data into a CIF monitor file:

      Show/Hide Code

        1###################################################################
        2# This code is provided by: I. Pison and R. Plauchu
        3# contact: isabelle.pison@lsce.ipsl.fr, robin.plauchu@lsce.ipsl.fr
        4# date: June 2022
        5# project: H2020 VERIFY, ANR ARGONAUT, CNES Tosca ARGOS
        6# relation to CIF: uses plugins, not embedded in.
        7# It is an example, to be used as a basis
        8# for developing other work.
        9# No guarantee as to its compatibility or results
       10# is attached since it is not part of the CIF
       11# Please don't forget acknowledgements if you use it.
       12####################################################################
       13
       14from netCDF4 import Dataset
       15from pycif.utils.datastores import dump
       16from pycif.utils.netcdf import readnc
       17from pycif.utils.datastores import empty
       18from pycif.utils.datastores.dump import dump_datastore
       19import calendar
       20import os
       21import pandas as pd
       22import datetime
       23import numpy as np
       24import xarray as xr
       25import netCDF4 as nc
       26
       27outdir = '/home/chimereicos/VERIFYdataWP4/ObsTROPOMI'
       28# Time period
       29year   = 2018
       30yearn = 2021
       31month  = 7
       32monthn = 2
       33day    = 1
       34dayn   = 28
       35
       36##### raw data to read
       37dir_sat = '/home/chimereges2/ipison/testdownload/'
       38spec = 'CH4'
       39# Satellite data
       40mission_name       = 'S5P'
       41processing_stream  = '*'
       42product_identifier = 'L2__{}____'.format(spec)
       43sat_data_id = mission_name+'_'+processing_stream+'_'+product_identifier
       44file_list_name    = 'filelist_{}_'.format(spec)
       45fixed_nlevels = 12 
       46# selec time slot when the satellite goes over the targeted area
       47hh_min = '14'
       48hh_max = '10'
       49
       50####### domain for filtering spatially
       51domain = 'EUROCOM'
       52print('---------------------------------------------------------------')
       53print("Start read COORD")
       54dirdom = '/home/users/afortems/PYTHON/PYVAR_CHIMERE_svn_Europe_ok/reposcode/reposcode/project_example/domains/'
       55coord_file        = '{}/HCOORD/COORD_{}'.format(dirdom,domain)
       56coord_corner_file = '{}/HCOORD//COORDcorner_{}'.format(dirdom,domain)
       57# Domain size
       58listdom = '{}/domainlist.nml'.format(dirdom)
       59f = open(listdom, 'r')
       60availabledom = f.readlines()
       61for line in availabledom:
       62    tmp = line.split() #str.split(line)
       63    if tmp[0] == domain:
       64        nlon_target = int(tmp[1])
       65        nlat_target = int(tmp[2])
       66#print('nlon, nlat :',nlon_target,nlat_target)
       67# Reading coord files
       68f1 = open(coord_corner_file, 'r')
       69coord1 = f1.readlines()
       70f1.close()
       71# Putting coords into matrix of coord corners
       72zlonc_target = np.zeros((nlon_target + 1, nlat_target + 1), float)
       73zlatc_target = np.zeros((nlon_target + 1, nlat_target + 1), float)
       74k = 0
       75for i in range(nlat_target + 1):
       76    for j in range(nlon_target + 1):
       77        tmp = coord1[k].split() #str.split(coord1[k])
       78        k += 1
       79        zlonc_target[j, i] = float(tmp[0])
       80        zlatc_target[j, i] = float(tmp[1])
       81lon_min, lon_max = zlonc_target.min(), zlonc_target.max()
       82lat_min, lat_max = zlatc_target.min(), zlatc_target.max()
       83print("End read COORD")
       84
       85##### datastore to dump in netcdf file for use by the CIF
       86# list of infos required for any type of obs
       87list_basic_cols = [ 'date', 'duration', 'station' , 'network', 'parameter', 'lon', 'lat', 'obs', 'obserror' ]
       88# fixed information for station, network and parameter + duration columns
       89stationID = 'TROPOMI'
       90networkID = 'KNMISRON'
       91spec = 'CH4'
       92obsduration = 1./60. 
       93
       94# loop over the available files
       95for yy in range(year, yearn+1):
       96  print(yy)
       97  endm = 12
       98  begm = 1
       99  if (yy == yearn): endm = monthn
      100  if (yy == year) : begm = month
      101  for mm in range(begm, endm+1):
      102    print(mm)
      103    endday = ndays = calendar.mdays[mm] + ((mm == 2) and calendar.isleap(yy))
      104    begday = 1
      105    if (mm == monthn): endday = dayn
      106    if( mm == month ) : begday = day
      107    for dd in range(begday,endday+1):
      108      print( mm, dd)
      109   
      110      # List of netcdf files for one day included in [hh_min; hh_max] as
      111      # an aera may be covered on several scanlines while the satellite is orbiting
      112      datetime0 = datetime.datetime(year=yy, month=mm, day=dd, hour=int(hh_min), minute=0)
      113      os.system('ls ' + dir_sat + sat_data_id + datetime0.strftime("%Y%m%dT") + '{' + hh_min + '..' + hh_max + '}*' +\
      114         datetime0.strftime("%Y%m%dT") + '*.nc > ' + dir_sat + file_list_name + str(yy)+"%2.2i"%mm+"%2.2i"%dd+'.txt')
      115      ff = open(dir_sat + file_list_name + str(yy)+"%2.2i"%mm+"%2.2i"%dd+'.txt','r')
      116      sat_files = ff.readlines()
      117      ff.close()
      118      print('#-----------------------------------------------------------------------')
      119      for theline in sat_files:
      120        print(theline)
      121      print('#-----------------------------------------------------------------------')
      122      if len(sat_files) > 1:
      123        print('More than 1 orbit is selected')
      124      if len(sat_files) == 0: continue
      125      # Read each netcdf file and select variables of interest
      126      count = 0
      127      nb_obs_trop = 0
      128      for theline in sat_files:
      129        count += 1
      130        print('#-----------------------------------------------------------------------')
      131        print('File n°:',count,'/',len(sat_files))
      132        #print(str.split(theline,'/')[-1])
      133
      134        s = str.split(theline)
      135
      136        filein = nc.Dataset(str.strip(theline), 'r')
      137        # read groups and such in netcdf4 file
      138        product_grp = filein.groups['PRODUCT'] 
      139        support_data_grp = product_grp.groups['SUPPORT_DATA']
      140        detailed_results_grp = support_data_grp.groups['DETAILED_RESULTS']
      141        geoloc_grp  = support_data_grp.groups['GEOLOCATIONS']
      142        input_data_grp       = support_data_grp.groups['INPUT_DATA']
      143        # lon-lat to create spatial mask for the domain
      144        # avoid too large tables in the following + shape = list of valid data + levels if required
      145        longitude = product_grp.variables['longitude'][:,:,:] 
      146        latitude = product_grp.variables['latitude'][:,:,:] 
      147        mask = (   (longitude < lon_max) & \
      148                   (longitude > lon_min) & \
      149                   ( latitude < lat_max) & \
      150                   ( latitude > lat_min) )
      151        
      152        ### fill in sub-datastore for the current file
      153        subdata = pd.DataFrame( columns = list_basic_cols )
      154    
      155        ## check number of levels 
      156        ## using: averaging kernels
      157        ak = detailed_results_grp.variables['column_averaging_kernel'][:,:,:,:] # time, scanline, ground_pixel, layer
      158        nlevread = ak.shape[3]
      159        if(nlevread != fixed_nlevels) : 
      160            print('ERROR NUMBER OF LEVELS')
      161            break
      162        
      163        ##date: starting date
      164        # must be a datetime object
      165        time_utc     = product_grp.variables['time_utc'][:,:] # time, scanline with time =1 always for S5P
      166        time_utc[time_utc=='']= '1870-01-01T01:01:00.00Z' #case of files with void dates!!
      167        # Trick to expand time_utc to time, scanline, ground_pixel shape
      168        time_utc_exp = np.repeat(time_utc[:, :, np.newaxis], ak.shape[2], axis=2)
      169        time_utc_exp = time_utc_exp[:,:,:][mask]
      170        dates = [datetime.datetime.strptime(time_utc_exp[i], \
      171                    '%Y-%m-%dT%H:%M:%S.%fZ').replace(microsecond=0) \
      172                     for i in range(len(time_utc_exp))]
      173        subdata['date'] = dates
      174        ## duration    duration in hours
      175        subdata = subdata.assign(duration = obsduration) 
      176        ## lon    longitude of the measurement
      177        ## lat    latitude of the measurement
      178        subdata['lon'] = longitude[mask]
      179        subdata['lat'] = latitude[mask]
      180        ## obs    observed value
      181        obs = product_grp.variables['methane_mixing_ratio'][:,:,:]
      182        subdata['obs'] = np.float64(obs[mask] )
      183        ## obserror  error on the observation
      184        obs_err = np.float64(product_grp.variables['methane_mixing_ratio_precision'][:,:,:]) # 1-sigma
      185        subdata['obserror'] = obs_err[mask] 
      186        subdata = subdata.assign(station = stationID)
      187        subdata = subdata.assign(network = networkID)
      188        subdata = subdata.assign(parameter = spec)
      189
      190        ## pressure levels: 
      191        # particular case of TROPOM CH4: create pressure grid from psurf and delta pressure
      192        psurf = input_data_grp.variables['surface_pressure'][:,:,:] # time, scanline, ground_pixel
      193        deltap = input_data_grp.variables['pressure_interval'][:,:,:] # time, scanline, ground_pixel
      194        pavg0 = np.array([psurf - deltap * lev for lev in range(fixed_nlevels + 1)])
      195        pavg = np.transpose(pavg0, (1,2,3,0))
      196        pavg = pavg[mask]
      197        pavg = np.float64(pavg)
      198        ## prior profile
      199        qa0 = input_data_grp.variables['methane_profile_apriori'][:,:,:,:] # time, scanline, ground_pixel, layer
      200        conv = input_data_grp.variables['methane_profile_apriori'].multiplication_factor_to_convert_to_molecules_percm2
      201        qa0 = qa0[mask] * conv
      202        qa0 = np.float64(qa0)
      203        ## averaging kernels
      204        ak = ak[mask]
      205        ak = np.float64(ak)
      206
      207        ## information particular to TROPOMI
      208        dvair = input_data_grp.variables['dry_air_subcolumns'][:,:,:,:]
      209        conv = input_data_grp.variables['dry_air_subcolumns'].multiplication_factor_to_convert_to_molecules_percm2 
      210        dvair = dvair[mask] * conv
      211        dvair = np.float64(dvair)
      212
      213        ## put everything in a xarray
      214        # satellite specific info
      215        ds = xr.Dataset({'qa0': (['index', 'level'], qa0),
      216                     'ak': (['index', 'level'], ak),
      217                     'pavg0': (['index', 'level_pressure'], pavg),
      218                     'dryair': (['index', 'level'], dvair)},
      219                     coords={'index': subdata.index,
      220                             'level': range(fixed_nlevels),
      221                              'level_pressure': range(fixed_nlevels + 1)})
      222
      223        # put together basic info and satellite info
      224        subdata = subdata.to_xarray()
      225        subdata = xr.merge([ds, subdata])
      226        ## filtering according to information provided with data
      227        Qa = 0.5 # strictly superior to 0.5
      228        Qa_read = product_grp.variables['qa_value'][:]
      229        # filter on Qa, already filtered over land and for zenithal angles < 60 deg
      230        subfilter = pd.DataFrame(columns = ['Qa'])
      231        subfilter['Qa'] = Qa_read[mask]
      232        mask_filter = np.where((subfilter['Qa'] > Qa))[0]
      233        subdata = subdata.isel(index=mask_filter)
      234        
      235        ## index = number of the observation among the remaining ones 
      236        subdata["index"] = np.arange(len(subdata["index"]))
      237        if len(subdata["index"] > 0 :
      238          ## dump to netCDF - here, one per day
      239          subdata.to_netcdf('{}/monitor_{}{}_{}_{}{}{}_orbit{}.nc'.format(outdir,stationID, networkID, spec, yy, "%2.2i"%mm, "%2.2i"%dd, count))
      240
      241
      242
      
    3. filter monitor.nc after a first forward simulation:

      Show/Hide Code

       1#########################################################
       2# This code is provided by: I. Pison
       3# contact: isabelle.pison@lsce.ipsl.fr
       4# date: June 2022
       5# project: H2020 VERIFY, ANR ARGONAUT, CNES Tosca ARGOS
       6# relation to CIF: uses plugins, not embedded in.
       7# It is an example, to be used as a basis
       8# for developing other work.
       9# No guarantee as to its compatibility or results
      10# is attached since it is not part of the CIF
      11# Please don't forget acknowledgements if you use it.
      12#########################################################
      13
      14from pycif.utils.datastores import dump
      15import copy
      16import numpy as np
      17import pandas as pd
      18from pycif.utils.netcdf import readnc
      19import xarray as xr
      20from itertools import zip_longest
      21
      22def retrieve_aks(fin,index,list_var): 
      23
      24   # fin is the original input monitor 
      25   # containing satellite-specific information: ak, qa0, etc
      26   # fin is read to retrieve this info for data number index
      27   index = int(index)
      28   qa0, ak, pavg0, dryair = readnc(fin, list_var)
      29
      30   return qa0[index].tolist(), ak[index].tolist(), pavg0[index].tolist(), dryair[index].tolist()
      31
      32# WORKDIR of the forward simulation with all "raw" data
      33refdir = '/home/chimereicos/fwd_all_data/obsoperator/fwd_0000/'
      34# directory for the info file
      35dirinfo = 'chain/satellites/default_00001'
      36# directory for the output monitor file
      37dirmonit = 'obsvect/satellites/CH4/'
      38# files to use
      39monitor_in = '{}{}/monitor.nc'.format(refdir,dirmonit)
      40infos = '{}{}/infos_201901010000.nc'.format(refdir,dirinfo)
      41
      42# basic data to provide in an input monitor
      43list_basic_cols = [ 'date', 'duration', 'station', 'network', 'parameter', 'lon', 'lat', 'obs', 'obserror' ]
      44
      45# Monitor and info data to use
      46ds = dump.read_datastore(monitor_in)
      47dsinf = dump.read_datastore(infos,col2dump= ['orig_index', 'file_aks'])
      48
      49# paste ds and dsinf together to get full information for each obs
      50ds3 = pd.concat( [ds,dsinf],axis = 1)
      51
      52# example filter: filter out obs outside the domain
      53ds3 = ds3.loc[ ~np.isnan(ds3['orig_index']) ] 
      54
      55# example change: set obserror at a given value
      56ds3 = ds3.assign(obserror = 20)
      57
      58# generate a satellite input monitor with these new data:
      59# for each line in ds2, retrieve qa0,ak,pavg0,dryair from the right file
      60list_var = ['qa0', 'ak', 'pavg0', 'dryair']
      61ds4 = ds3.apply(lambda x: retrieve_aks(x.file_aks,x.orig_index,list_var), axis = 1)
      62# reformat ds4 to put satellite specific info in a xarray
      63# keep index to keep track of selection - no impact on the CIF
      64ds5 = pd.DataFrame((item for item in ds4), columns = list_var, index=ds3.index)
      65fixed_nlevels = len(ds5['ak'].iloc[0])
      66list_tab = []
      67for k,var in enumerate(list_var):
      68  list_tab.append(pd.DataFrame((item for item in ds5[var])).values)
      69# satellite specific info
      70ds_sat = xr.Dataset({'qa0': (['index', 'level'], list_tab[0]),
      71                     'ak': (['index', 'level'], list_tab[1]),
      72                     'pavg0': (['index', 'level_pressure'], list_tab[2]),
      73                     'dryair': (['index', 'level'], list_tab[3])},
      74                     coords={'index': ds5.index,
      75                             'level': range(fixed_nlevels),
      76                              'level_pressure': range(fixed_nlevels + 1)})
      77# basic data 
      78basic_data = ds3[list_basic_cols].to_xarray()
      79
      80# merge basic and satellite-specific data
      81alldata = xr.merge([ds_sat, basic_data])
      82
      83# create new clean input monitor
      84alldata.to_netcdf('new_monitor.nc')
      85
      86
      
    4. run the job for year 2019:

    Show/Hide Code

      1#####################
      2# PYCIF config file #
      3#####################
      4# PYCIF parameters
      5
      6# Verbose level
      7verbose : 1
      8
      9# Log file (to be saved in $wordkir)
     10logfile: first-CH4-simu
     11
     12# Execution directory
     13workdir: /home/chimereicos/VERIFYoutputs/inv-CH4-withstrato-m1qn3-still-better-correrr/
     14
     15# Dates
     16datei : 2019-01-01 00:00:00
     17datef : 2020-01-01 00:00:00
     18
     19#####################################################################
     20# Running mode for PYCIF
     21mode:
     22  plugin:
     23    name: 4dvar
     24    version: std
     25    type: mode
     26  save_out_netcdf: True
     27  minimizer:
     28    plugin:
     29      name: M1QN3
     30      version: std
     31      type: minimizer
     32    maxiter: 10
     33    epsg: 0.1
     34    df1: 0.01
     35    simulator:
     36      plugin:
     37        name: gausscost
     38        version: std
     39      reload_from_previous: True
     40
     41#####################################################################
     42# Observation operator
     43obsoperator:
     44  plugin:
     45    name: standard
     46    version: std
     47  autorestart: True
     48#####################################################################
     49controlvect:
     50  plugin:
     51    name: standard
     52    version: std
     53  save_out_netcdf: True
     54
     55#####################################################################
     56# Transport model
     57model :
     58  plugin:
     59    name    : CHIMERE
     60    version : std
     61
     62  # Executable
     63  direxec: /home/users/ipison/cif/model_sources/chimereGES/
     64  #  Number of physical steps per hour
     65  nphour_ref : 6
     66  # Number of chemical refined iterations
     67  ichemstep : 2
     68  # Deep convection
     69  ideepconv: 0
     70  #  Number of vertical layers in output files
     71  nivout: 17
     72  # Number of levels for emissions
     73  nlevemis: 1
     74
     75  # Number of MPI subdomains in zonal and meridian directions for // use
     76  nzdoms : 5
     77  nmdoms : 2
     78
     79  ignore_input_dates: True
     80  force_clean_run: True
     81  useRAMonly: True
     82
     83####################################################################
     84obsvect:
     85  plugin:
     86    name: standard
     87    version: std 
     88
     89#####################################################################
     90platform:
     91  plugin:
     92    name: LSCE
     93    version: obelix
     94
     95#####################################################################
     96domain :
     97  plugin:
     98    name    : CHIMERE
     99    version : std
    100  repgrid: /home/satellites14/mkouyate/CHIMERE/domains/domains/
    101  domid : EUROCOM
    102  nlev: 17
    103  p1: 997
    104  pmax: 200
    105
    106#####################################################################
    107chemistry :
    108  plugin:
    109    name: CHIMERE
    110    version: gasJtab
    111  schemeid: ges.espigrad
    112  dir_precomp :  /home/satellites14/mkouyate/CHIMERE/CIF_chimere_fwd_test2/
    113
    114
    115#####################################################################
    116datavect:
    117  plugin:
    118    name: standard
    119    version: std
    120
    121  components:
    122    meteo:
    123      plugin:
    124        name: CHIMERE
    125        version: std
    126        type: meteo
    127      dir:  /home/satellites14/mkouyate/CHIMERE/PYCIF_DATA/%Y/CHIMERE/EUROCOM/CH4/
    128      file: METEO.%Y%m%d%H.24.nc
    129      file_freq: 1D
    130    flux:   
    131      plugin:
    132        name: CHIMERE
    133        version: AEMISSIONS
    134        type: flux
    135      dir: /home/satellites14/mkouyate/CHIMERE/PYCIF_DATA/%Y/CHIMERE/EUROCOM/CH4/
    136      file: AEMISSIONS.%Y%m%d%H.24.nc
    137      file_freq: 1D
    138      parameters:
    139        CH4:
    140          hresol: hpixels
    141          vresol: vpixels
    142          tresol: 7D
    143          nlev: 1
    144          err: 1.
    145          hcorrelations: 
    146             sigma_land: 50
    147             sigma_sea: 50
    148             landsea: True
    149             filelsm: /home/satellites14/mkouyate/CHIMERE/PYCIF_DATA/eurocom_lndseamask.nc
    150             dump_hcorr: True
    151          tcorrelations:
    152             sigma_t: 7D
    153             dump_tcorr: True
    154    inicond:
    155      plugin:
    156        name: CHIMERE
    157        version: icbc
    158        type: field
    159      dir: /home/satellites14/mkouyate/CHIMERE/PYCIF_DATA/2017/CHIMERE/EUROCOM/CH4/ 
    160      file: INI_CONCS.20170101_20170101.00_EUROCOM.nc 
    161      comp_type: inicond    
    162      parameters:
    163        CH4:
    164          hresol: hpixels
    165          vresol: vpixels
    166          nlev: 17
    167          err: 0.02
    168    latcond:
    169      plugin:
    170        name: CHIMERE
    171        version: icbc
    172        type: field
    173      comp_type: latcond
    174      dir: /home/satellites14/mkouyate/CHIMERE/PYCIF_DATA/2017/CHIMERE/EUROCOM/CH4/
    175      file: BOUN_CONCS.2017%m%d%H.24.nc 
    176      file_freq: 1D
    177      parameters:
    178        CH4:
    179          hresol : bands 
    180          bands_lat : [ 30, 53, 75 ]
    181          bands_lon : [-18, 10, 38 ]
    182          is_lbc: True
    183          vresol: vpixels
    184          nlev: 17 
    185          tresol: 2D
    186          err: 0.02
    187          hcorrelations: 
    188             sigma: 100 
    189          tcorrelations: 
    190             sigma_t: 120H
    191    topcond:
    192      plugin:
    193        name: CHIMERE
    194        version: icbc
    195        type: field
    196      comp_type: topcond
    197      dir: /home/satellites14/mkouyate/CHIMERE/PYCIF_DATA/2017/CHIMERE/EUROCOM/CH4/ 
    198      file: BOUN_CONCS.2017%m%d%H.24.nc
    199      file_freq: 1D
    200      parameters:
    201        CH4:
    202          hresol: global
    203          tresol: 2D
    204          err: 0.02
    205          tcorrelations:
    206            sigma_t: 120H
    207   
    208    satellites:
    209      parameters:
    210        CH4:
    211          dir: /home/chimereicos/VERIFYoutputs/
    212          file: secondfilter_monitor_TROPOMI.nc
    213          formula: 5
    214          pressure: Pa
    215          product: level
    216          cropstrato: True
    217          correct_pthick: False 
    218          fill_strato: True
    219          molmass: 16.
    220          nchunks: 100
    221
    222    stratosphere:
    223      parameters:
    224        CH4:
    225          dir: /home/chimereges/ecmwf/camsCH4v17r1/
    226          file: z_cams_l_tno-sron_2017%m_v17r1_ra_ml_6h_ch4.nc
    227          varname: CH4
    228          plugin:
    229            name: CAMS
    230            type: field
    231            version: netcdf          
    232          regrid:
    233            method: bilinear
    234          unit_conversion: # CAMS = for MMR*1e-9 to ppb for CHIMERE
    235            scale: 1.
    236          file_freq: 1MS
    237          vresol: column
    238          tresol: 2D
    239          tcorrelations:
    240            sigma_t: 120H
    241          hresol: hpixels
    242          hcorrelations: 
    243             sigma: 100
    244          err: 0.02
    
    1. post-processing includes:

    1. maps of fluxes and increments: see example for surface inversions

    2. comparison of columns (obs to prior and obs to post) as scatter plots:

    Show/Hide Code

     1#####################################################
     2# This code is provided by: I. Pison
     3# contact: isabelle.pison@lsce.ipsl.fr
     4# date: June 2022
     5# project: H2020 VERIFY, ANR ARGONAUT, CNES ARGOS
     6# relation to CIF: uses plugins, not embedded in.
     7# It is an example, to be used as a basis
     8# for developing other work.
     9# No guarantee as to its compatibility or results
    10# is attached since it is not part of the CIF
    11# Please don't forget acknowledgements if you use it.
    12#####################################################
    13
    14import numpy as np
    15import pandas as pd
    16from pycif.utils.netcdf import readnc
    17from pycif.utils.datastores import dump
    18import matplotlib.pyplot as plt
    19import datetime
    20
    21dir_results = '/home/chimereicos/VERIFYoutputs/inv-CH4-withstrato-m1qn3-still-better-correrr_2/'
    22dir_out = '/obsoperator/fwd_'
    23dir_monit = '/obsvect/satellites/CH4/'
    24
    25year_beg = 2019
    26year_end = 2019
    27
    28# read monitors of all years
    29ds_all_years = pd.DataFrame()
    30
    31for yy in range(year_beg, year_end + 1, 1):
    32  print(yy)
    33  fmonitpr = '{}{}-001{}/monitor.nc'.format(dir_results, dir_out, dir_monit)
    34  fmonitpo = '{}{}-002{}/monitor.nc'.format(dir_results, dir_out, dir_monit)
    35  dspr = dump.read_datastore(fmonitpr)
    36  dspo = dump.read_datastore(fmonitpo)
    37  dsprpo = pd.concat( [dspr, dspo['maindata']['sim']], axis = 1)
    38  ds_all_years = pd.concat( [dsprpo, ds_all_years], axis = 0)
    39
    40# scatter plot sim vs obs
    41obs = ds_all_years[('maindata','obs')]
    42simpr = ds_all_years[('maindata', 'sim')]
    43simpo = ds_all_years['sim']
    44
    45coef_pr = np.polyfit(obs, simpr, 1)
    46poly1d_fn_pr = np.poly1d(coef_pr) 
    47corr_pr = np.corrcoef(obs, simpr)[0, 1]
    48coef_po = np.polyfit(obs, simpo, 1)
    49poly1d_fn_po = np.poly1d(coef_po) 
    50corr_po = np.corrcoef(obs, simpo)[0, 1]
    51
    52xmin, xmax, ymin, ymax = 1000, 2500, 1000, 2500
    53
    54fig, axs = plt.subplots(ncols=2, sharey=True, figsize=(21, 11))
    55fig.subplots_adjust(hspace=0.5, left=0.07, right=0.93)
    56fig.suptitle('Methane columns: TROPOMI data and CHIMERE simulated equivalents')
    57ax = axs[0]
    58hb = ax.hexbin(obs, simpr, gridsize=100, bins="log", cmap='cividis')
    59ax.plot([xmin, xmax], poly1d_fn_pr([xmin, xmax]), '--k', linewidth=5,
    60        label="y={:.2f}+{:.2f}x; r={:.2f}".format(coef_pr[1], coef_pr[0], corr_pr) )
    61ax.legend()
    62ax.axis([xmin, xmax, ymin, ymax])
    63ax.set_title("Prior vs observations: density of points")
    64ax.set_ylabel('Prior (ppb)')
    65ax.set_xlabel('TROPOMI data (ppb)')
    66cb = fig.colorbar(hb, ax=ax)
    67cb.set_label('log10(N)')
    68
    69ax = axs[1]
    70hb = ax.hexbin(obs, simpo, gridsize=100, bins='log', cmap='cividis')
    71ax.plot([xmin, xmax], poly1d_fn_po([xmin, xmax]), '--k', linewidth=5,
    72        label="y={:.2f}+{:.2f}x; r={:.2f}".format(coef_po[1], coef_po[0], corr_po))
    73ax.legend()
    74ax.axis([xmin, xmax, ymin, ymax])
    75ax.set_title("Posterior vs observations: density of points")
    76ax.set_ylabel('Analysis (ppb)')
    77ax.set_xlabel('TROPOMI data (ppb)')
    78cb = fig.colorbar(hb, ax=ax)
    79cb.set_label('log10(N)')
    80
    81plt.savefig('hexbin_TROPOMI_2019_OK.png')
    82plt.close()
    83
    
    1. maps of columns data/prior/posterior and differences:

    Show/Hide Code

      1#####################################################
      2# This code is provided by: I. Pison
      3# contact: isabelle.pison@lsce.ipsl.fr
      4# date: June 2022
      5# project: H2020 VERIFY, ANR ARGONAUT, CNES ARGOS
      6# relation to CIF: uses plugins, not embedded in.
      7# It is an example, to be used as a basis
      8# for developing other work.
      9# No guarantee as to its compatibility or results
     10# is attached since it is not part of the CIF
     11# Please don't forget acknowledgements if you use it.
     12#####################################################
     13
     14import pandas as pd
     15import matplotlib.pyplot as plt
     16import cartopy.crs as crs
     17import numpy as np
     18from pycif.utils.datastores import dump
     19from mpl_toolkits.axes_grid1.inset_locator import InsetPosition
     20from cartopy.feature import ShapelyFeature
     21from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
     22
     23dir_results = '/home/chimereicos/VERIFYoutputs/inv-CH4-withstrato-m1qn3-still-better-correrr/'
     24dir_out = '/obsoperator/fwd_'
     25dir_monit = '/obsvect/satellites/CH4/'
     26
     27year_beg = 2019
     28year_end = 2019
     29
     30# read monitors of all years
     31ds_all_years = pd.DataFrame()
     32
     33for yy in range(year_beg, year_end + 1, 1):
     34  print(yy)
     35  fmonitpr = '{}{}-001{}/monitor.nc'.format(dir_results, dir_out, dir_monit)
     36  fmonitpo = '{}{}-002{}/monitor.nc'.format(dir_results, dir_out, dir_monit)
     37  print(fmonitpo)
     38  dspr = dump.read_datastore(fmonitpr)
     39  dspo = dump.read_datastore(fmonitpo)
     40  dsprpo = pd.concat( [dspr, dspo['maindata']['sim']], axis = 1)
     41  ds_all_years = pd.concat( [dsprpo, ds_all_years], axis = 0)
     42
     43
     44lon = ds_all_years[('metadata', 'lon')]
     45lat = ds_all_years[('metadata', 'lat')]
     46conc_obs = ds_all_years[('maindata', 'obs')]
     47conc_prior = ds_all_years[('maindata', 'sim')]
     48conc_post = ds_all_years['sim']
     49
     50#############
     51# plot maps #
     52#############
     53# color scale from existing colormaps
     54Ncolor = 15
     55cmap = plt.cm.jet
     56cmaplist = cmap(np.linspace(0, 1, Ncolor))
     57cmap_plot = cmap.from_list('custom', cmaplist, Ncolor)
     58
     59Ncolor2 = 11
     60cmap2 = plt.cm.bwr
     61cmaplist2 = cmap2(np.linspace(0, 1, Ncolor2))
     62cmap2 = cmap2.from_list('custom', cmaplist2, Ncolor2)
     63
     64cmin, cmax = conc_obs.min(), conc_obs.max()
     65
     66fig = plt.figure(figsize=(21, 10))
     67gs = fig.add_gridspec(nrows = 2, ncols = 4, width_ratios = [2,2,2,0.1], 
     68                 left=0., right=0.97, bottom=0, top=1,
     69                  height_ratios = [1, 1 ] ,hspace = 0.01, wspace = 0.25)
     70
     71subobs = fig.add_subplot(gs[0,0], projection=crs.PlateCarree())
     72subobs.scatter(lon, lat, c=conc_obs, s=5,
     73                        vmin=cmin, vmax=cmax, cmap = cmap_plot,
     74                        transform=crs.PlateCarree())
     75subobs.coastlines()
     76subobs.set_title("TROPOMI data (ppb)")
     77
     78subsim = fig.add_subplot(gs[0,1], projection=crs.PlateCarree())
     79im = subsim.scatter(lon, lat, c=conc_prior, s=5,
     80                        vmin=cmin, vmax=cmax, cmap = cmap_plot,
     81                        transform=crs.PlateCarree())
     82subsim.coastlines()
     83subsim.set_title("Prior simulations (ppb)")
     84
     85subpost = fig.add_subplot(gs[0,2], projection=crs.PlateCarree())
     86subpost.scatter(lon, lat, c=conc_post, s=5,
     87                        vmin=cmin, vmax=cmax, cmap = cmap_plot,
     88                        transform=crs.PlateCarree())
     89subpost.coastlines()
     90subpost.set_title("Posterior simulations (ppb)")
     91
     92echelle = fig.add_subplot(gs[0,3])
     93ip = InsetPosition(echelle, [0,0,1,1]) 
     94echelle.set_axes_locator(ip)
     95p00 = echelle.get_position()
     96p01 = echelle.get_position()
     97p02 = subpost.get_position()
     98p00_new = [p00.x0, p02.y0, p00.width, p02.height]
     99echelle.set_position(p00_new)
    100cbar = plt.colorbar(im, cax=echelle, ax = [subobs , subpost] )
    101cbar.ax.set_title('ppb') 
    102
    103dcmin, dcmax = -20, 20
    104dcmin,dcmax = (conc_post - conc_obs).min(), (conc_post - conc_obs).max()
    105dcmin, dcmax = -100, 100
    106dsubsim = fig.add_subplot(gs[1,1], projection=crs.PlateCarree())
    107im2 = dsubsim.scatter(lon, lat, c=conc_prior - conc_obs, s=5,
    108                        vmin=dcmin, vmax=dcmax,
    109                        cmap = cmap2,
    110                        transform=crs.PlateCarree())
    111dsubsim.coastlines()
    112dsubsim.set_title("Prior - TROPOMI data")
    113
    114dsubpost = fig.add_subplot(gs[1,2], projection=crs.PlateCarree())
    115dsubpost.scatter(lon, lat, c=conc_post - conc_obs, s=5,
    116                        vmin=dcmin, vmax=dcmax,
    117                        cmap = cmap2,
    118                        transform=crs.PlateCarree())
    119dsubpost.coastlines()
    120dsubpost.set_title("Post - TROPOMI data ")
    121echelle2 = fig.add_subplot(gs[1,3])
    122ip = InsetPosition(echelle2, [0,0,1,1]) 
    123echelle2.set_axes_locator(ip)
    124p00 = echelle2.get_position()
    125p01 = echelle2.get_position()
    126p02 = dsubpost.get_position()
    127p00_new = [p00.x0, p02.y0, p00.width, p02.height]
    128echelle2.set_position(p00_new)
    129cbar2 = plt.colorbar(im2, cax=echelle2, ax = dsubpost)
    130cbar2.ax.set_title('%') 
    131
    132
    133plt.savefig('obsandanalysis.png')
    134
    
    1. time series of emissions, budgets per country: see example for surface inversions.

Nitrous oxide net total emissions in Europe (WP4)