Inversions in the project VERIFY¶
Methane net total emissions in Europe (WP4)¶
inversions with surface stations:
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
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
post-processing of the results include:
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
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
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
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() 77inversions with TROPOMI data:
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('#-----------------------------------------------------------------------') 325put 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 242filter 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 86run 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
post-processing includes:
maps of fluxes and increments: see example for surface inversions
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
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
time series of emissions, budgets per country: see example for surface inversions.
