Changes to MOPITT pairing

melodies_monet/driver.py

@@ -259,7 +259,7 @@ def open_sat_obs(self, time_interval=None, control_dict=None):
         None
         """
         from .util import time_interval_subset as tsub
-
+        from glob import glob
         try:
             if self.sat_type == 'omps_l3':
                 print('Reading OMPS L3')
@@ -287,6 +287,13 @@ def open_sat_obs(self, time_interval=None, control_dict=None):
                 else: flst = self.file
                 self.obj = mio.sat._mopitt_l3_mm.open_dataset(flst, ['column','pressure_surf','apriori_col',
                                                                           'apriori_surf','apriori_prof','ak_col'])
+
+                # Determine if monthly or daily product and set as attribute
+                if 'MOP03JM' in glob(self.file)[0]:
+                    self.obj.attrs['monthly'] = True
+                else:
+                    self.obj.attrs['monthly'] = False
+
             elif self.sat_type == 'modis_l2':
                 # from monetio import modis_l2
                 print('Reading MODIS L2')
@@ -440,6 +447,7 @@ def __init__(self):
         """Initialize a :class:`model` object."""
         self.model = None
         self.apply_ak = False
+        self.mod_to_overpass = False
         self.radius_of_influence = None
         self.mod_kwargs = {}
         self.file_str = None
@@ -936,6 +944,8 @@ def open_models(self, time_interval=None,load_files=True):
                 # set the model label in the dictionary and model class instance
                 if "apply_ak" in self.control_dict['model'][mod].keys():
                     m.apply_ak = self.control_dict['model'][mod]['apply_ak']
+                if "mod_to_overpass" in self.control_dict['model'][mod].keys():
+                    m.mod_to_overpass = self.control_dict['model'][mod]['mod_to_overpass']
                 if 'radius_of_influence' in self.control_dict['model'][mod].keys():
                     m.radius_of_influence = self.control_dict['model'][mod]['radius_of_influence']
                 else:
@@ -1395,11 +1405,20 @@ def pair_data(self, time_interval=None):
 
                     elif obs.sat_type == 'mopitt_l3':
                         from .util import satellite_utilities as sutil
+
                         if mod.apply_ak: 
                             model_obj = mod.obj[keys+['pres_pa_mid']]
                             # trim to only data within analysis window, as averaging kernels can't be applied outside it
                             obs_dat = obs.obj.sel(time=slice(self.start_time.date(),self.end_time.date()))#.copy()
                             model_obj = model_obj.sel(time=slice(self.start_time.date(),self.end_time.date()))#.copy()
+
+                            # Sample model to observation overpass time
+                            if mod.mod_to_overpass:
+                                print('sampling model to 10:30 local overpass time')
+                                overpass_datetime = pd.date_range(self.start_time.replace(hour=10,minute=30),
+                                                                  self.end_time.replace(hour=10,minute=30),freq='D')
+                                model_obj = sutil.mod_to_overpasstime(model_obj,overpass_datetime)
+
                             # interpolate model to observation, calculate column with averaging kernels applied
                             paired = sutil.mopitt_l3_pairing(model_obj,obs_dat,keys[0])
                             p = pair()

melodies_monet/util/sat_l2_swath_utility.py

@@ -143,7 +143,7 @@ def trp_interp_swatogrd_ak(obsobj, modobj):
     modlat = modobj.coords['latitude']
     modlon = modobj.coords['longitude']
 
-    tmpvalue = np.zeros([ny, nx], dtype = np.float)
+    tmpvalue = np.zeros([ny, nx], dtype = np.float64)
 
     time   = [ datetime.strptime(x,'%Y-%m-%d') for x in obsobj.keys()]
     ntime  = len(list(obsobj.keys()))
@@ -291,18 +291,16 @@ def cal_amf_wrfchem(scatw, wrfpreslayer, tpreslev, troppres, wrfno2layer_molec,
     vertical_pres = []
     vertical_scatw = []
     vertical_wrfp = []
-
+    
     for llb in range(len(lb[0])):
         yy = lb[0][llb]
         xx = lb[1][llb]
         vertical_pres = tpreslev[:,yy,xx] # mli, update to new dimension
         vertical_scatw = scatw[yy,xx,:]
         vertical_wrfp = wrfpreslayer[yy,xx,:]
-
         f = interpolate.interp1d(np.log10(vertical_pres[:]),vertical_scatw[:], fill_value="extrapolate")# relationship between pressure to avk
         wrfavk[yy,xx,:] = f(np.log10(vertical_wrfp[:])) #wrf-chem averaging kernel
 
-
     for l in range(nz-1):
         # check if it's within tropopause
         preminus[:,:]         = wrfpreslayer[:,:,l] - troppres[:,:]

melodies_monet/util/satellite_utilities.py

@@ -33,6 +33,56 @@ def vertical_regrid(input_press, input_values, output_press):
             out_array[y,x,:] = f(xnew)
     return out_array
 
+def mod_to_overpasstime(modobj,opass_tms):
+    '''
+    Interpolate model to satellite overpass time.
+
+    Parameters
+    ----------
+    modobj : xarray, model data
+    opass_tms : datetime64, satellite overpass local time
+
+    Output
+    ------
+    outmod : revised model data at local overpass time
+    '''
+    import pandas as pd
+    import xarray as xr
+
+    nst, = opass_tms.shape
+    nmt, = modobj.time.shape
+    ny,nx = modobj.longitude.shape
+
+    # Determine local time offset
+    local_utc_offset = np.zeros([ny,nx],dtype='timedelta64[ns]')
+    # pandas timedelta calculation doesn't work on ndarrays
+    for xi in np.arange(nx):
+        local_utc_offset[:,xi] = pd.to_timedelta((modobj['longitude'].isel(x=xi)/15).astype(np.int64),unit='h')
+
+    # initialize local time as variable
+    modobj['localtime'] = (['time','y','x'],np.zeros([nmt,ny,nx],dtype='datetime64[ns]'))
+    # fill
+    for ti in np.arange(nmt):
+        modobj['localtime'][ti] = modobj['time'][ti].data + local_utc_offset
+
+    # initalize new model object with satellite datetimes
+    outmod = []
+
+    for ti in np.arange(nst):
+        # Apply filter to select model data within +/- 1 output time step of the overpass time
+        tempmod = modobj.where(np.abs(modobj['localtime'] - opass_tms[ti].to_datetime64()) < (modobj.time[1] - modobj.time[0]))
+
+        # determine factors for linear interpolation in time
+        tfac = 1 - (np.abs(tempmod['localtime'] - opass_tms[ti].to_datetime64())/(modobj.time[1] - modobj.time[0]))
+        tempmod = tempmod.drop_vars('localtime')
+        # Carry out time interpolation
+        ## Note regarding current behavior: will only carry out time interpolation if at least 2 model timesteps
+        outmod.append((tfac*tempmod).sum(dim='time', min_count=2,keep_attrs=True))
+    #print(outmod)
+    outmod = xr.concat(outmod,dim='time')
+    outmod['time'] = (['time'],opass_tms)
+    return outmod
+
 def check_timestep(model_data,obs_data):
     ''' When pairing to level 3 data, model data may need to be aggregated to observation timestep.
         This function checks if the model data and observation data have the same timestep. Model data 
@@ -68,17 +118,31 @@ def mopitt_l3_pairing(model_data,obs_data,co_ppbv_varname):
     ## Check if same number of timesteps:
     if obs_data.time.size == model_data.time.size:
         model_obstime = model_data
-    elif obs_data.time.size < model_data.time.size and obs_data.time.size >2:
-        model_obstime = check_timestep(model_data,obs_data)
-    elif obs_data.time.size < model_data.time.size and obs_data.time.size < 3:
-        print('Model data and obs data timesteps do not match, and there are not enough observation timesteps to infer the spacing from.')
-        raise
-    elif obs_data.time.size > mod_data.time.size:
+        filtstr = '%Y-%m'
+    elif obs_data.time.size > model_data.time.size:
         print('Observation data appears to be a finer time resolution than model data')
         raise
+    elif obs_data.attrs['monthly']:
+        # if obs_data is monthly, take montly mean of model data
+        model_obstime = model_data.resample(time='MS').mean()
+        filtstr = '%Y-%m'
+    elif obs_data.attrs['monthly'] == False:
+        # obs_data is daily, so model and obs seem to be on same time step
+        model_obstime = model_data
+        filtstr = '%Y-%m-%d'
+    else:
+        # check frequency of model data 
+        tstep = xr.infer_freq(model_data.time.dt.round('D'))
+        if tstep == 'MS' or tstep == 'M':
+            model_obstime = model_data
+            filtstr = '%Y-%m'
+        else:
+            print('Time resolution of model data and MOPITT data is incompatible')
+            raise
+
     # initialize regridder for horizontal interpolation 
     # from model grid to MOPITT grid
-    grid_adjust = xe.Regridder(model_obstime[['latitude','longitude']],obs_data[['lat','lon']],'bilinear')
+    grid_adjust = xe.Regridder(model_obstime[['latitude','longitude']],obs_data[['lat','lon']],'bilinear',periodic=True)
     co_model_regrid = grid_adjust(model_obstime[co_ppbv_varname])
     pressure_model_regrid = grid_adjust(model_obstime['pres_pa_mid']/100.)
 
@@ -90,7 +154,8 @@ def mopitt_l3_pairing(model_data,obs_data,co_ppbv_varname):
     co_regrid = xr.full_like(obs_data['pressure'], np.nan)
     # MEB: loop over time outside of regrid lowers memory usage
     for t in range(obs_data.time.size):
-        co_regrid[t] = vertical_regrid(pressure_model_regrid[t].values, co_model_regrid[t].values, obs_data['pressure'][t].values)
+        obs_day = obs_data.time[t].dt.strftime(filtstr)
+        co_regrid[t] = vertical_regrid(pressure_model_regrid.sel(time=obs_day).values.squeeze(), co_model_regrid.sel(time=obs_day).values.squeeze(), obs_data['pressure'][t].values)
 
     # apply AK
     ## log apriori and model data
@@ -137,7 +202,7 @@ def omps_l3_daily_o3_pairing(model_data,obs_data,ozone_ppbv_varname):
     column = (du_fac*(model_data['dp_pa']/100.)*model_data[ozone_ppbv_varname]).sum('z')
 
     # initialize regrid and apply to column data
-    grid_adjust = xe.Regridder(model_data[['latitude','longitude']],obs_data[['latitude','longitude']],'bilinear')
+    grid_adjust = xe.Regridder(model_data[['latitude','longitude']],obs_data[['latitude','longitude']],'bilinear',periodic=True)
     mod_col_obsgrid = grid_adjust(column)
     # Aggregate time-step to daily means
     daily_mean = mod_col_obsgrid.resample(time='1D').mean()