diff --git a/conda_recipe/conda_build_config.yaml b/conda_recipe/conda_build_config.yaml index 7f18f5e..7bdf815 100644 --- a/conda_recipe/conda_build_config.yaml +++ b/conda_recipe/conda_build_config.yaml @@ -1,4 +1,3 @@ python: - - 3.5 - 3.6 - 3.7 \ No newline at end of file diff --git a/conda_recipe/meta.yaml b/conda_recipe/meta.yaml index 23926c8..b0b38a2 100644 --- a/conda_recipe/meta.yaml +++ b/conda_recipe/meta.yaml @@ -25,3 +25,4 @@ requirements: - metpy - basemap-data-hires - tqdm + - bottleneck diff --git a/eurec4a_snd/L1_bufr.py b/eurec4a_snd/L1_bufr.py index 9de1525..f4db29d 100644 --- a/eurec4a_snd/L1_bufr.py +++ b/eurec4a_snd/L1_bufr.py @@ -22,12 +22,15 @@ import numpy as np import netCDF4 from netCDF4 import Dataset, default_fillvals, num2date, date2num +import xarray as xr sys.path.append(os.path.dirname(os.path.abspath(__file__))) from config import cfg_creator as configupdater from _helpers import * import _thermo as thermo +json_config_fn = '/'.join([os.path.dirname(os.path.abspath(__file__)),'config/mwx_config.json']) + def get_args(): parser = argparse.ArgumentParser(formatter_class=argparse.RawTextHelpFormatter) @@ -78,6 +81,24 @@ def get_args(): 'not be included in the converted files.', required=False, default=[], nargs='+', type=int) + parser.add_argument("--instrument_id", metavar='Instrument identifier', + help="Instrument identifier e.g. Vaisala-RS or Meteomodem-RS", + default='RS', + required=False, + type=str) + + parser.add_argument("--platform_id", metavar='Platform identifier', + help="Platform identifier as used in config e.g. Atalante or BCO", + default=None, + required=False, + type=str) + + parser.add_argument("--campaign", metavar='Campaign', + help="Campaign name as used in config e.g. EUREC4A", + default='EUREC4A', + required=False, + type=str) + parser.add_argument('-v', '--verbose', metavar="DEBUG", help='Set the level of verbosity [DEBUG, INFO,' ' WARNING, ERROR]', @@ -119,7 +140,7 @@ def main(args={}): try: git_module_version = sp.check_output( - ["git", "describe", "--always", "--dirty"], stderr=sp.STDOUT).strip() + ["git", "describe", "--always", "--dirty"], stderr=sp.STDOUT).strip().decode() git_version_set = True except (sp.CalledProcessError, FileNotFoundError): logging.debug('No git-version could be found.') @@ -163,7 +184,8 @@ def main(args={}): logging.debug('Create filelist') if args['inputfile'] is None: - filelist = args['inputpath'].glob('*.bfr') + filelist = list(args['inputpath'].glob('*.bfr')) + filelist.extend(list(args['inputpath'].glob('*.BFR'))) else: filelist = [args['inputfile']] filelist = sorted(filelist) @@ -207,25 +229,17 @@ def main(args={}): sounding.direction = get_sounding_direction(sounding.meta_data['bufr_msg']) if sounding.direction == 1: # Upward - direction_str = 'Ascent' + direction_str = 'ascent' elif sounding.direction == -1: # Downward - direction_str = 'Descent' + direction_str = 'descent' sounding = expected_unit_check(sounding) - # after all needed header information is read, the reduced data field - # is masked for NaN values and an output file produced afterward: - logging.debug('Mask invalid sounding data') - sounding.time = np.ma.masked_invalid(sounding.time) - sounding.gpm = np.ma.masked_invalid(sounding.gpm) - sounding.pressure = np.ma.masked_invalid(sounding.pressure) - sounding.temperature = np.ma.masked_invalid(sounding.temperature) - sounding.dewpoint = np.ma.masked_invalid(sounding.dewpoint) - sounding.windspeed = np.ma.masked_invalid(sounding.windspeed) - sounding.winddirection = np.ma.masked_invalid(sounding.winddirection) - sounding.latitude = np.ma.masked_invalid(sounding.latitude) - sounding.longitude = np.ma.masked_invalid(sounding.longitude) + # Correct surface values of MeteoModem soundings + if sondetype == 177: + sounding = correct_meteomodem_surface(sounding, bufr_file) + # Calculate additional variables e = thermo.es(sounding.dewpoint, sounding.pressure*100) @@ -242,7 +256,7 @@ def main(args={}): sounding = exclude_specific_extendedVerticalSoundingSignificance_levels(sounding, args['significant_levels']) # Remove 1000hPa reduced gpm - sounding = exclude_1000hPa_gpm(sounding) + #sounding = exclude_1000hPa_gpm(sounding) # Ascent rate sounding = calc_ascentrate(sounding) @@ -251,7 +265,19 @@ def main(args={}): sounding = sort_sounding_by_time(sounding) # Find temporal resolution - time_resolution = calc_temporal_resolution(sounding) + resolution = calc_temporal_resolution(sounding) + + # Get global attributes + campaign = args['campaign'] + platform_id = args['platform_id'] + instrument_id = args['instrument_id'] + level = 'L1' + if package_version_set: + version = 'v{}'.format(__version__) + else: + version = git_module_version + glob_attrs_dict = get_global_attrs(json_config_fn, f'{campaign}_{platform_id}_{instrument_id}_{level}') + platform_location = glob_attrs_dict['platform_location'] # Create outputfile with time information from file sounding_date = sounding.sounding_start_time @@ -261,226 +287,136 @@ def main(args={}): outpath = Path(sounding.sounding_start_time.strftime(outpath_fmt.as_posix())) if outpath.suffix == '.nc': - outfile = Path(outpath.as_posix().format(platform=config['PLATFORM']['platform_name_short'], - location=config['PLATFORM']['platform_location']. - replace(' ', ''). - replace(',', ''). - replace(';', ''), - direction='{}Profile'.format(direction_str), - date=sounding_date.strftime('%Y%m%d_%H%M'))) + outfile = Path(outpath.as_posix().format(campaign=campaign, + platform=platform_id, + instrument=instrument_id, + level=level, + direction='{}'.format(direction_str), + date=sounding_date.strftime('%Y%m%dT%H%M'), + version=version)) else: outfile = Path(os.path.join(outpath, \ - "{platform}_Sounding{direction}_{location}_{date}.nc".\ - format(platform=config['PLATFORM']['platform_name_short'], - location=config['PLATFORM']['platform_location']. - replace(' ', ''). - replace(',', ''). - replace(';', ''), - direction='{}Profile'.format(direction_str), - date=sounding_date.strftime('%Y%m%d_%H%M')))) + "{campaign}_{platform}_{instrument}_{level}-{direction}_{date}_{version}.nc".\ + format(campaign=campaign, + platform=platform_id, + instrument=instrument_id, + level=level, + direction='{}'.format(direction_str), + date=sounding_date.strftime('%Y%m%dT%H%M'), + version=version))) if not outfile.parent.exists(): os.makedirs(outfile.parent) - # Creation of output NetCDF file - fo = Dataset(outfile, 'w', format='NETCDF4') - - # assign NetCDF file attributes from meta data - fo.title = 'Sounding data containing temperature, pressure, humidity,'\ - ' latitude, longitude, wind direction, wind speed, and time' - # Platform information - fo.platform_name = '{long} ({short})'.format( - long=config['PLATFORM']['platform_name_long'], - short=config['PLATFORM']['platform_name_short']) - fo.location = config['PLATFORM']['platform_location'] - fo.surface_altitude = config['PLATFORM']['platform_altitude'] - - # Instrument metadata - fo.instrument = config['INSTRUMENT']['instrument_description'] - fo.number_of_Probe = serial - fo.sonde_type = str(sondetype) - fo.sonde_frequency = sondefreq - - # Information about launch - fo.date_YYYYMMDD = sounding_date.strftime('%Y%m%d') - fo.time_of_launch_HHmmss = sounding_date.strftime('%H%M%S') - fo.launch_unixtime = date2num(sounding.sounding_start_time, date_unit) - fo.latitude_of_launch_location = '{0:5.2f} deg N'.\ - format(sounding.station_lat) - fo.longitude_of_launch_location = '{0:6.2f} deg E'.\ - format(sounding.station_lon) - - # Information about output - fo.resolution = "{:g} sec".format(time_resolution) - fo.source = Path(PureWindowsPath(bufr_file)).absolute().as_posix() - fo.git_version = git_module_version - fo.created_with = '{file} with its last modifications on {time}'.\ - format(time=time.ctime(os.path.getmtime(os.path.realpath(__file__))), - file=os.path.basename(__file__)) - fo.created_on = str(time.ctime(time.time())) - fo.contact_person = '{name} ({mail})'.format( + xr_output = sounding.to_dataset() + ## Global + xr_output.attrs['title'] = "EUREC4A level 1 sounding data".format(campaign) + xr_output.attrs['campaign_id'] = campaign + xr_output.attrs['platform_id'] = f'{platform_id}' + xr_output.attrs['instrument_id'] = f'{instrument_id}' + xr_output.attrs['platform_location'] = platform_location + xr_output.attrs['contact_person'] = '{name} ({mail})'.format( name=config['OUTPUT']['contact_person_name'], mail=config['OUTPUT']['contact_person_email']) - fo.institution = config['OUTPUT']['institution'] - fo.converted_by = '{name} ({mail})'.format( + xr_output.attrs['institution'] = config['OUTPUT']['institution'] + xr_output.attrs['converted_by'] = '{name} ({mail})'.format( name=config['OUTPUT']['executive_person_name'], mail=config['OUTPUT']['executive_person_email']) - fo.python_version = "{} (with numpy:{}, netCDF4:{}, eurec4a_snd:{})".\ - format(sys.version, np.__version__, netCDF4.__version__, - __version__) - fo.Conventions = 'CF-1.7' - fo.featureType = "trajectory" - - # Define Dimension (record length) from ASCII record counter - fo.createDimension('levels', len(sounding.pressure)) - fo.createDimension('sounding', None) - fo.createDimension('str_dim', 1000) - fillval = default_fillvals['f4'] - - # Creation of NetCDF Variables, including description and unit - nc_prof = fo.createVariable( - 'sounding', 'S1', ('sounding', 'str_dim'), - fill_value='', - zlib=True) - nc_prof.cf_role = "sounding_id" - nc_prof.long_name = 'sounding identifier' - nc_prof.description = 'unique string describing the soundings origin' - - nc_launchtime = fo.createVariable('launch_time', 'f8', ('sounding'), - zlib=True) - nc_launchtime.long_name = "time at which the sonde has been launched" - nc_launchtime.units = 'seconds since 1970-01-01 00:00:00 UTC' - nc_launchtime.calendar = 'gregorian' - nc_launchtime.standard_name = 'time' - - nc_tindex = fo.createVariable( - 'flight_time', 'f4', ('sounding', 'levels'), - fill_value=fillval, - zlib=True) - nc_tindex.long_name = 'time passed since launch' - nc_tindex.standard_name = 'time' - nc_tindex.units = 'seconds since {launch}'.format( - launch=sounding_date.strftime('%Y-%m-%d %H:%M:%S UTC')) - nc_tindex.axis = 'T' - nc_tindex.calendar = "gregorian" - nc_vvert = fo.createVariable( - 'ascentRate', 'f4', ('sounding', 'levels'), - fill_value=fillval, - zlib=True) - nc_vvert.long_name = 'ascent/descent rate of balloon or other measuring device' - nc_vvert.description = 'ascent rate is positive/ descent rate is negative' - nc_vvert.units = 'm/s' - nc_vvert.coordinates = "flight_time longitude latitude pressure" - nc_alti = fo.createVariable( - 'altitude', 'f4', ('sounding', 'levels'), - fill_value=fillval, - zlib=True) - nc_alti.standard_name = 'geopotential_height' - nc_alti.units = 'gpm' - nc_alti.coordinates = "flight_time longitude latitude pressure" - nc_pres = fo.createVariable( - 'pressure', 'f4', ('sounding', 'levels'), - fill_value=fillval, - zlib=True) - nc_pres.standard_name = 'air_pressure' - nc_pres.units = 'hPa' - nc_pres.axis = 'Z' - nc_pres.positive = 'down' - nc_temp = fo.createVariable( - 'temperature', 'f4', ('sounding', 'levels'), - fill_value=fillval, - zlib=True) - nc_temp.standard_name = 'air_temperature' - nc_temp.units = 'degrees_Celsius' - nc_temp.coordinates = "flight_time longitude latitude pressure" - nc_rh = fo.createVariable( - 'humidity', 'f4', ('sounding', 'levels'), - fill_value=fillval, - zlib=True) - nc_rh.standard_name = 'relative_humidity' - nc_rh.units = '%' - nc_rh.coordinates = "flight_time longitude latitude pressure" - nc_dewp = fo.createVariable( - 'dewPoint', 'f4', ('sounding', 'levels'), - fill_value=fillval, - zlib=True) - nc_dewp.standard_name = 'dew_point_temperature' - nc_dewp.units = 'degrees_Celsius' - nc_dewp.coordinates = "flight_time longitude latitude pressure" - nc_mix = fo.createVariable( - 'mixingRatio', 'f4', ('sounding', 'levels'), - fill_value=fillval, - zlib=True) - nc_mix.long_name = 'water vapor mixing ratio' - nc_mix.standard_name = 'humidity_mixing_ratio' - nc_mix.units = 'g/kg' - nc_mix.coordinates = "flight_time longitude latitude pressure" - nc_vhori = fo.createVariable( - 'windSpeed', 'f4', ('sounding', 'levels'), - fill_value=fillval, - zlib=True) - nc_vhori.standard_name = 'wind_speed' - nc_vhori.units = 'm/s' - nc_vhori.coordinates = "flight_time longitude latitude pressure" - nc_vdir = fo.createVariable( - 'windDirection', 'f4', ('sounding', 'levels'), - fill_value=fillval, - zlib=True) - nc_vdir.standard_name = 'wind_from_direction' - nc_vdir.units = 'degrees' - nc_vdir.coordinates = "flight_time longitude latitude pressure" - nc_lat = fo.createVariable( - 'latitude', 'f4', ('sounding', 'levels'), - fill_value=fillval, - zlib=True) - nc_lat.long_name = 'latitude' - nc_lat.standard_name = 'latitude' - nc_lat.units = 'degrees_north' - nc_lat.axis = 'Y' - nc_long = fo.createVariable( - 'longitude', 'f4', ('sounding', 'levels'), - fill_value=fillval, - zlib=True) - nc_long.long_name = 'longitude' - nc_long.standard_name = 'longitude' - nc_long.units = 'degrees_east' - nc_long.axis = 'X' + if 'surface_altitude' in glob_attrs_dict.keys(): + altitude = glob_attrs_dict['surface_altitude'] + else: + altitude = config['PLATFORM']['platform_altitude'] + xr_output.attrs['surface_altitude'] = altitude + if sondetype == 177: + xr_output.attrs['instrument'] = f'Radiosonde GPSonde M10 by MeteoModem' + elif sondetype == 123: + xr_output.attrs['instrument'] = f'Radiosonde RS41-SGP by Vaisala' + xr_output.attrs['number_of_probe'] = serial + xr_output.attrs['sonde_type'] = str(sondetype) + xr_output.attrs['sonde_frequency'] = sondefreq + xr_output.attrs['date_YYYYMMDD'] = sounding_date.strftime('%Y%m%d') + xr_output.attrs['time_of_launch_HHmmss'] = sounding_date.strftime('%H%M%S') + xr_output.attrs['launch_unixtime'] = float(xr_output['launch_time'].values[0]) + xr_output.attrs[ + 'latitude_of_launch_location'] = '{0:5.2f} deg N'.\ + format(sounding.station_lat) + xr_output.attrs[ + 'longitude_of_launch_location'] = '{0:6.2f} deg E'.\ + format(sounding.station_lon) + xr_output.attrs['source'] = str(Path(PureWindowsPath(bufr_file))) + xr_output.attrs['git_version'] = git_module_version + xr_output.attrs['created_with'] = '{file} with its last modifications on {time}'. \ + format(time=time.ctime(os.path.getmtime(os.path.realpath(__file__))), + file=os.path.basename(__file__)) + xr_output.attrs['created_on'] = str(time.ctime(time.time())) + xr_output.attrs['python_version'] = "{} (with numpy:{}, netCDF4:{}, eurec4a_snd:{})". \ + format(sys.version, np.__version__, netCDF4.__version__, __version__) + xr_output.attrs['resolution'] = f'{resolution} sec' + xr_output.attrs['Conventions'] = "CF-1.7" + xr_output.attrs['featureType'] = "trajectory" + + # Overwrite standard attrs with those defined in config file + # Get global meta data from mwx_config.json + glob_attrs_dict = get_global_attrs(json_config_fn, f'{campaign}_{platform_id}_{instrument_id}_{level}') + for attrs, value in glob_attrs_dict.items(): + xr_output.attrs[attrs] = value + if 'extendedVerticalSoundingSignificance' in args['additional_variables']: - nc_evss = fo.createVariable( - 'extendedVerticalSoundingSignificance', 'i4', ('sounding', 'levels'), - fill_value=0, - zlib=True) - nc_evss.long_name = 'extended vertical soudning significance' - nc_evss.description = 'see BUFR code flag table to decode' - nc_evss.units = '-' - - sounding_name = '{platform}__{lat:5.2f}_{lon:5.2f}__{launchtime}'.\ - format(platform=config['PLATFORM']['platform_name_short'], + xr_output['extendedVerticalSoundingSignificance'] = xr.DataArray([sounding.extendedVerticalSoundingSignificance], dims=['sounding','levels']) + + sounding_name = '{platform}__{direction}__{lat:05.2f}_{lon:06.2f}__{launchtime}'.\ + format(platform=platform_id, + direction=direction_str, lat=sounding.station_lat, lon=sounding.station_lon, launchtime=str(YYYYMMDDHHMM)) - sounding_name_parts = [] - for char in sounding_name: - sounding_name_parts.extend(char) - - nc_prof[0, 0:len(sounding_name_parts)] = sounding_name_parts - nc_launchtime[0] = date2num(sounding.sounding_start_time, date_unit) - - nc_tindex[0, :] = sounding.time - nc_vvert[0, :] = sounding.ascentrate - nc_alti[0, :] = sounding.gpm - nc_pres[0, :] = sounding.pressure - nc_temp[0, :] = sounding.temperature - nc_rh[0, :] = sounding.relativehumidity - nc_dewp[0, :] = sounding.dewpoint - nc_mix[0, :] = sounding.mixingratio - nc_vhori[0, :] = sounding.windspeed - nc_vdir[0, :] = sounding.winddirection - nc_lat[0, :] = sounding.latitude - nc_long[0, :] = sounding.longitude - if 'extendedVerticalSoundingSignificance' in args['additional_variables']: - nc_evss[0, :] = sounding.extendedVerticalSoundingSignificance - fo.close() + + xr_output['sounding_id'] = xr.DataArray([sounding_name], dims = ['sounding']) + + with open(json_config_fn, 'r') as f: + j = json.load(f) + meta_data_dict = j['meta_data'] + for variable in xr_output.data_vars: + if variable in meta_data_dict.keys(): + variable_meta_data = meta_data_dict[variable] + for attr, value in variable_meta_data.items(): + xr_output[variable].attrs[attr] = value + + # Reduce dtype to float instead of double + xr_output.sounding_id.encoding = {'dtype': 'S1000', 'char_dim_name': 'str_dim'} + for variable in ['altitude', 'ascentRate', 'dewPoint', 'humidity', 'latitude', 'longitude', + 'mixingRatio', 'pressure', 'temperature', 'windDirection', 'windSpeed']: + xr_output[variable].encoding['dtype'] = 'f4' + xr_output['flight_time'].encoding['dtype'] = 'f8' + + for variable in xr_output.data_vars: + xr_output[variable].encoding['zlib'] = True + + xr_output.to_netcdf(outfile, unlimited_dims=['sounding']) + + # sounding_name_parts = [] + # for char in sounding_name: + # sounding_name_parts.extend(char) + + # nc_prof[0, 0:len(sounding_name_parts)] = sounding_name_parts + # nc_launchtime[0] = date2num(sounding.sounding_start_time, date_unit) + # + # nc_tindex[0, :] = sounding.time + # nc_vvert[0, :] = sounding.ascentrate + # nc_alti[0, :] = sounding.gpm + # nc_pres[0, :] = sounding.pressure + # nc_temp[0, :] = sounding.temperature + # nc_rh[0, :] = sounding.relativehumidity + # nc_dewp[0, :] = sounding.dewpoint + # nc_mix[0, :] = sounding.mixingratio + # nc_vhori[0, :] = sounding.windspeed + # nc_vdir[0, :] = sounding.winddirection + # nc_lat[0, :] = sounding.latitude + # nc_long[0, :] = sounding.longitude + # if 'extendedVerticalSoundingSignificance' in args['additional_variables']: + # nc_evss[0, :] = sounding.extendedVerticalSoundingSignificance + # fo.close() + logging.info('DONE: {input} converted to {output}'.format( input=filelist[ifile], output=outfile)) diff --git a/eurec4a_snd/L1_mwx.py b/eurec4a_snd/L1_mwx.py index c5a6e62..09b5070 100644 --- a/eurec4a_snd/L1_mwx.py +++ b/eurec4a_snd/L1_mwx.py @@ -9,13 +9,9 @@ import os from pathlib import Path, PureWindowsPath import subprocess as sp -import configparser -from configparser import ExtendedInterpolation import argparse import logging -import tempfile import glob -from xml.dom import minidom from datetime import timedelta import datetime as dt import netCDF4 @@ -24,19 +20,18 @@ import tqdm import json import numpy as np -import pandas as pd import xarray as xr import sys -sys.path.append('.') +sys.path.append(os.path.dirname(os.path.abspath(__file__))) from _mwx_helpers import * -from _helpers import calc_saturation_pressure, unixpath, setup_logging, load_configuration +from _helpers import calc_saturation_pressure, unixpath, setup_logging, load_configuration, get_global_attrs # User-configuration -campaign = 'EUREC4A' -output_filename_format = '{campaign}_{platform_short}_sounding_{direction}_%Y%m%d_%H%M.nc' # including path +level = 'L1' +output_filename_format = '{campaign}_{platform_short}_{instrument_id}_{level}-{direction}_%Y%m%dT%H%M_{version}.nc' # including path -json_config_fn = 'mwx_config.json' +json_config_fn = '/'.join([os.path.dirname(os.path.abspath(__file__)),'config/mwx_config.json']) def get_args(): parser = argparse.ArgumentParser(formatter_class=argparse.RawTextHelpFormatter) @@ -93,6 +88,21 @@ def get_args(): default=None, required=False, type=str) + parser.add_argument("--instrument_id", metavar='Instrument identifier', + help="Instrument identifier e.g. Vaisala-RS or Meteomodem-RS", + default='RS', + required=False, + type=str) + parser.add_argument("--platform_id", metavar='Platform identifier', + help="Platform identifier as used in config e.g. Atalante or BCO", + default=None, + required=False, + type=str) + parser.add_argument("--campaign", metavar='Campaign', + help="Campaign name as used in config e.g. EUREC4A", + default='EUREC4A', + required=False, + type=str) parser.add_argument("--round-like-BUFR", metavar='BOOLEAN', help="Switch on/off rounding of output values as output in BUFR files", default=False, @@ -172,7 +182,7 @@ def main(args={}): args["inputpath"] = config["FILES"]["INPUT_MWX"] if args["outputfolder"] is None: args["outputfolder"] = config["FILES"]["OUTPUT_MWX2NC"] - + output_folder = args['outputfolder'] basename = os.path.basename(output_folder) dirname = os.path.dirname(output_folder) @@ -193,12 +203,15 @@ def main(args={}): # Loading standard config file with open(json_config_fn, 'r') as f: j=json.load(f) - sync_sounding_values = j['sync_sounding_items'] std_sounding_values = j['std_sounding_items'] radiosondes_values = j['radiosondes_sounding_items'] meta_data_dict = j['meta_data'] + campaign = args['campaign'] + platform_id = args['platform_id'] + instrument_id = args['instrument_id'] + # Function definitions f_flighttime = lambda radiorx: begin_time_dt + timedelta(seconds=radiorx-float(launch_time)) @@ -226,74 +239,30 @@ def main(args={}): for mwx_file in tqdm.tqdm(mwx_files): # Decompress/open mwx file - tmpdir, tmpdir_obj = getTmpDir() - decompressed_files = np.array(decompress(mwx_file, tmpdir+'/')) - - # Get the files SynchronizedSoundingData.xml, Soundings.xml, ... - sync_mask = [f_sync(file) for file in decompressed_files] - if np.sum(sync_mask) == 0: - logging.warning('No sounding data found in {}. Skipped'.format(mwx_file)) - continue - sync_filename = decompressed_files[sync_mask][0] - snd_mask = [f_snd(file) for file in decompressed_files] - if np.sum(snd_mask) == 0: - logging.warning('No sounding data found in {}. Skipped'.format(mwx_file)) - continue - snd_filename = decompressed_files[snd_mask][0] - std_mask = [f_std(file) for file in decompressed_files] - if np.sum(std_mask) == 0: - logging.warning('No sounding data found in {}. Skipped'.format(mwx_file)) - continue - std_filename = decompressed_files[std_mask][0] - radio_mask = [f_radio(file) for file in decompressed_files] - radio_filename = decompressed_files[radio_mask][0] - - # Read Soundings.xml to get base time - itemlist = read_xml(snd_filename) - for i, item in enumerate(itemlist): - begin_time = item.attributes['BeginTime'].value - launch_time = item.attributes['LaunchTime'].value - altitude = item.attributes['Altitude'].value - begin_time_dt = dt.datetime.strptime(begin_time,'%Y-%m-%dT%H:%M:%S.%f') - - # Read SynchronizedSoundingData.xml with processed sounding data - itemlist = read_xml(sync_filename) - sounding_dict = {} - try: + with MWX(mwx_file) as mwx: + decompressed_files = mwx.decompressed_files + + # Get the files SynchronizedSoundingData.xml, Soundings.xml, ... + a1, sync_filename = check_availability(decompressed_files, 'SynchronizedSoundingData.xml', True) + a2, snd_filename = check_availability(decompressed_files, 'Soundings.xml', True) + a3, radio_filename = check_availability(decompressed_files, 'Radiosondes.xml', True) + if np.any([not a1, not a2, not a3]): + logging.warning('No sounding data found in {}. Skipped'.format(mwx_file)) + continue + + # Read Soundings.xml to get base time + itemlist = read_xml(snd_filename) for i, item in enumerate(itemlist): - level_dict = {} - for var in sync_sounding_values: - level_dict[var] = item.attributes[var].value - sounding_dict[i] = level_dict - except: - logging.warning('Key {} not found in file {}'.format(var, mwx_file)) - continue - pd_snd = pd.DataFrame.from_dict(sounding_dict, orient='index', dtype=float) - - # Read StdPressureLevels.xml to include measurements interpolated to std-levels - itemlist = read_xml(std_filename) - sounding_std_dict = {} - for i, item in enumerate(itemlist): - level_dict = {} - for var in std_sounding_values: - level_dict[var] = item.attributes[var].value - sounding_std_dict[i] = level_dict - pd_std = pd.DataFrame.from_dict(sounding_std_dict, orient='index', dtype=float) - - # Read Radiosounding.xml to get sounding metadata - itemlist = read_xml(radio_filename) - sounding_meta_dict = {} - for i, item in enumerate(itemlist): - assert i == 0, 'further entries were found, meaning soundings meta data could be mixed up' - for var in radiosondes_values: - try: - sounding_meta_dict[var] = item.attributes[var].value - except KeyError: - logging.error('Attribute {} could not found and is assumed to be RS41-SGP'.format(var)) - sounding_meta_dict[var] = 'RS41-SGP' - - # Convert missing values (-32768) to nan - pd_snd = pd_snd.replace(-32768, np.nan) + begin_time = item.attributes['BeginTime'].value + launch_time = item.attributes['LaunchTime'].value + altitude = item.attributes['Altitude'].value + begin_time_dt = dt.datetime.strptime(begin_time,'%Y-%m-%dT%H:%M:%S.%f') + + # Read sounding data + pd_snd = get_sounding_profile(sync_filename, sync_sounding_values) + + # Read Radiosounding.xml to get sounding metadata + sounding_meta_dict = get_sounding_metadata(radio_filename, radiosondes_values) # Create flight time pd_snd['flight_time'] = pd_snd.RadioRxTimePk.apply(f_flighttime) @@ -315,7 +284,7 @@ def main(args={}): # Write output for s,sounding in enumerate([pd_snd_asc, pd_snd_dsc]): if len(sounding) < 2: - logging.warning('Sounding ({}) does not contain data. Skip sounding-direction of{}'.format(direction_dict[s], mwx_file)) #, direction_dict[sounding.Dropping.values[0]])) + logging.warning('Sounding ({}) does not contain data. Skip sounding-direction of {}'.format(direction_dict[s], mwx_file)) #, direction_dict[sounding.Dropping.values[0]])) continue xr_snd = xr.Dataset.from_dataframe(sounding) @@ -325,7 +294,7 @@ def main(args={}): ## Dew point temperature dewpoint = convert_RH_to_dewpoint(xr_snd.Temperature.values, xr_snd.Humidity.values) ## Mixing ratio - e_s = calc_saturation_pressure(xr_snd.Temperature.values) # calc_vapor_pressure(xr_snd, method='hardy') + e_s = calc_saturation_pressure(xr_snd.Temperature.values) mixing_ratio = calc_wv_mixing_ratio(xr_snd, e_s)*xr_snd.Humidity.values/100. ## Launch time as type(datetime) flight_time_unix = sounding.flight_time.values.astype(float)/1e9 @@ -334,7 +303,8 @@ def main(args={}): ## Resolution resolution = calc_temporal_resolution(flight_time_unix) ## Sounding ID - sounding_id = '{platform}__{lat:3.2f}_{lon:4.2f}__{time}'.format(platform = platform_name_short, + sounding_id = '{platform}__{direction}__{lat:3.2f}_{lon:4.2f}__{time}'.format(platform = platform_id, + direction = direction_dict[s], lat = xr_snd.Latitude.values[0], lon = xr_snd.Longitude.values[0], time = launch_time_dt.strftime('%Y%m%d%H%M')) @@ -360,6 +330,7 @@ def main(args={}): xr_output['windDirection'] = xr.DataArray([xr_snd.WindDir.values], dims = ['sounding', 'levels']) xr_output['latitude'] = xr.DataArray([xr_snd.Latitude.values], dims = ['sounding', 'levels']) xr_output['longitude'] = xr.DataArray([xr_snd.Longitude.values], dims = ['sounding', 'levels']) + xr_output['altitude_WGS84'] = xr.DataArray([xr_snd.Altitude.values], dims=['sounding', 'levels']) # xr_output['standard_level_flag'] = xr.DataArray() # Write attributes @@ -371,9 +342,11 @@ def main(args={}): xr_output[variable].attrs[attr] = value ## Global - xr_output.attrs['title'] = "Sounding data during the EUREC4A campaign" - xr_output.attrs['platform_name'] = f'{platform_name_long} ({platform_name_short})' - xr_output.attrs['location'] = platform_location + xr_output.attrs['title'] = "Sounding data during the {} campaign (level 1)".format(campaign) + xr_output.attrs['campaign_id'] = campaign + xr_output.attrs['platform_id'] = f'{platform_id}' + xr_output.attrs['instrument_id'] = f'{instrument_id}' + xr_output.attrs['platform_location'] = platform_location xr_output.attrs['surface_altitude'] = '{:3.1f} m'.format(float(altitude)) xr_output.attrs['instrument'] = f'Radiosonde {sounding_meta_dict["SondeTypeName"]} by Vaisala' xr_output.attrs['number_of_probe'] = sounding_meta_dict['SerialNbr'] @@ -396,23 +369,36 @@ def main(args={}): xr_output.attrs['Conventions'] = "CF-1.7" xr_output.attrs['featureType'] = "trajectory" + # Overwrite standard attrs with those defined in config file + # Get global meta data from mwx_config.json + glob_attrs_dict = get_global_attrs(json_config_fn, f'{campaign}_{platform_id}_{instrument_id}_{level}') + for attrs, value in glob_attrs_dict.items(): + xr_output.attrs[attrs] = value + # Reduce dtype to float instead of double xr_output.sounding_id.encoding = {'dtype': 'S1000', 'char_dim_name': 'str_dim'} for variable in ['altitude', 'ascentRate', 'dewPoint', 'humidity', 'latitude', 'longitude', 'mixingRatio', 'pressure', 'temperature', 'windDirection', 'windSpeed']: - xr_output[variable].encoding = {'dtype': 'f4'} + xr_output[variable].encoding['dtype'] = 'f4' for variable in xr_output.data_vars: xr_output[variable].encoding['zlib'] = True + if package_version_set: + version = 'v{}'.format(__version__) + else: + version = git_module_version filename = output_format.format(campaign=campaign, - platform_short=platform_name_short, - direction=direction_dict[sounding.Dropping.values[0]]) + platform_short=platform_id, + direction=direction_dict[sounding.Dropping.values[0]], + instrument_id=args["instrument_id"], + version=version, + level=level + ) filename = launch_time_dt.strftime(filename) xr_output.to_netcdf(filename, unlimited_dims=['sounding']) logging.info('File converted to {}'.format(filename)) - tmpdir_obj.cleanup() if __name__ == "__main__": main() diff --git a/eurec4a_snd/_helpers.py b/eurec4a_snd/_helpers.py index 43144a5..c92d5a6 100644 --- a/eurec4a_snd/_helpers.py +++ b/eurec4a_snd/_helpers.py @@ -3,6 +3,7 @@ """ import tempfile import os +import re import inspect import platform import json @@ -10,6 +11,10 @@ import numpy as np import logging from pathlib import Path, PureWindowsPath +import configparser +from configparser import ExtendedInterpolation +from netCDF4 import date2num +import xarray as xr class UnitChangedError(Exception): @@ -19,6 +24,38 @@ class UnitChangedError(Exception): class UnexpectedUnit(Exception): pass +class RegexDict(dict): + """ + Dictionary with capability of taking regular xpressions + """ + def get_matching(self, event): + return (self[key] for key in self if re.match(key, event)) + + def get_matching_combined(self, event): + """ + Find matching keys and return combined dictionary + + >>> d = {'EUREC4A_*':{'a':0}, 'EUREC4A_BCO':{'p':1}} + >>> rd = RegexDict(d) + >>> rd.get_matching_combined('EUREC4A_BCO') + {'a': 0, 'p': 1} + """ + matching = self.get_matching(event) + dall = {} + for d in matching: + dall.update(d) + return dall + + +def get_global_attrs(cfg_file, key): + """ + Get global attributes from configuration file + """ + with open(cfg_file, 'r') as f: + j=json.load(f) + rd = RegexDict(j['global_meta_data']) + return rd.get_matching_combined(key) + def unixpath(path_in): """ @@ -188,6 +225,36 @@ def __init__(self): self.extendedVerticalSoundingSignificance = [] self.meta_data = {} + def to_dataset(self): + """ + Convert sounding to dataset + """ + date_unit = "seconds since 1970-01-01 00:00:00 UTC" + var_dict = { + 'launch_time': date2num(self.sounding_start_time, date_unit), + 'flight_time': self.time, + 'altitude': self.gpm, + 'ascentRate': self.ascentrate, + 'dewPoint': self.dewpoint, + 'humidity': self.relativehumidity, + 'latitude': self.latitude, + 'longitude': self.longitude, + 'mixingRatio': self.mixingratio, + 'pressure': self.pressure, + 'temperature': self.temperature, + 'windDirection': self.winddirection, + 'windSpeed': self.windspeed + } + + ds = xr.Dataset() + for var, value in var_dict.items(): + if isinstance(value, np.ndarray): + ds[var] = xr.DataArray([value], dims=['sounding', 'levels']) + else: + ds[var] = xr.DataArray([value], dims=['sounding']) + return ds + + def _ensure_measurement_integrity(self): """ Test integrity of each measurement unit @@ -345,7 +412,7 @@ def _ensure_measurement_integrity(self): msg_fmt = search_bufr_msg_format(descriptors) s.meta_data['bufr_msg'] = msg_fmt elif json_flat[key_key+'_key'] == 'radiosondeOperatingFrequency': - s.meta_data['sonde_frequency'] = str(json_flat[key_key+'_value']) + json_flat[key_key+'_units'] + s.meta_data['sonde_frequency'] = str(json_flat[key_key+'_value']) +' '+ json_flat[key_key+'_units'] _ensure_measurement_integrity(s) @@ -867,6 +934,7 @@ def exclude_sounding_level(sounding, nan_mask): return sounding + def exclude_specific_extendedVerticalSoundingSignificance_levels(sounding, significance_bits): """ Exclude levels with specific extendedVerticalSoundingSignificance @@ -896,7 +964,71 @@ def exclude_specific_extendedVerticalSoundingSignificance_levels(sounding, signi for significance_level in significance_levels: current_level = sounding.extendedVerticalSoundingSignificance[significance_level] current_level = set(decode_extendedVerticalSoundingSignificance(current_level)) - to_delete_mask[significance_level] = current_level.issubset(significance_bits) + #to_delete_mask[significance_level] = current_level.issubset(significance_bits) + to_delete_mask[significance_level]= ((current_level == set(significance_bits)) or (current_level == set([4]))) sounding = exclude_sounding_level(sounding, ~to_delete_mask) return sounding + + + +def correct_meteomodem_surface(sounding, bufr_file): + """ + Correct specific meteo modem soundings + """ + d = {'ATALANTE_2020012512_1.309057.BFR': 21.676335267277707, + 'ATALANTE_2020012612_1.309057.BFR': 20.299818275809564, + 'ATALANTE_2020012612_2.309057.BFR': 20.490235079280602, + 'ATALANTE_2020012618_1.309057.BFR': 20.395027551835543, + 'ATALANTE_2020012618_2.309057.BFR': 20.72719422108181, + 'ATALANTE_2020012618_3.309057.BFR': 20.72719422108181, + 'ATALANTE_2020012618_4.309057.BFR': 20.490235079280602, + 'ATALANTE_2020012618_5.309057.BFR': 20.961170412362186, + 'ATALANTE_2020012618_6.309057.BFR': 21.096516728956811, + 'ATALANTE_2020012700_1.309057.BFR': 21.228696671077643, + 'ATALANTE_2020012700_2.309057.BFR': 21.645982812790148, + 'ATALANTE_2020012700_3.309057.BFR': 21.676335267277707, + 'ATALANTE_2020012800_1.309057.BFR': 22.734527932378281, + 'ATALANTE_2020012800_2.309057.BFR': np.nan, + 'ATALANTE_2020020218_1.309057.BFR': 20.631811617269207, + 'ATALANTE_2020020218_2.309057.BFR': 21.192243399966955, + 'ATALANTE_2020020300_1.309057.BFR': 23.250129401644031, + 'ATALANTE_2020020300_2.309057.BFR': 20.770057518368084, + 'ATALANTE_2020020300_3.309057.BFR': 21.228696671077643, + 'ATALANTE_2020020300_4.309057.BFR': 21.132797893968846, + 'ATALANTE_2020020300_5.309057.BFR': 21.594661125167445, + 'ATALANTE_2020020300_6.309057.BFR': 22.347032338755625, + 'ATALANTE_2020020300_7.309057.BFR': 21.12715294391278, + 'ATALANTE_2020020306_1.309057.BFR': 21.16565351341211, + 'ATALANTE_2020020306_2.309057.BFR': 21.350899212666299, + 'ATALANTE_2020020306_3.309057.BFR': 21.114158685504709, + 'ATALANTE_2020020306_4.309057.BFR': 21.047255593956379, + 'ATALANTE_2020020306_5.309057.BFR': 20.653281422508336, + 'ATALANTE_2020020306_6.309057.BFR': 20.927958315808677, + 'ATALANTE_2020020306_7.309057.BFR': 21.357437545602338, + 'ATALANTE_2020020312_1.309057.BFR': 21.510617824089007, + 'ATALANTE_2020020312_2.309057.BFR': 21.195165907041719, + 'ATALANTE_2020020312_3.309057.BFR': np.nan, + 'ATALANTE_2020020312_4.309057.BFR': 21.000788492386096, + 'ATALANTE_2020020312_5.309057.BFR': 21.000788492386096, + 'ATALANTE_2020020312_6.309057.BFR': 21.13786737261098, + 'ATALANTE_2020020312_7.309057.BFR': 21.877213391027734, + 'ATALANTE_2020020312_8.309057.BFR': 21.469277211341332, + 'ATALANTE_2020020318_1.309057.BFR': 21.93987011060981, + 'ATALANTE_2020020318_3.309057.BFR': 21.543164982699349, + 'ATALANTE_2020020318_4.309057.BFR': 21.431461248468246, + 'ATALANTE_2020020318_5.309057.BFR': 21.31908609493534, + 'ATALANTE_2020021100_1.309057.BFR': 21.534170247830687, + 'ATALANTE_2020021800_1.309057.BFR': 21.216973566186809, + 'ATALANTE_2020021800_2.309057.BFR': 21.483871685146806, + 'ATALANTE_2020021800_4.309057.BFR': 22.123277349708609, + 'ATALANTE_2020021800_5.309057.BFR': 21.931339923645719} + if os.path.basename(bufr_file) in d.keys(): + Td_corrected = d[os.path.basename(bufr_file)] + print('Correct Td from {} to {}'.format(sounding.dewpoint[0], Td_corrected)) + sounding.dewpoint[0] = Td_corrected + if np.isnan(Td_corrected): + sounding.temperature[0] = np.nan + else: + print('Unknown BFR file. No correction will be applied.') + return sounding diff --git a/eurec4a_snd/_mwx_helpers.py b/eurec4a_snd/_mwx_helpers.py index b250464..0b042ce 100644 --- a/eurec4a_snd/_mwx_helpers.py +++ b/eurec4a_snd/_mwx_helpers.py @@ -5,10 +5,24 @@ import numpy as np import subprocess import shutil +import pandas as pd +import warnings class SondeTypeNotImplemented(Exception): pass +class VariableNotFoundInSounding(Warning): + pass + +class SondeTypeNotIdentifiable(Warning): + pass + +def custom_formatwarning(msg, *args, **kwargs): + # ignore everything except the message + return str(msg) + '\n' + +warnings.formatwarning = custom_formatwarning + def getTmpDir(): """ @@ -44,6 +58,84 @@ def compress(folder, compressed_file): return +def open_mwx(mwx_file): + """ + Open Vaisala MWX41 archive file (.mwx) + + Input + ----- + mwx_file : str + Vaisala MW41 archive file + + Returns + ------- + decompressed_files : list + List of temporarily decompressed .xml files + within the archive file + """ + tmpdir, tmpdir_obj = getTmpDir() + decompressed_files = np.array(decompress(mwx_file, tmpdir + '/')) + return decompressed_files + +class MWX(object): + """ + Open Vaisala MWX41 archive file (.mwx) + + Input + ----- + mwx_file : str + Vaisala MW41 archive file + + Returns + ------- + decompressed_files : list + List of temporarily decompressed .xml files + within the archive file + """ + def __init__(self, mwx_file): + self.tmpdir, self.tmpdir_obj = getTmpDir() + self.decompressed_files = np.array(decompress(mwx_file, self.tmpdir + '/')) + def __enter__(self): + return self + def __exit__(self, type, value, traceback): + self.tmpdir_obj.cleanup() + def get_decompressed_files(self): + return self.decompressed_files + + +def check_availability(decomp_files, file, return_name=False): + """ + Check whether xml file exist in decompressed + file list + + Returns + ------- + avail : bool + Availability of file + filename : str (optional) + Full filename of requested file + """ + basenames = [os.path.basename(decomp_file) for decomp_file in decomp_files] + + # Availability + availability_mask = np.in1d(basenames, file) + + if np.sum(availability_mask) > 0: + avail = True + else: + avail = False + + if return_name: + if avail: + idx = np.where(availability_mask)[0][0] + fullname = decomp_files[idx] + else: + fullname = None + return avail, fullname + else: + return avail + + def read_xml(filename, return_handle=False): xmldoc = minidom.parse(filename) itemlist = xmldoc.getElementsByTagName('Row') @@ -53,6 +145,53 @@ def read_xml(filename, return_handle=False): return itemlist +def get_sounding_profile(file, keys): + """ + Get sounding profile from provided xml file + + Input + ----- + file : str + XML file containing sounding data e.g. + SynchronizedSoundingData.xml + keys : list + list of variables to look for + + Returns + ------- + pd_snd : pandas.DataFrame + sounding profile + """ + itemlist = read_xml(file) + sounding_dict = {} + try: + for i, item in enumerate(itemlist): + level_dict = {} + for var in keys: + level_dict[var] = item.attributes[var].value + sounding_dict[i] = level_dict + except: + warnings.warn('Key {} not found.'.format(var), VariableNotFoundInSounding) + pd_snd = pd.DataFrame.from_dict(sounding_dict, orient='index', dtype=float) + + # Set missing values to NaN + pd_snd = pd_snd.replace(-32768, np.nan) + return pd_snd + + +def get_sounding_metadata(file, keys): + itemlist = read_xml(file, keys) + sounding_meta_dict = {} + for i, item in enumerate(itemlist): + assert i == 0, 'further entries were found, meaning soundings meta data could be mixed up' + for var in keys: + try: + sounding_meta_dict[var] = item.attributes[var].value + except KeyError: + warnings.warn('Attribute {} could not found and is assumed to be RS41-SGP'.format(var), SondeTypeNotIdentifiable) + sounding_meta_dict[var] = 'RS41-SGP' + return sounding_meta_dict + def calc_ascent_rate(sounding): """ Calculate ascent rate @@ -117,7 +256,7 @@ def calc_temporal_resolution(time): """ time_differences = np.abs(np.diff(np.ma.compressed(time))) time_differences_counts = np.bincount(time_differences.astype(int)) - most_common_diff = time_differences[np.argmax(time_differences_counts)] + most_common_diff = np.argmax(time_differences_counts) temporal_resolution = most_common_diff return temporal_resolution diff --git a/eurec4a_snd/mwx_config.json b/eurec4a_snd/config/mwx_config.json similarity index 53% rename from eurec4a_snd/mwx_config.json rename to eurec4a_snd/config/mwx_config.json index baa26fe..394e967 100644 --- a/eurec4a_snd/mwx_config.json +++ b/eurec4a_snd/config/mwx_config.json @@ -13,92 +13,91 @@ "units": "seconds since 1970-01-01 00:00:00 UTC", "axis": "T", "calendar": "gregorian", - "coordinates": "flight_time longitude latitude pressure", + "coordinates": "sounding_id flight_time longitude latitude pressure", "_FillValue": 9.969209968386869e+36 }, "launch_time": { + "standard_name": "time", "long_name": "time at which the sounding started", "units": "seconds since 1970-01-01 00:00:00 UTC", "calendar": "gregorian", - "standard_name": "time", "_FillValue": 9.969209968386869e+36 }, "ascentRate": { "long_name": "ascent/desent rate of measuring device", - "coordinates": "flight_time longitude latitude pressure", + "coordinates": "sounding_id flight_time longitude latitude pressure", "description": "ascent rate is positive/ descent rate is negative", "units": "m/s", "_FillValue": 9.969209968386869e+36 }, "pressure": { - "long_name": "pressure", - "coordinates": "flight_time longitude latitude pressure", "standard_name": "air_pressure", + "long_name": "pressure", "units": "hPa", - "axis": "Z", - "positive": "down", + "coordinates": "sounding_id flight_time longitude latitude", "_FillValue": 9.969209968386869e+36 }, "temperature": { - "long_name": "dry bulb temperature", - "coordinates": "flight_time longitude latitude pressure", "standard_name": "air_temperature", + "long_name": "dry bulb temperature", "units": "degrees_Celsius", + "coordinates": "sounding_id flight_time longitude latitude pressure", "_FillValue": 9.969209968386869e+36 }, "humidity": { - "long_name": "relative_humidity", "standard_name": "relative_humidity", - "coordinates": "flight_time longitude latitude pressure", + "long_name": "relative humidity", "units": "%", + "coordinates": "sounding_id flight_time longitude latitude pressure", "_FillValue": 9.969209968386869e+36 }, "specific_humidity": { - "long_name": "specific humidity", "standard_name": "specific_humidity", + "long_name": "specific humidity", "units": "g/kg", + "coordinates": "sounding_id flight_time longitude latitude pressure", "_FillValue": 9.969209968386869e+36 }, "dewPoint": { - "long_name": "dew point temperature", - "coordinates": "flight_time longitude latitude pressure", "standard_name": "dew_point_temperature", + "long_name": "dew point temperature", "units": "degrees_Celsius", + "coordinates": "sounding_id flight_time longitude latitude pressure", "_FillValue": 9.969209968386869e+36 }, "mixingRatio": { + "standard_name": "humidity_mixing_ratio", "long_name": "water vapor mixing ratio", - "coordinates": "flight_time longitude latitude pressure", "units": "g/kg", - "standard_name": "humidity_mixing_ratio", + "coordinates": "sounding_id flight_time longitude latitude pressure", "_FillValue": 9.969209968386869e+36 }, "windSpeed": { - "long_name": "wind speed", - "coordinates": "flight_time longitude latitude pressure", "standard_name": "wind_speed", + "long_name": "wind speed", "units": "m/s", + "coordinates": "sounding_id flight_time longitude latitude pressure", "_FillValue": 9.969209968386869e+36 }, "windDirection": { - "long_name": "wind direction", - "coordinates": "flight_time longitude latitude pressure", "standard_name": "wind_from_direction", - "units": "degrees", + "long_name": "wind direction", + "units": "degree", + "coordinates": "sounding_id flight_time longitude latitude pressure", "_FillValue": 9.969209968386869e+36 }, "wind_u": { - "long_name": "u-component of the wind", "standard_name": "eastward_wind", + "long_name": "u-component of the wind", "units": "m/s", "coordinates": "flight_time longitude latitude altitude", "_FillValue": 9.969209968386869e+36 }, "wind_v": { - "long_name": "v-component of the wind", "standard_name": "northward_wind", + "long_name": "v-component of the wind", "units": "m/s", - "coordinates": "flight_time longitude latitude altitude", + "coordinates": "sounding_id flight_time longitude latitude altitude", "_FillValue": 9.969209968386869e+36 }, "latitude": { @@ -109,34 +108,95 @@ "_FillValue": 9.969209968386869e+36 }, "longitude": { - "long_name": "longitude", "standard_name": "longitude", + "long_name": "longitude", "units": "degrees_east", "axis": "X", "_FillValue": 9.969209968386869e+36 }, "ascent_flag": { "long_name": "indicator of vertical flight direction", - "flag_values": "1, 0", + "flag_values": "True False", "flag_meanings": "ascending descending", - "valid_range": "0, 1" + "valid_range": "False True" }, "platform": { "long_name": "platform identifier", - "units": "-", - "description": "1: BCO, 2: Meteor, 3: RH-Brown, 4: MS-Merian, 5: Atalante" + "units": "1", + "description": "1: BCO, 2: Meteor, 3: RonBrown, 4: MS-Merian, 5: Atalante" }, "sounding_id": { "cf_role": "trajectory_id", "long_name": "sounding identifier", - "description": "unique string describing the soundings origin" + "description": "unique string describing the soundings origin (PLATFORM_SND-DIRECTION_LAT_LON_TIME)" }, "altitude": { - "standard_name": "geopotenial_height", + "standard_name": "geopotential_height", "long_name": "geopotential height retrieved from PTU", - "units": "gpm", - "coordinates": "flight_time longitude latitude pressure", + "units": "m", + "positive": "up", + "coordinates": "sounding_id flight_time longitude latitude pressure", + "_FillValue": 9.969209968386869e+36 + }, + "altitude_WGS84": { + "standard_name": "height_above_reference_ellipsoid", + "positive": "up", + "coordinates": "sounding_id flight_time longitude latitude pressure", + "units": "m", "_FillValue": 9.969209968386869e+36 + }, + "extendedVerticalSoundingSignificance":{ + "long_name": "extended vertical soudning significance", + "description": "see BUFR code flag table to decode", + "units": "" + } + }, + "global_meta_data":{ + "EUREC4A_.*?_.*?_L1":{ + "title": "EUREC4A level 1 sounding data" + }, + "EUREC4A_.*?_.*?_L2":{ + "title": "EUREC4A level 2 sounding data" + }, + "EUREC4A_.*?_.*?_.*?":{ + "campaign_id": "EUREC4A", + "doi": "10.25326/62", + "institution": "Max Planck Institute for Meteorology, Hamburg, Germany", + "references": "Stephan et al. (2020): Ship- and island-based atmospheric soundings from the 2020 EUREC4A field campaign, ESSD", + "acknowledgement": "The MPI-M is listed as the institute of first contact. Investigators, Institutions and Agencies who contributed to this dataset are numerous and are acknowledged in the reference publication" + }, + "EUREC4A_Atalante_.*?_.*?":{ + "platform_id": "Atalante", + "platform_location": "N/A" + }, + "EUREC4A_Atalante_Vaisala-RS_.*?":{ + "platform_id": "Atalante", + "instrument_id": "Vaisala-RS" + }, + "EUREC4A_Atalante_Meteomodem-RS_.*?":{ + "platform_id": "Atalante", + "instrument_id": "Meteomodem-RS", + "surface_altitude": "13 m" + }, + "EUREC4A_Meteor_.*?_.*?":{ + "platform_id": "Meteor", + "platform_location": "N/A", + "instrument_id": "Vaisala-RS" + }, + "EUREC4A_BCO_.*?_.*?":{ + "platform_id": "BCO", + "platform_location": "Deebles Point, Barbados, West Indies", + "instrument_id": "Vaisala-RS" + }, + "EUREC4A_RonBrown_.*?_.*?":{ + "platform_id": "RonBrown", + "platform_location": "N/A", + "instrument_id": "Vaisala-RS" + }, + "EUREC4A_MS-Merian_.*?_.*?":{ + "platform_id": "MS-Merian", + "platform_location": "N/A", + "instrument_id": "Vaisala-RS" } } } \ No newline at end of file diff --git a/eurec4a_snd/interpolate/batch_interpolate_soundings.py b/eurec4a_snd/interpolate/batch_interpolate_soundings.py index ff67390..aec0ff2 100644 --- a/eurec4a_snd/interpolate/batch_interpolate_soundings.py +++ b/eurec4a_snd/interpolate/batch_interpolate_soundings.py @@ -16,10 +16,14 @@ from metpy.units import units import netCDF4 from netCDF4 import num2date, default_fillvals +import pyproj import eurec4a_snd -sys.path.append(os.path.dirname(os.path.abspath(__file__))) +pwd = os.path.dirname(os.path.abspath(__file__)) +sys.path.append(pwd) from postprocessing import * +sys.path.append(pwd+'/../') +from _helpers import get_global_attrs, setup_logging def get_args(): @@ -75,146 +79,202 @@ def get_args(): return parsed_args -def setup_logging(verbose): - assert verbose in ["DEBUG", "INFO", "WARNING", "ERROR"] - logging.basicConfig( - level=logging.getLevelName(verbose), - format="%(levelname)s - %(name)s - %(funcName)s - %(message)s", - handlers=[ - logging.FileHandler("{}.log".format(__file__)), - logging.StreamHandler() - ]) - - variables_dict = {'launch_time': 'launch_time', 'flight_time': 'flight_time', 'pressure':'pressure', 'latitude': 'latitude', 'longitude': 'longitude', 'ascentRate':'ascent_rate', 'temperature': 'temperature', 'dewPoint': 'dew_point', 'windSpeed': 'wind_speed', - 'wind_u': 'wind_u', 'wind_v': 'wind_v'} + 'wind_u': 'wind_u', 'wind_v': 'wind_v' + } output_variables = ['altitude', 'temperature', 'pressure', 'dew_point', 'wind_u', 'wind_v', 'wind_speed', 'longitude', 'latitude', 'mixing_ratio', 'launch_time', 'flight_time', 'ascent_rate'] meta_data_dict = {'flight_time': {'long_name': 'time at pressure level', + 'standard_name': 'time', 'units': 'seconds since 1970-01-01 00:00:00 UTC', 'coordinates': 'flight_time longitude latitude altitude', - '_FillValue': default_fillvals['f8']}, + '_FillValue': default_fillvals['f8'], + 'ancillary_variables': 'N_ptu m_ptu', + 'cell_methods': 'altitude: mean (interval: 10 m comment: m_ptu)' + }, 'launch_time': {'long_name': 'time at which the sounding started', 'units': 'seconds since 1970-01-01 00:00:00 UTC', '_FillValue': default_fillvals['f8'] }, 'ascent_rate': {'long_name': 'ascent/desent rate of measuring device', - 'coordinates': 'flight_time longitude latitude altitude', + 'coordinates': 'longitude latitude flight_time sounding_id', 'description': 'calculated from interpolated geopotential height changes', - 'units': 'gpm/s', - '_FillValue': default_fillvals['f4']}, + 'units': 'm/s', + '_FillValue': default_fillvals['f4'], + 'ancillary_variables': 'N_ptu m_ptu', + 'cell_methods': 'altitude: mean (interval: 10 m comment: m_ptu)' + }, 'altitude': {'long_name': 'geopotential height', 'standard_name': 'geopotential_height', - 'units': 'gpm', - 'axis': 'Y', - '_FillValue': default_fillvals['f4'] + 'units': 'm', + 'axis': 'Z', + 'positive': 'up', + 'bounds': 'alt_bnds', + '_FillValue': False }, + 'alt_bnds': { + '_FillValue': False, + 'comment': '(lower bound, upper bound]', + }, 'pressure': {'long_name': 'pressure', 'standard_name': 'air_pressure', 'units': 'hPa', - 'coordinates': 'flight_time longitude latitude altitude', - '_FillValue': default_fillvals['f4'] + 'coordinates': 'longitude latitude flight_time sounding_id', + '_FillValue': default_fillvals['f4'], + 'ancillary_variables': 'N_ptu m_ptu', + 'cell_methods': 'altitude: mean (interval: 10 m comment: m_ptu)' }, 'temperature': {'long_name': 'dry bulb temperature', 'standard_name': 'air_temperature', - 'units': 'degrees Celsius', - 'coordinates': 'flight_time longitude latitude altitude', - '_FillValue': default_fillvals['f4'] + 'units': 'degrees_Celsius', + 'coordinates': 'longitude latitude flight_time sounding_id', + '_FillValue': default_fillvals['f4'], + 'cell_methods': 'altitude: point (derived from averaged theta)' }, 'theta': {'long_name': 'potential temperature', 'standard_name': 'air_potential_temperature', 'units': 'K', - 'coordinates': 'flight_time longitude latitude altitude', - '_FillValue': default_fillvals['f4'] + 'coordinates': 'flight_time longitude latitude sounding_id', + '_FillValue': default_fillvals['f4'], + 'ancillary_variables': 'N_ptu m_ptu', + 'cell_methods': 'altitude: mean (interval: 10 m comment: m_ptu)' }, 'relative_humidity': {'long_name': 'relative_humidity', 'standard_name': 'relative_humidity', - 'coordinates': 'flight_time longitude latitude altitude', + 'coordinates': 'flight_time longitude latitude sounding_id', 'units': '%', - '_FillValue': default_fillvals['f4'] + '_FillValue': default_fillvals['f4'], + 'cell_methods': 'altitude: point (derived from averaged q, T, p)' }, 'specific_humidity': {'long_name': 'specific humidity', 'standard_name': 'specific_humidity', 'units': 'g/kg', - 'coordinates': 'flight_time longitude latitude altitude', - '_FillValue': default_fillvals['f4'] + 'coordinates': 'flight_time longitude latitude sounding_id', + '_FillValue': default_fillvals['f4'], + 'ancillary_variables': 'N_ptu m_ptu', + 'cell_methods': 'altitude: mean (interval: 10 m comment: m_ptu)' }, 'dew_point': {'long_name': 'dew point temperature', 'standard_name': 'dew_point_temperature', - 'units': 'degrees Celsius', - 'coordinates': 'flight_time longitude latitude altitude', - '_FillValue': default_fillvals['f4']}, + 'units': 'degrees_Celsius', + 'coordinates': 'flight_time longitude latitude sounding_id', + '_FillValue': default_fillvals['f4'], + 'ancillary_variables': 'N_ptu m_ptu', + 'cell_methods': 'altitude: mean (interval: 10 m comment: m_ptu)' + }, 'mixing_ratio': {'long_name': 'water vapor mixing ratio', - 'coordinates': 'flight_time longitude latitude altitude', + 'coordinates': 'flight_time longitude latitude sounding_id', 'units': 'g/kg', 'standard_name': 'humidity_mixing_ratio', - '_FillValue': default_fillvals['f4'] + '_FillValue': default_fillvals['f4'], + 'ancillary_variables': 'N_ptu', + 'cell_methods': 'altitude: point (derived from averaged q)' }, 'wind_speed': {'long_name': 'wind speed', 'standard_name': 'wind_speed', 'units': 'm/s', - 'coordinates': 'flight_time longitude latitude altitude', - '_FillValue': default_fillvals['f4'] + 'coordinates': 'flight_time longitude latitude sounding_id', + '_FillValue': default_fillvals['f4'], + 'cell_methods': 'altitude: point (derived from averaged u, v)' }, 'wind_direction': {'long_name': 'wind direction', - 'coordinates': 'flight_time longitude latitude altitude', + 'coordinates': 'flight_time longitude latitude sounding_id', 'standard_name': 'wind_from_direction', - 'units': 'degrees', - '_FillValue': default_fillvals['f4'] + 'units': 'degree', + '_FillValue': default_fillvals['f4'], + 'cell_methods': 'altitude: point (derived from averaged u, v)' }, 'wind_u': {'long_name': 'u-component of the wind', 'standard_name': 'eastward_wind', 'units': 'm/s', - 'coordinates': 'flight_time longitude latitude altitude', - '_FillValue': default_fillvals['f4'] + 'coordinates': 'flight_time longitude latitude sounding_id', + '_FillValue': default_fillvals['f4'], + 'ancillary_variables': 'N_gps m_gps', + 'cell_methods': 'altitude: mean (interval: 10 m comment: m_gps)' }, 'wind_v': {'long_name': 'v-component of the wind', 'standard_name': 'northward_wind', 'units': 'm/s', - 'coordinates': 'flight_time longitude latitude altitude', - '_FillValue': default_fillvals['f4'] + 'coordinates': 'flight_time longitude latitude sounding_id', + '_FillValue': default_fillvals['f4'], + 'ancillary_variables': 'N_gps m_gps', + 'cell_methods': 'altitude: mean (interval: 10 m comment: m_gps)' }, 'latitude': {'long_name': 'latitude', 'standard_name': 'latitude', - 'units': 'degree_north', - '_FillValue': default_fillvals['f4'] + 'units': 'degrees_north', + '_FillValue': default_fillvals['f4'], + 'ancillary_variables': 'N_gps m_gps', + 'cell_methods': 'altitude: mean (interval: 10 m comment: m_gps)' }, 'longitude': {'long_name': 'longitude', 'standard_name': 'longitude', - 'units': 'degree_east', - '_FillValue': default_fillvals['f4'] + 'units': 'degrees_east', + '_FillValue': default_fillvals['f4'], + 'ancillary_variables': 'N_gps m_gps', + 'cell_methods': 'altitude: mean (interval: 10 m comment: m_gps)', }, 'ascent_flag': {'long_name': 'indicator of vertical flight direction', - 'flag_values': '1, 0', + 'flag_values': np.array([1, 0],dtype=np.int16), 'flag_meanings': 'ascending descending', - 'valid_range': '0, 1' + 'valid_range': np.array([0, 1],dtype=np.int16) }, - 'platform': {'long_name': 'platform identifier', - 'units': '-', - 'description': '1: BCO, 2: Meteor, 3: RH-Brown, 4: MS-Merian, 5: Atalante' + 'platform_id': {'long_name': 'platform identifier', + 'units': '1', + 'flag_values': np.array([1, 2, 3, 4, 5], dtype='int16'), + 'flag_meanings': 'BCO Meteor RonBrown MS-Merian Atalante', + 'description': '1: BCO, 2: Meteor, 3: RonBrown, 4: MS-Merian, 5: Atalante' }, - 'data_count': {'description': 'number of measurements that have been used to derive level 2 data point average', - 'flag_values':'np.nan, 0 , 1+', - 'flag_meanings':'no data, linear interpolation, averaging'} + 'N_ptu': {'standard_name': 'number_of_observations', + 'description': 'number of observations used to derive level 2 PTU-data average', + 'units': '1', + 'coordinates': 'flight_time longitude latitude sounding_id', + '_FillValue': default_fillvals['f4'] + }, + 'N_gps': {'standard_name': 'number_of_observations', + 'description': 'number of observations used to derive level 2 GPS-data average', + 'units': '1', + 'coordinates': 'flight_time longitude latitude sounding_id', + '_FillValue': default_fillvals['f4'] + }, + 'm_ptu': {'long_name': 'bin method', + 'description': 'method used to derive level 2 PTU-data average', + 'flag_values': np.array([0, 1 , 2], dtype='int16'), + 'flag_meanings':'no_data interpolation averaging' + }, + 'm_gps': {'long_name': 'bin_method', + 'description': 'method used to derive level 2 GPS-data average', + 'flag_values': np.array([0, 1 , 2], dtype='int16'), + 'flag_meanings':'no_data interpolation averaging' + }, + 'x':{'long_name':'WGS84_x'}, + 'y':{'long_name':'WGS84_y'}, + 'z':{'long_name':'WGS84_z'} } platform_rename_dict = {'Atalante (ATL)': 'Atalante', 'Barbados Cloud Observatory (BCO)': 'BCO', 'Maria S Merian (MER)': 'MS-Merian', 'FS_METEOR (MET)': 'Meteor', - 'Research Vessel Ronald H. Brown (WTEC) (RHB)': 'RH-Brown'} -platform_number_dict = {'Atalante (ATL)': 5, - 'Barbados Cloud Observatory (BCO)': 1, - 'Maria S Merian (MER)': 4, - 'FS_METEOR (MET)': 2, - 'Research Vessel Ronald H. Brown (WTEC) (RHB)': 3} + 'Research Vessel Ronald H. Brown (WTEC) (RHB)': 'RonBrown'} +platform_number_dict = {'Atalante': 5, + 'BCO': 1, + 'MS-Merian': 4, + 'Meteor': 2, + 'RonBrown': 3} std_pressures = [1000.00, 925.00, 850.00, 700.00, 500.00, 400.00, 300.00, 250.00, 200.00, 150.00, 100.00, 70.00, 50.00] +interpolation_grid = np.arange(0, 31000, 10) # meters +interpolation_bins = np.arange(-5,31005,10).astype('int') # Bins len(interpolation_grid)+1; (a,b]; (meters) +max_gap_fill = 50 # Maximum data gap size that should be filled by interpolation (meters) + +json_config_fn = pwd+'/../config/mwx_config.json' + def main(args={}): """ @@ -291,15 +351,15 @@ def main(args={}): ds = ds.isel({'levels': uniq_altitude_idx}) # Check if platform is known - if ds.platform_name not in platform_rename_dict.keys(): - logging.error('The platform {} is not known. Please choose one of {}'.format(ds.platform_name, platform_rename_dict.keys())) - sys.exit() + # if ds.platform_id not in platform_rename_dict.keys(): + # logging.error('The platform {} is not known. Please choose one of {}'.format(ds.platform_id, platform_rename_dict.keys())) + # sys.exit() # Consistent platform test if f == 0: - platform = ds.platform_name + platform = ds.platform_id else: - assert ds.platform_name == platform, 'The platform seems to change from {} to {}'.format(platform, ds.platform_name) + assert ds.platform_id == platform, 'The platform seems to change from {} to {}'.format(platform, ds.platform_id) # Unique levels test if len(ds.altitude) != len(np.unique(ds.altitude)): @@ -307,24 +367,31 @@ def main(args={}): break # Prepare some data that cannot be linearly interpolated - wind_dir_rad = np.deg2rad(ds.windDirection.values) - ds['wind_u'] = -1*ds.windSpeed.values * np.sin(wind_dir_rad) - ds['wind_v'] = -1*ds.windSpeed.values * np.cos(wind_dir_rad) - - dewPoint_K = ds['dewPoint'].values+273.15 - pressure_Pa = ds['pressure'].values*100 + u, v = get_wind_components(ds.windDirection.values, ds.windSpeed.values) + ds['wind_u'] = xr.DataArray(u, dims=['levels']) + ds['wind_v'] = xr.DataArray(v, dims=['levels']) + + if 'altitude_WGS84' in ds.keys(): + # Convert lat, lon, alt to cartesian coordinates + ecef = pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84') + lla = pyproj.Proj(proj='latlong', ellps='WGS84', datum='WGS84') + x, y, z = pyproj.transform(lla, ecef, + ds.longitude.values, + ds.latitude.values, + ds.altitude_WGS84.values, + radians=False) + variables_dict.update({'x':'x', 'y':'y', 'z':'z'}) + for var, val in {'x':x, 'y':y, 'z':z}.items(): + ds[var] = xr.DataArray(val, dims=['levels']) + else: + logging.warning('No WGS84 altitude could be found. The averaging of the position might be faulty especially at the 0 meridian and close to the poles') theta = calc_theta_from_T(ds['temperature'].values, ds['pressure'].values) - #dp = mpcalc.dewpoint_from_relative_humidity( - #ds.temperature.values * units.degC, ds.humidity.values / 100).magnitude - q = calc_q_from_rh(ds['dewPoint'].values, ds['pressure'].values) - e_s = calc_saturation_pressure(ds['temperature'].values+273.15) - w_s = mpcalc.mixing_ratio(e_s*units.Pa, ds['pressure'].values*units.hPa).magnitude + e_s = calc_saturation_pressure(ds['temperature'].values + 273.15) + w_s = mpcalc.mixing_ratio(e_s * units.Pa, ds['pressure'].values*units.hPa).magnitude w = ds['humidity'].values/100.*w_s q = w/(1+w) - #mixing_ratio = np.array(calc_mixing_ratio_hardy(dewPoint_K, - # pressure_Pa))*1000 - #q = mpcalc.specific_humidity_from_mixing_ratio(mixing_ratio/1000) + da_w = xr.DataArray([w*1000], dims=['sounding', 'altitude'], coords={'altitude': ds.altitude.values}) @@ -359,24 +426,50 @@ def main(args={}): ds_new = ds_new.dropna(dim='altitude', subset=output_variables, how='any') - ds_interp = ds_new.interp(altitude=np.arange(0, 31000, 10)) + ds_interp = ds_new.interp(altitude=interpolation_grid) elif args['method'] == 'bin': - ds_interp = ds_new.groupby_bins('altitude',np.arange(-5,31005,10), labels=np.arange(0,31000,10), restore_coord_dims=True).mean() + ds_interp = ds_new.groupby_bins('altitude', interpolation_bins, + labels=interpolation_grid, + restore_coord_dims=True).mean() + ds_interp = ds_interp.transpose() ds_interp = ds_interp.rename({'altitude_bins':'altitude'}) + # Create bounds variable + ds_interp['alt_bnds'] = xr.DataArray(np.array([interpolation_bins[:-1],interpolation_bins[1:]]).T, + dims=['altitude','nv'], + coords={'altitude': ds_interp.altitude.values} + ) + ds_interp['launch_time'] = ds_new['launch_time'] ## Interpolation NaN - ds_interp = ds_interp.interpolate_na('altitude', max_gap=50, use_coordinate=True) + ds_interp = ds_interp.interpolate_na('altitude', max_gap=max_gap_fill, use_coordinate=True) dims_2d = ['sounding', 'altitude'] coords_1d = {'altitude': ds_interp.altitude.values} wind_u = ds_interp.isel({'sounding': 0})['wind_u'] wind_v = ds_interp.isel({'sounding': 0})['wind_v'] - wind_direction = np.rad2deg(np.arctan2(-1*wind_u, -1*wind_v)) % 360 - ds_interp['wind_direction'] = xr.DataArray([np.array(wind_direction.values)], + dir, wsp = get_directional_wind(wind_u, wind_v) + ds_interp['wind_direction'] = xr.DataArray([np.array(dir)], + dims=dims_2d, + coords=coords_1d) + ds_interp['wind_speed'] = xr.DataArray([np.array(wsp)], dims=dims_2d, coords=coords_1d) + if 'altitude_WGS84' in ds.keys(): + lon, lat, alt = pyproj.transform(ecef, lla, + ds_interp['x'].values[0], + ds_interp['y'].values[0], + ds_interp['z'].values[0], + radians=False) + for var, val in {'latitude':lat, 'longitude':lon, 'altitude_WGS84': alt}.items(): + ds_interp[var] = xr.DataArray([val], dims=dims_2d, coords=coords_1d) + + del ds_interp['x'] + del ds_interp['y'] + del ds_interp['z'] + del ds_interp['altitude_WGS84'] + ds_input = ds_input.sortby('altitude') ds_input.altitude.load() ds_input.pressure.load() @@ -387,27 +480,12 @@ def main(args={}): dims=dims_2d, coords=coords_1d) - # Calculations after interpolation - #specific_humidity = metpy.calc.specific_humidity_from_mixing_ratio(ds_interp['mixing_ratio']/1000) - #relative_humidity = metpy.calc.relative_humidity_from_specific_humidity( - # specific_humidity, ds_interp.temperature.values * units.degC, - # ds_interp.pressure.values * units.hPa) - - #ds_interp['specific_humidity'] = xr.DataArray(np.array(specific_humidity*1000), - # dims=dims_2d, - # coords=coords_1d) - #ds_interp['relative_humidity'] = xr.DataArray(np.array(relative_humidity)*100, - # dims=dims_2d, - # coords=coords_1d) ds_interp['launch_time'] = xr.DataArray([ds_interp.isel({'sounding': 0}).launch_time.item()/1e9], dims=['sounding']) - ds_interp['platform'] = xr.DataArray([platform_number_dict[platform]], + ds_interp['platform_id'] = xr.DataArray([platform_number_dict[platform]], dims=['sounding']) - #ds_interp['ascent_rate'] = xr.DataArray([calc_ascentrate(ds_interp.isel({'sounding': 0}).altitude.values, - # ds_interp.isel({'sounding': 0}).flight_time.values)], - # dims=dims_2d, - # coords=coords_1d) + # Calculations after interpolation # Recalculate temperature and relative humidity from theta and q temperature = calc_T_from_theta(ds_interp.isel(sounding=0)['theta'].values, ds_interp.isel(sounding=0)['pressure'].values) ds_interp['temperature'] = xr.DataArray([np.array(temperature)], @@ -418,28 +496,39 @@ def main(args={}): e_s = calc_saturation_pressure(ds_interp.isel(sounding=0)['temperature'].values+273.15) w_s = mpcalc.mixing_ratio(e_s*units.Pa, ds_interp.isel(sounding=0)['pressure'].values*units.hPa).magnitude relative_humidity = w/w_s*100 - #relative_humidity = calc_rh_from_q(ds_interp['q'].values, ds_interp['temperature_re'].values, ds_interp['pressure'].values) - #saturation_p = calc_saturation_pressure(ds_interp['temperature_re'].values[0,:]+273.15) - #w_s = mpcalc.mixing_ratio(saturation_p * units.Pa, ds_interp['pressure'].values * units.hPa).magnitude - #relative_humidity = (ds_interp['specific_humidity']/1000)/((1-ds_interp['specific_humidity']/1000)*w_s)*100 ds_interp['relative_humidity'] = xr.DataArray([np.array(relative_humidity)], dims=dims_2d, coords=coords_1d) - # Interpolate NaNs - ## max_gap is the maximum gap of NaNs in meters that will be still interpolated - #ds_interp = ds_interp.interpolate_na('altitude', max_gap=50, use_coordinate=True) - ds_interp['data_count'] = xr.DataArray(ds_new.pressure.groupby_bins('altitude',np.arange(-5,31005,10), labels=np.arange(0,31000,10), restore_coord_dims=True).count().values, - dims=dims_2d, - coords=coords_1d) - data_avail_or_interp = np.where(~np.isnan(ds_interp.pressure), 0, np.nan) - stacked_data_counts = np.vstack([ds_interp.data_count.values[0,:], data_avail_or_interp[0,:]]) - nan_idx_both = np.logical_and(np.isnan(stacked_data_counts[0]), np.isnan(stacked_data_counts[1])) - data_counts_combined = np.empty(len(stacked_data_counts[0])) - data_counts_combined.fill(np.nan) - data_counts_combined[~nan_idx_both] = np.nanmax(stacked_data_counts[:,~nan_idx_both], axis=0) - ds_interp['data_count'].values[0,:] = data_counts_combined + # Count number of measurements within each bin + ds_interp['N_ptu'] = xr.DataArray( + ds_new.pressure.groupby_bins('altitude', interpolation_bins, labels=interpolation_grid, + restore_coord_dims=True).count().values, + dims=dims_2d, + coords=coords_1d) + ds_interp['N_gps'] = xr.DataArray( + ds_new.latitude.groupby_bins('altitude', interpolation_bins, labels=interpolation_grid, + restore_coord_dims=True).count().values, + dims=dims_2d, + coords=coords_1d) + + # Cell method used + data_exists = np.where(np.isnan(ds_interp.pressure), False, True) + data_mean = np.where(np.isnan(ds_interp.N_ptu), False, True) # no data or interp: nan; mean > 0 + data_method = np.zeros_like(data_exists, dtype='uint') + data_method[np.logical_and(data_exists, data_mean)] = 2 + data_method[np.logical_and(data_exists, ~data_mean)] = 1 + ds_interp['m_ptu'] = xr.DataArray(data_method, dims=dims_2d, coords=coords_1d) + ds_interp['N_ptu'].values[np.logical_and(~data_mean, data_method > 0)] = 0 + + data_exists = np.where(np.isnan(ds_interp.latitude), False, True) + data_mean = np.where(np.isnan(ds_interp.N_gps), False, True) # no data or interp: nan; mean > 0 + data_method = np.zeros_like(data_exists, dtype='uint') + data_method[np.logical_and(data_exists, data_mean)] = 2 + data_method[np.logical_and(data_exists, ~data_mean)] = 1 + ds_interp['m_gps'] = xr.DataArray(data_method, dims=dims_2d, coords=coords_1d) + ds_interp['N_gps'].values[np.logical_and(~data_mean, data_method > 0)] = 0 direction = get_direction(ds_interp, ds) if direction == 'ascending': @@ -447,17 +536,21 @@ def main(args={}): else: ds_interp['ascent_flag'] = xr.DataArray([0], dims=['sounding']) + # Copy trajectory id from level1 dataset + ds_interp['sounding_id'] = xr.DataArray([ds['sounding_id'].values], dims=['sounding']) + ds_interp.sounding_id.attrs = ds['sounding_id'].attrs + script_basename = os.path.basename(__file__) script_modification_time = time.ctime(os.path.getmtime(os.path.realpath(__file__))) - global_attrs_dict = {'title': 'EUREC4A interpolated sounding data', - 'platform_name': platform, + glob_attrs_dict = {'title': 'EUREC4A interpolated sounding data', + 'platform_id': platform, 'surface_altitude': ds.attrs['surface_altitude'], 'instrument': ds.instrument, 'doi': 'pending', 'created_with': '{file} with its last modifications on {time}'. format(time=script_modification_time, file=script_basename), - 'git-version': git_module_version, + 'git_version': git_module_version, 'python_version': "{} (with numpy:{}, netCDF4:{}, eurec4a_snd:{})". format(sys.version, np.__version__, netCDF4.__version__, __version__), 'created_on': str(time.ctime(time.time())), @@ -465,31 +558,52 @@ def main(args={}): 'Conventions': 'CF-1.7' } - ds_interp = set_global_attributes(ds_interp, global_attrs_dict) + # Overwrite standard attrs with those defined in config file + # Get global meta data from mwx_config.json + level='L2' + glob_attrs_dict2 = get_global_attrs(json_config_fn, f'{ds.campaign_id}_{ds.platform_id}_{ds.instrument_id}_{level}') + + for attrs, value in glob_attrs_dict2.items(): + glob_attrs_dict[attrs] = value + + ds_interp = set_global_attributes(ds_interp, glob_attrs_dict) ds_interp = set_additional_var_attributes(ds_interp, meta_data_dict) for variable in ['temperature', 'dew_point', 'wind_speed', 'pressure', 'wind_u', 'wind_v', 'latitude', 'longitude', - 'mixing_ratio', 'altitude', 'wind_direction', + 'mixing_ratio', 'wind_direction', 'specific_humidity', 'relative_humidity', 'ascent_rate' ]: # ds_interp[variable].values = np.round(ds_interp[variable].values, 2) - ds_interp[variable].encoding = {'dtype': 'f4'} - ds_interp['ascent_flag'].encoding = {'dtype': 'bool'} - ds_interp['platform'].encoding = {'dtype': 'uint8'} + ds_interp[variable].encoding['dtype'] = 'f4' + ds_interp['ascent_flag'].encoding['dtype'] = 'int16' + ds_interp['platform_id'].encoding['dtype'] = 'int16' + ds_interp['m_gps'].encoding['dtype'] = 'int16' + ds_interp['m_ptu'].encoding['dtype'] = 'int16' + ds_interp['N_gps'].encoding['dtype'] = 'f4' + ds_interp['N_ptu'].encoding['dtype'] = 'f4' + ds_interp['alt_bnds'].encoding['dtype'] = 'int16' + ds_interp['altitude'].encoding['dtype'] = 'int16' # Transpose dataset if necessary for variable in ds_interp.data_vars: + if variable == 'alt_bnds': continue dims = ds_interp[variable].dims if (len(dims) == 2) and (dims[0] != 'sounding'): ds_interp[variable] = ds_interp[variable].T time_dt = num2date(ds_interp.isel({'sounding': 0}).launch_time, "seconds since 1970-01-01 00:00:00") - time_fmt = time_dt.strftime('%Y%m%d%H%M') - platform_filename = platform_rename_dict[platform] - outfile = args['outputfolder']+'EUREC4A_{platform}_soundings_{date}.nc'.format(platform=platform_filename, date=time_fmt) + time_fmt = time_dt.strftime('%Y%m%dT%H%M') + platform_filename = platform # platform_rename_dict[platform] + outfile = args['outputfolder']+'{campaign}_{platform}_{instrument_id}_{level}_{date}_{version}.nc'.format(campaign=ds.campaign_id, + platform=platform_filename, + instrument_id=ds.instrument_id, + level=level, + date=time_fmt, + version=git_module_version + ) logging.info('Write output to {}'.format(outfile)) write_dataset(ds_interp, outfile) diff --git a/eurec4a_snd/interpolate/postprocessing.py b/eurec4a_snd/interpolate/postprocessing.py index a1b78ea..2a7dfea 100644 --- a/eurec4a_snd/interpolate/postprocessing.py +++ b/eurec4a_snd/interpolate/postprocessing.py @@ -147,6 +147,59 @@ def calc_rh_from_q(q, T, p): return rh +def get_wind_components(dir, spd): + """ + Convert directional windspeed and winddirection + to wind components + + Input + ----- + dir : array-like + wind from direction (deg) + spd : array-like + wind speed along dir + + Return + ------ + u : array-like + eastward wind component + (positive when directed eastward) + v : array-like + northward wind component + (positive when directed northward) + """ + dir_rad = np.deg2rad(dir) + u = -1*spd * np.sin(dir_rad) + v = -1*spd * np.cos(dir_rad) + return u, v + + +def get_directional_wind(u, v): + """ + Convert wind-components to directional + wind + + Input + ----- + u : array-like + eastward wind component + (positive when directed eastward) + v : array-like + northward wind component + (positive when directed northward) + + Return + ------ + dir : array-like + wind from direction (deg) + spd : array-like + wind speed along dir + """ + dir = np.rad2deg(np.arctan2(-1*u, -1*v)) % 360 + spd = np.sqrt(u**2+v**2) + return dir, spd + + def set_global_attributes(ds, global_attrs_dict): """ diff --git a/postprocessing/batch_reprocess.sh b/postprocessing/batch_reprocess.sh index e572ce8..03ad59b 100755 --- a/postprocessing/batch_reprocess.sh +++ b/postprocessing/batch_reprocess.sh @@ -11,10 +11,13 @@ module unload netcdf_c module load anaconda3 -source activate new_campaign +source activate /home/mpim/m300408/conda-envs/new_campaign cd ~/GITHUB/eurec4a_snd/eurec4a_snd/ -git checkout master +git checkout final_changes + +rm -r /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/ +rm -r /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/ # Correct mwx files from the Meteor ## Get surface values from the DSHIP data @@ -42,87 +45,100 @@ git checkout master # Convert Level0 to Level1 -git checkout postprocessing +#git checkout postprocessing ## Original METEOR data #python L1_mwx.py -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v3/level0/MET/MWX/DWD/ -o /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level1_mwx/ --platform_name_long FS_METEOR --platform_name_short MET --platform_location NA #python L1_mwx.py -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v3/level0/MET/ -o /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level1_mwx/ --platform_name_long FS_METEOR --platform_name_short MET --platform_location NA ## Corrected METEOR data -#python L1_mwx.py -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level0_corrected/MET/ -o /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level1_mwx/ --platform_name_long FS_METEOR --platform_name_short MET --platform_location NA -#python L1_mwx.py -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v3/level0/MER/ -o /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level1_mwx/ --platform_name_long "Maria S Merian" --platform_name_short MER --platform_location NA -#python L1_mwx.py -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v3/level0/ATL/Vaisala/ -o /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level1_mwx/ --platform_name_long "Atalante" --platform_name_short ATL --platform_location NA -#python L1_mwx.py -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v3/level0/RHB/ -o /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level1_mwx/ --platform_name_long "Research Vessel Ronald H. Brown (WTEC)" --platform_name_short RHB --platform_location NA -#python L1_mwx.py -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v3/level0/BCO/ -o /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level1_mwx/ --platform_name_long "Barbados Cloud Observatory" --platform_name_short BCO --platform_location "Deebles Point, Barbados, West Indies" +## Version_5 +python L1_mwx.py -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level0_corrected/MET/ -o /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level1_mwx/ --platform_id Meteor --campaign EUREC4A --instrument_id Vaisala-RS +python L1_mwx.py -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v3/level0/MER/ -o /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level1_mwx/ --platform_id MS-Merian --campaign EUREC4A --instrument_id Vaisala-RS +python L1_mwx.py -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v3/level0/ATL/Vaisala/ -o /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level1_mwx/ --platform_id Atalante --campaign EUREC4A --instrument_id Vaisala-RS +python L1_mwx.py -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v3/level0/RHB/ -o /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level1_mwx/ --platform_id RonBrown --campaign EUREC4A --instrument_id Vaisala-RS +python L1_mwx.py -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v3/level0/BCO/ -o /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level1_mwx/ --platform_id BCO --campaign EUREC4A --instrument_id Vaisala-RS +python L1_bufr.py -p "/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v3/level0/ATL/MeteoModem/BFR/" -o /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level1_bufr/ --platform_id Atalante --campaign EUREC4A --instrument_id Meteomodem-RS -s 1 5 6 7 14 + # Remove short soundings (less than 30 levels) -python remove_short_soundings.py -i '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level1_mwx/EUREC4A*sounding*.nc' -t 30 -d True +python ../postprocessing/remove_short_soundings.py -i '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level1_mwx/EUREC4A*L1*.nc' -t 30 -d True +python ../postprocessing/remove_short_soundings.py -i '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level1_bufr/EUREC4A*L1*.nc' -t 30 -d True # Convert Level1 to Level2 (interpolate) -git checkout interpolation_method -python ./interpolate/batch_interpolate_soundings.py -m bin -i '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level1_mwx/EUREC4A_MET_sounding*.nc' -o '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/' -python ./interpolate/batch_interpolate_soundings.py -m bin -i '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level1_mwx/EUREC4A_RHB_sounding*.nc' -o '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/' -python ./interpolate/batch_interpolate_soundings.py -m bin -i '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level1_mwx/EUREC4A_MER_sounding*.nc' -o '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/' -python ./interpolate/batch_interpolate_soundings.py -m bin -i '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level1_mwx/EUREC4A_BCO_sounding*.nc' -o '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/' -python ./interpolate/batch_interpolate_soundings.py -m bin -i '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level1_mwx/EUREC4A_ATL_sounding*.nc' -o '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/' +mkdir -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level1_mwx/ +python ./interpolate/batch_interpolate_soundings.py -m bin -i '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level1_mwx/EUREC4A_Meteor_*Vaisala-RS*.nc' -o '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/' +python ./interpolate/batch_interpolate_soundings.py -m bin -i '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level1_mwx/EUREC4A_RonBrown_*Vaisala-RS*.nc' -o '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/' +python ./interpolate/batch_interpolate_soundings.py -m bin -i '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level1_mwx/EUREC4A_MS-Merian_*Vaisala-RS*.nc' -o '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/' +python ./interpolate/batch_interpolate_soundings.py -m bin -i '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level1_mwx/EUREC4A_BCO_*Vaisala-RS*.nc' -o '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/' +python ./interpolate/batch_interpolate_soundings.py -m bin -i '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level1_mwx/EUREC4A_Atalante_*Vaisala-RS*.nc' -o '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/' +python ./interpolate/batch_interpolate_soundings.py -m bin -i '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level1_bufr/EUREC4A_Atalante_Meteomodem-RS_L1*_20200*.nc' -o '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_bufr/' # Concatenate Level2 files by platform module load nco -rm /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_Meteor_soundings.nc -ncrcat -h /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/*Meteor_soundings_20*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_Meteor_soundings.nc -rm /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_BCO_soundings.nc -ncrcat -h /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/*BCO_soundings_20*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_BCO_soundings.nc -rm /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_Atalante_soundings_Vaisala.nc -ncrcat -h /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/*Atalante_soundings_20*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_Atalante_soundings_Vaisala.nc -rm /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_MS-Merian_soundings.nc -ncrcat -h /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/*Merian_soundings_20*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_MS-Merian_soundings.nc -rm /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_RH-Brown_soundings.nc -ncrcat -h /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/*RH-Brown_soundings_20*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_RH-Brown_soundings.nc - -# Checks -python ./tests/sounding_graphical_comparison.py - -# Export to different directory -mkdir /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/ -mkdir -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_2/ -cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_RH-Brown_soundings.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_2/EUREC4A_RH-Brown_soundings.nc -cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_Atalante_soundings_Vaisala.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_2/EUREC4A_Atalante_soundings_Vaisala.nc -cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_Meteor_soundings.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_2/EUREC4A_Meteor_soundings.nc -cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_BCO_soundings.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_2/EUREC4A_BCO_soundings.nc -cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_MS-Merian_soundings.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_2/EUREC4A_MS-Merian_soundings.nc -cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v3/level2/EUREC4A_Atalante_soundings_MeteoModem.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_2/EUREC4A_Atalante_soundings_MeteoModem.nc - -mkdir -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_1/ - -mkdir -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_1/MET/Vaisala/ -cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level1_mwx/EUREC4A_MET_sounding*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_1/MET/Vaisala/ -mkdir -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_1/MER/Vaisala/ -cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level1_mwx/EUREC4A_MER_sounding*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_1/MER/Vaisala/ -mkdir -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_1/ATL/Vaisala/ -cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level1_mwx/EUREC4A_ATL_sounding*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_1/ATL/Vaisala/ -mkdir -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_1/RHB/Vaisala/ -cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level1_mwx/EUREC4A_RHB_sounding*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_1/RHB/Vaisala/ -mkdir -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_1/BCO/Vaisala/ -cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level1_mwx/EUREC4A_BCO_sounding*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_1/BCO/Vaisala/ -mkdir -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_1/ATL/MeteoModem/ -cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v3/level1_bufr/ATL/MeteoModem/EUREC4A_Atalante_sounding*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_1/ATL/MeteoModem/ +rm /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/EUREC4A_Meteor_soundings.nc +ncrcat -h /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/*Meteor_*Vaisala-RS*_20*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/EUREC4A_Meteor_soundings.nc +rm /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/EUREC4A_BCO_soundings.nc +ncrcat -h /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/*BCO_*Vaisala-RS*_20*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/EUREC4A_BCO_soundings.nc +rm /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/EUREC4A_Atalante_soundings_Vaisala.nc +ncrcat -h /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/*Atalante_*Vaisala-RS*_20*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/EUREC4A_Atalante_soundings_Vaisala.nc +rm /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/EUREC4A_MS-Merian_soundings.nc +ncrcat -h /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/*Merian_*Vaisala-RS*_20*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/EUREC4A_MS-Merian_soundings.nc +rm /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/EUREC4A_RonBrown_soundings.nc +ncrcat -h /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/*RonBrown_*Vaisala-RS*_20*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/EUREC4A_RonBrown_soundings.nc +rm /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_bufr/EUREC4A_Atalante_soundings_Meteomodem.nc +ncrcat -h /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_bufr/EUREC4A_Atalante_Meteomodem-RS*_20*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_bufr/EUREC4A_Atalante_soundings_Meteomodem.nc -mkdir -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_0/ -cp -r /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level0_* /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_0/ - - -python ../eurec4a_snd/interpolate/batch_interpolate_soundings.py -m bin -i '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v3/level1_bufr/ATL/MeteoModem/EUREC4A_Atalante_sounding_ascent_20200*.nc' -o '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/' -ncrcat -h /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/EUREC4A_Atalante_soundings_20200* /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_2/EUREC4A_Atalante_soundings_MeteoModem.nc - -python remove_lower_levels.py - -cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_RH-Brown_soundings.nc2 /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_2/EUREC4A_RH-Brown_soundings.nc -cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_Atalante_soundings_Vaisala.nc2 /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_2/EUREC4A_Atalante_soundings_Vaisala.nc -cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_Meteor_soundings.nc2 /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_2/EUREC4A_Meteor_soundings.nc -cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_BCO_soundings.nc2 /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_2/EUREC4A_BCO_soundings.nc -cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_MS-Merian_soundings.nc2 /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_2/EUREC4A_MS-Merian_soundings.nc - -mv /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_2/EUREC4A_Atalante_soundings_MeteoModem.nc2 /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_2/EUREC4A_Atalante_soundings_MeteoModem.nc - -python rename_variable_names.py +# Checks +#python ./tests/sounding_graphical_comparison.py + +# Copy level0 data to export folder +mkdir -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_0/ +cp -r /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level0_* /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_0/ + +# Copy level1 data to export folder +mkdir -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_1/ +mkdir -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_1/MET/Vaisala/ +cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level1_mwx/EUREC4A_Meteor_Vaisala-RS_L1*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_1/MET/Vaisala/ +mkdir -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_1/MER/Vaisala/ +cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level1_mwx/EUREC4A_MS-Merian_Vaisala-RS_L1*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_1/MER/Vaisala/ +mkdir -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_1/ATL/Vaisala/ +cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level1_mwx/EUREC4A_Atalante_Vaisala-RS_L1*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_1/ATL/Vaisala/ +mkdir -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_1/RHB/Vaisala/ +cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level1_mwx/EUREC4A_RonBrown_Vaisala-RS_L1*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_1/RHB/Vaisala/ +mkdir -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_1/BCO/Vaisala/ +cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level1_mwx/EUREC4A_BCO_Vaisala-RS_L1*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_1/BCO/Vaisala/ +mkdir -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_1/ATL/MeteoModem/ +cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level1_bufr/EUREC4A_Atalante_Meteomodem-RS_L1*.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_1/ATL/MeteoModem/ + +# Export level2 data to export folder +mkdir -p /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/ +cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/EUREC4A_RonBrown_soundings.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_RonBrown_soundings.nc +cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/EUREC4A_Atalante_soundings_Vaisala.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_Atalante_soundings_Vaisala.nc +cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/EUREC4A_Meteor_soundings.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_Meteor_soundings.nc +cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/EUREC4A_BCO_soundings.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_BCO_soundings.nc +cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_mwx/EUREC4A_MS-Merian_soundings.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_MS-Merian_soundings.nc +cp /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v5/level2_bufr/EUREC4A_Atalante_soundings_Meteomodem.nc /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_Atalante_soundings_Meteomodem.nc + + +python ../postprocessing/remove_lower_levels.py + +mv /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_RonBrown_soundings.nc2 /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_RonBrown_soundings.nc +mv /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_Atalante_soundings_Vaisala.nc2 /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_Atalante_soundings_Vaisala.nc +mv /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_Meteor_soundings.nc2 /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_Meteor_soundings.nc +mv /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_BCO_soundings.nc2 /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_BCO_soundings.nc +mv /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_MS-Merian_soundings.nc2 /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_MS-Merian_soundings.nc +mv /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_Atalante_soundings_Meteomodem.nc2 /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_Atalante_soundings_Meteomodem.nc + +python ../postprocessing/rename_variable_names.py +python ../postprocessing/change_units.py -i '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/*.nc' +python ../postprocessing/change_units.py -i '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_1/*/*/*.nc' + +for file in `ls /mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_1/*/*/*.nc`: +do +ncatted -O -h -a coordinates,global,d,, $file +done + + +echo "Please include the version number in the filename of the level 2 data. This is not done automatically." diff --git a/postprocessing/change_units.py b/postprocessing/change_units.py new file mode 100644 index 0000000..dbe69be --- /dev/null +++ b/postprocessing/change_units.py @@ -0,0 +1,81 @@ +""" +Script to change units +""" +import os +import glob +import tqdm +import argparse +import numpy as np +import xarray as xr +import logging + +logging.basicConfig(filename='change_units.log', level=logging.INFO) + + +parser = argparse.ArgumentParser() +parser.add_argument('-i', '--inputfilefmt', required=True, + help='Input filename format (level 1 files), e.g. EUREC4A_*soundings_*.nc') +args = vars(parser.parse_args()) + +files = sorted(glob.glob(args['inputfilefmt'])) + +logging.info('Total number of soundings found: {}'.format(len(files))) + + +def hPa2Pa(hPa): + """ + Convert hPa to Pa + """ + assert np.all(hPa[~np.isnan(hPa)] < 1100), 'Is this really hPa?' + return hPa * 100. + + +def gkg2kgkg(gkg): + return gkg / 1000. + + +def degC2K(degC): + assert np.all(degC[~np.isnan(degC)] < 150), 'Is this really deg C?' + return degC + 273.15 + + +def percent2frac(percent): + return percent/100. + + +units_convert = {'p': {'units_old': 'hPa', + 'units_new': 'Pa', + 'converter': hPa2Pa}, + 'mr': {'units_old': 'g/kg', + 'units_new': 'kg/kg', + 'converter': gkg2kgkg}, + 'q': {'units_old': 'g/kg', + 'units_new': 'kg/kg', + 'converter': gkg2kgkg}, + 'dp': {'units_old': 'degrees_Celsius', + 'units_new': 'K', + 'converter': degC2K}, + 'ta': {'units_old': 'degrees_Celsius', + 'units_new': 'K', + 'converter': degC2K}, + 'rh': {'units_old': '%', + 'units_new': '1', + 'converter': percent2frac}, + 'dp': {'units_old': 'degrees_Celsius', + 'units_new': 'K', + 'converter': degC2K} + } + + +for f, file in enumerate(tqdm.tqdm(files)): + ds = xr.load_dataset(file) + for var in ds.data_vars: + if var in units_convert.keys(): + var_dict = units_convert[var] + if ds[var].attrs['units'] == var_dict['units_old']: + ds[var].values = var_dict['converter'](ds[var].values) + ds[var].attrs['units'] = var_dict['units_new'] + if 'coordinates' in ds.attrs: + del ds.attrs['coordinates'] + ds.to_netcdf(file, unlimited_dims=['sounding']) + diff --git a/postprocessing/remove_lower_levels.py b/postprocessing/remove_lower_levels.py index ab5c3a4..1308111 100644 --- a/postprocessing/remove_lower_levels.py +++ b/postprocessing/remove_lower_levels.py @@ -5,30 +5,33 @@ """ level2_files = [ - '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_Meteor_soundings.nc', - '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_BCO_soundings.nc', - '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_RH-Brown_soundings.nc', - '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_MS-Merian_soundings.nc', - '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_Atalante_soundings_Vaisala.nc', - '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_2/EUREC4A_Atalante_soundings_MeteoModem.nc' + '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_Meteor_soundings.nc', + '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_BCO_soundings.nc', + '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_RonBrown_soundings.nc', + '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_MS-Merian_soundings.nc', + '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_Atalante_soundings_Vaisala.nc', + '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_Atalante_soundings_Meteomodem.nc' ] +import os import numpy as np import xarray as xr +import tqdm -for file in level2_files: +for file in tqdm.tqdm(level2_files): print(file) - ds_in = xr.open_dataset(file) - if file == '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_v4/level2_mwx/EUREC4A_BCO_soundings.nc': + ds_in = xr.load_dataset(file) + #os.remove(file) + if file == '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EUREC4A_BCO_soundings.nc': ds_out = ds_in else: ds_ = ds_in lowest_mask = np.where(ds_in.altitude<40, True, False) for var in ds_in.data_vars: - if len(ds_in[var].dims) == 2: + if (len(ds_in[var].dims) == 2) and ('_FillValue' in ds_in[var].encoding.keys()): if ds_in[var].dims[0] == 'sounding': - ds_[var][0,lowest_mask] = np.nan - else: - ds_[var][lowest_mask,0] = np.nan + ds_[var].values[:,lowest_mask] = np.nan + elif ds_in[var].dims[1] == 'soudning': + ds_[var].values[lowest_mask,:] = np.nan ds_out = ds_ for var in ds_out.data_vars: @@ -40,4 +43,4 @@ del ds_out.attrs['python_version'] except: pass - ds_out.to_netcdf(file+'2') + ds_out.to_netcdf(file+'2', unlimited_dims=['sounding']) diff --git a/postprocessing/rename_variable_names.py b/postprocessing/rename_variable_names.py index c1320a1..21272cd 100644 --- a/postprocessing/rename_variable_names.py +++ b/postprocessing/rename_variable_names.py @@ -8,12 +8,13 @@ import numpy as np import xarray as xr -files_l2 = sorted(glob.glob('/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_2/EURE*.nc')) -files_l1 = sorted(glob.glob('/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_1/*/*/*.nc')) +files_l2 = sorted(glob.glob('/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_2/EURE*.nc')) +files_l1 = sorted(glob.glob('/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export2/level_1/*/*/*.nc')) rename_dict = { 'altitude': 'alt', - 'ascentRate': 'ascent_rate', + 'ascentRate': 'dz', + 'ascent_rate': 'dz', 'dewPoint': 'dp', 'flight_time': 'flight_time', 'humidity': 'rh', @@ -24,7 +25,7 @@ 'relative_humidity': 'rh', 'specific_humidity': 'q', 'pressure': 'p', - 'temperature': 'tc', + 'temperature': 'ta', 'windDirection': 'wdir', 'windSpeed':'wspd', 'wind_direction': 'wdir', @@ -34,9 +35,16 @@ 'wind_v': 'v' } -attrs_to_delete = ['surface_altitude', 'latitude_of_launch_location', 'longitude_of_launch_location', 'python_version', 'converted_by', 'contact_person', 'institution', 'location'] +rename_attrs_dict = { + 'git-version': 'version', + 'git_version': 'version' + } -files = np.hstack([files_l1, files_l2]) +attrs_to_delete = ['platform_location', 'surface_altitude', 'latitude_of_launch_location', 'longitude_of_launch_location', 'python_version', 'converted_by', 'contact_person', 'institution', 'location'] + +vars_to_delete = ['altitude_WGS84'] + +files = np.hstack([files_l2, files_l1]) for file in tqdm.tqdm(files): #if file == '/mnt/lustre02/work/mh0010/m300408/EUREC4Asoundings_export/level_1/ATL/Vaisala/EUREC4A_ATL_sounding_ascent_20200127_1059.nc': @@ -47,6 +55,11 @@ if var in rename_dict.keys(): ds = ds.rename({var: rename_dict[var]}) + for attr in list(ds.attrs.keys()): + if attr in rename_attrs_dict.keys(): + ds.attrs[rename_attrs_dict[attr]] = ds.attrs[attr] + del ds.attrs[attr] + for var in list(ds.variables): ds[var].encoding['zlib'] = True if '_FillValue' in ds[var].encoding.keys(): @@ -56,16 +69,19 @@ ds[var].encoding['_FillValue'] = 9.96921e36 if 'coordinates' in ds[var].attrs.keys(): del ds[var].attrs['coordinates'] - ds[var].encoding['coordinates'] = "flight_time lat lon" + ds[var].encoding['coordinates'] = "sounding_id flight_time lat lon" attrs = list(ds.attrs.keys()) for attr in attrs: if attr in attrs_to_delete: del ds.attrs[attr] + for var in vars_to_delete: + if var in ds.data_vars: + del ds[var] for var in list(ds.variables): if 'coordinates' in ds[var].encoding.keys(): - ds[var].encoding['coordinates'] = "flight_time lat lon" + ds[var].encoding['coordinates'] = "sounding_id flight_time lat lon" ds.flight_time.encoding['units'] = "seconds since 1970-01-01 00:00:00 UTC" ds.launch_time.encoding['units'] = "seconds since 1970-01-01 00:00:00 UTC" - ds.to_netcdf(file) + ds.to_netcdf(file, unlimited_dims=['sounding']) diff --git a/requirements.txt b/requirements.txt index ec6b2be..a5ddc4c 100644 --- a/requirements.txt +++ b/requirements.txt @@ -7,3 +7,4 @@ pillow>=5.0.0 eccodes>=2.13.0 metpy>=0.10.0 tqdm>=4.0.0 +bottleneck>=1.3.2 diff --git a/setup.py b/setup.py index 220a247..61c29cf 100644 --- a/setup.py +++ b/setup.py @@ -29,7 +29,7 @@ 'sounding_visualize=eurec4a_snd.make_quicklooks_rs41:main', 'sounding_skewT=eurec4a_snd.visualize.make_skewT_metpy:main', 'sounding_interpolate=eurec4a_snd.interpolate.batch_interpolate_soundings:main', - 'sounding_converter_mwx=eurec4a_snd.L1_mwx.py:main' + 'sounding_converter_mwx=eurec4a_snd.L1_mwx:main' ]}, - package_data={"eurec4a_snd": ["examples/data/*", "config/meta_information_template.ini"]} + package_data={"eurec4a_snd": ["examples/data/*", "config/meta_information_template.ini", "config/mwx_config.json"]} )