Inversions with TROPOMI data

1. Plot raw TROPOMI data

Show/Hide Code

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

2. Put TROPOMI data into a CIF monitor file

Show/Hide Code

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

3. Filter monitor.nc after a first forward simulation

Show/Hide Code

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

4. Run the job for year 2019

Show/Hide Code

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

5. Post-processing includes:

a. Maps of fluxes and increments

See example for surface inversions.

b. 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

c. 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

d. Time series of emissions, budgets per country

See example for surface inversions.