Inversions with surface stations¶
1. Prepare the data to assimilate¶
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
2. Automatically make one job per year¶
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
To fill-in a .sed file, you may use cat to replace the variables, as illustrated in the script for a job in 2015:
Show/Hide Code
1#!/bin/sh -login 2 3echo "Machine on which the job is running:" 4hostname 5 6# User' choices: 7# Give year to run 8year=2015 9# Give directory for the .yml.sed file 10dirsed=$HOME/VERIFY/ 11 12# next year for end of simulation 13ynext=$[ year + 1 ] 14 15# Path of the yml.sed file to use 16fyml=${dirsed}/config_chimere_inv_1Y 17 18# Fill-in the yaml.sed into a .yml 19cat << EOD > ${fyml}${year}.yml.sed-commands 20s,_YYYY_,${year}, 21s,_YYNN_,${ynext}, 22s,_INI_,INI_CONCS.${year}0101_${year}0101.00_EUROCOM.nc, 23EOD 24sed -f ${fyml}${year}.yml.sed-commands ${fyml}.yaml.sed > ${fyml}${year}.yml 25 26# Run the CIF 27time python -m pycif ${fyml}${year}.yml
3. Post-processing¶
Post-processing of the results include:
a. 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
b. Time series of mixing ratios and maps of statistical indicators at the stations¶
Time series of mixing ratios at the stations make it easy to compare the observation, the prior and the analysis (obtained with the posterior emissions).
The maps of statistical indicators at the stations show the ratios of the posterior to the prior values of the bias, RMS, correlation for example.
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
c. Yearly emission budgets for various regions¶
Examples here of regions are 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
d. Time series of emissions for various regions¶
Examples here of regions are 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
