from logging import info, debug
import xarray as xr
import pandas as pd
import numpy as np
import datetime
import netCDF4 as nc
from ....utils.datastores.empty import init_empty
[docs]
def do_parse(
self,
obs_file,
**kwargs
):
"""Parse function for a file from template observations
Args:
obs_file (str) :
Path to input file
Returns:
pandas.DataFrame :
Dataframe from input file df[parameter][station]
"""
info("Parsing official TROPOMI data "
"from observation file: {}".format(obs_file))
ddi, ddf = self.datei, self.datef
product_data = xr.open_dataset(obs_file)
obs_dates = product_data["time_utc"].values.astype(
np.datetime64).astype(datetime.datetime)
info(obs_dates.size)
if obs_dates.size :
if ddi > obs_dates.max() or ddf < obs_dates.min():
debug(f"File {obs_file} with date {obs_dates.min()} "
f"is outside simulation time period {ddi} / {ddf}")
return init_empty()
else :
debug(f"Empty {obs_file}")
return init_empty()
self.data_date_max = obs_dates.max()
# Filter out data outside domain
longitude = product_data.variables['longitude'].values
latitude = product_data.variables['latitude'].values
longitude_bounds = product_data.variables["longitude_bounds"].values
latitude_bounds = product_data.variables["latitude_bounds"].values
# better to use corners to crop domain ...
mask_domain = np.ones_like(longitude, dtype=bool)
if hasattr(self, "crop_domain"):
crop_domain = self.crop_domain
mask_domain = ((longitude < crop_domain.xmax) &
(longitude > crop_domain.xmin) &
(latitude < crop_domain.ymax) &
(latitude > crop_domain.ymin))
if mask_domain.sum() == 0:
return init_empty()
# Check number of levels
nlevread = product_data.dims["phony_dim_2"]
if nlevread != self.nlayers:
raise Exception(
f"The file {obs_file} does not have the expected number of "
f"averaging kernel layers: {nlevread} instead of "
f"{self.nlayers}. \n"
f"Please check your files or change the parameter 'nlayers' "
f"in the Yaml (expert users only).")
# Valid data
scanline = product_data.dims["phony_dim_0"]
ground_pixel = product_data.dims["phony_dim_1"]
obs = product_data.variables['nitrogendioxide_tropospheric_column'].values
conv = product_data.variables[
'nitrogendioxide_tropospheric_column'
].attrs['multiplication_factor_to_convert_to_molecules_percm2']
obs*=conv
obserror = product_data.variables[
'nitrogendioxide_tropospheric_column_precision'].values
conv = product_data.variables[
'nitrogendioxide_tropospheric_column_precision'].attrs['multiplication_factor_to_convert_to_molecules_percm2']
obserror*=conv
mask_valid = product_data.variables['qa_value'].values >= self.qa_value
#mask_valid = product_data.variables['qa_value'].values >= self.qa_value
mask_valid = mask_valid & ~np.isnan(obs)
#mask_valid = ~np.isnan(obs)
if (mask_domain & mask_valid).sum() == 0:
return init_empty()
# Filter out data not in simulation window
mask_time = (obs_dates >= ddi) & (obs_dates <= ddf)
mask_all = mask_time[..., np.newaxis] & mask_valid & mask_domain
if mask_all.sum() == 0:
return init_empty()
# Now process actual data
debug(f"Data found in file {obs_file}")
list_basic_cols = [
'date', 'duration', 'station', 'network', 'parameter',
'lon', 'lat', 'obs', 'obserror']
subdata = pd.DataFrame(columns=list_basic_cols)
subdata['date'] = np.repeat(
obs_dates[..., np.newaxis], ground_pixel, axis=1)[mask_all]
subdata = subdata.assign(duration=self.obs_duration)
subdata['lon'] = longitude[mask_all]
subdata['lat'] = latitude[mask_all]
subdata = subdata.assign(station=self.stationID)
subdata = subdata.assign(network=self.networkID)
subdata = subdata.assign(parameter=self.parameter)
subdata["obs"] = obs[mask_all].astype(float)
subdata["obserror"] = obserror[mask_all].astype(float)
# Load aks and obs
detailed_results_grp = product_data
aktot = np.float64(
detailed_results_grp.variables['averaging_kernel']
.values[mask_all, :]
)
# Pressure levels
# particular case of TROPOM CH4: create pressure grid from psurf and delta
psurf = product_data.variables['surface_pressure'].values
a = product_data.variables['tm5_constant_a'].values[:,1].T
b = product_data.variables['tm5_constant_b'].values[:,1].T
plev = np.array([a[lev] + b[lev]*psurf for lev in range(self.nlayers)])
pavg = np.concatenate((psurf[np.newaxis,:],plev),axis=0)
pavg = np.float64(pavg[:,mask_all]).T
# Load aks
trop_id = product_data.variables['tm5_tropopause_layer_index'].values
aktot = product_data.variables['averaging_kernel'].values
amftot = product_data.variables['air_mass_factor_total'].values
amftrop = product_data.variables['air_mass_factor_troposphere'].values
aktrop = np.array([ aktot[:,:,lev]*amftot[:,:]/amftrop[:,:] for lev in range(self.nlayers)])
lat_bounds = latitude_bounds[mask_all,:]
lon_bounds = longitude_bounds[mask_all,:]
if False:
for index, idtrop in np.ndenumerate(trop_id):
if np.isnan(idtrop): continue
aktrop[int(idtrop):,index[0],index[1]]=0
aktrop = np.float64(aktrop[:,mask_all]).T
ds = xr.Dataset({'ak': (['index', 'level'], aktrop),
'pavg0': (['index', 'level_pressure'], pavg),
'lat0': (['index'],lat_bounds[:,0]),
'lat1': (['index'],lat_bounds[:,1]),
'lat2': (['index'],lat_bounds[:,2]),
'lat3': (['index'],lat_bounds[:,3]),
'lon0': (['index'],lon_bounds[:,0]),
'lon1': (['index'],lon_bounds[:,1]),
'lon2': (['index'],lon_bounds[:,2]),
'lon3': (['index'],lon_bounds[:,3])},
coords={'index': subdata.index,
'level': range(self.nlayers),
'level_pressure': range(self.nlayers + 1)})
subdata = subdata.to_xarray()
subdata = xr.merge([ds, subdata])
debug(f"Observations after filter: {len(subdata.index)}")
return subdata
