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