diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 6f1b59c18..cae0d55ec 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -146,3 +146,5 @@ add_subdirectory(atm) # gdas aerosol tests add_subdirectory(aero) + +add_subdirectory(marine) diff --git a/test/marine/CMakeLists.txt b/test/marine/CMakeLists.txt new file mode 100644 index 000000000..4dda602ca --- /dev/null +++ b/test/marine/CMakeLists.txt @@ -0,0 +1,161 @@ +# Copy the bkg files +# ------------------ +set(TESTDATA ${PROJECT_BINARY_DIR}/test/testdata ) + + +# Symlink test input yaml files +# ----------------------------- +# create testinput dir +set (TESTINPUT_DIR ${PROJECT_BINARY_DIR}/test/marine/testinput) + +file(MAKE_DIRECTORY ${TESTINPUT_DIR}) + +# symlink +foreach(FILENAME ${test_input}) + get_filename_component(filename ${FILENAME} NAME) + execute_process( COMMAND ${CMAKE_COMMAND} -E create_symlink + ${FILENAME} + ${TESTINPUT_DIR} + ) +endforeach(FILENAME) + +# install +install(FILES ${test_input} + DESTINATION "test/testinput/") + + +########################################################################### +# bufr to ioda tests: +########################################################################### + +# prepare a test.yaml file from test.yaml.in by replacing +# placeholder patterns __BUFRINPUTDIR__ and __IODAOUTPUTDIR__ and __OCEANBASIN__ +# with actual directory paths +function(CREATE_CONFIG_FILE + test_config_in + test_config_out + bufr_input_dir + ioda_output_dir + ocean_basin_file +) + file(READ "${test_config_in}" file_content) + string(REPLACE "__BUFRINPUTDIR__" "\"${bufr_input_dir}\"" temp_content "${file_content}") + string(REPLACE "__IODAOUTPUTDIR__" "\"${ioda_output_dir}\"" temp_content2 "${temp_content}") + string(REPLACE "__OCEANBASIN__" "\"${ocean_basin_file}\"" temp_content3 "${temp_content2}") + file(WRITE "${test_config_out}" "${temp_content3}") +endfunction() + + +set(TEST_WORKING_DIR ${PROJECT_BINARY_DIR}/test/marine) + +set(MARINE_BUFR2IODA_DIR ${PROJECT_SOURCE_DIR}/ush/ioda/bufr2ioda/marine) +set(MARINE_BUFR2IODA_DIR ${MARINE_BUFR2IODA_DIR}/b2i) +set(CONFIG_DIR ${PROJECT_SOURCE_DIR}/test/marine/testinput) +set(TESTREF_DIR ${PROJECT_SOURCE_DIR}/test/marine/testref) + + +function(CHECK_AND_SET_PATH PATH1 PATH2 RESULT_VAR) + # Check if PATH1 exists + if(EXISTS ${PATH1}) + set(${RESULT_VAR} ${PATH1} PARENT_SCOPE) + set(${RESULT_VAR}_EXISTS TRUE PARENT_SCOPE) + return() + endif() + + # Check if PATH2 exists + if(EXISTS ${PATH2}) + set(${RESULT_VAR} ${PATH2} PARENT_SCOPE) + set(${RESULT_VAR}_EXISTS TRUE PARENT_SCOPE) + return() + endif() + + # If neither path exists + set(${RESULT_VAR} "" PARENT_SCOPE) + set(${RESULT_VAR}_EXISTS FALSE PARENT_SCOPE) +endfunction() + + +# find the input files +set(BUFR_TEST_DIR_ORION + "/work/noaa/da/marineda/gfs-marine/data/obs/ci/bufr" +) +set(BUFR_TEST_DIR_HERA + "/scratch1/NCEPDEV/da/common/ci/bufr" +) +CHECK_AND_SET_PATH( + ${BUFR_TEST_DIR_ORION} + ${BUFR_TEST_DIR_HERA} + BUFR_TEST_DIR +) +if (NOT BUFR_TEST_DIR_EXISTS) + message(STATUS "BUFR test file directory not found -- bufr to ioda tests not generted.") + set(GENERATE_BUFR2IODA_TESTS FALSE) +else() + # message(STATUS "Found bufr test directory: ${BUFR_TEST_DIR}") + + set(OCEAN_BASIN_FILE_HERA + "/scratch2/NCEPDEV/ocean/Guillaume.Vernieres/data/static/common/RECCAP2_region_masks_all_v20221025.nc" + ) + set(OCEAN_BASIN_FILE_ORION + "/work/noaa/global/glopara/fix/gdas/soca/20240802/common/RECCAP2_region_masks_all_v20221025.nc" + ) + CHECK_AND_SET_PATH( + ${OCEAN_BASIN_FILE_ORION} + ${OCEAN_BASIN_FILE_HERA} + OCEAN_BASIN_FILE + ) + if (NOT OCEAN_BASIN_FILE_EXISTS) + message("Ocean basin data file not found -- bufr to ioda tests not generated.") + set(GENERATE_BUFR2IODA_TESTS FALSE) + endif() + # message(STATUS "Found ocean basin data in ${OCEAN_BASIN_FILE}") + set(GENERATE_BUFR2IODA_TESTS TRUE) +endif() + + +function(ADD_INSITU_TEST testname testbufr) + # set(CONFIG_TYPE "json") + set(CONFIG_TYPE "yaml") + set(DATE "2021063006") + set(TEST "bufr2ioda_insitu_${testname}") + + set(TESTREF_FILE "${TEST}_${DATE}.ref") + + # stage the input file to directory ${BUFR_INPUT_DIR} + set(BUFR_INPUT_DIR ${TEST_WORKING_DIR}) + set(BUFR_TEST_FILE "${DATE}-gdas.t06z.${testbufr}.tm00.bufr_d") + set(BUFR_FILE "${BUFR_TEST_DIR}/${BUFR_TEST_FILE}") + if (NOT EXISTS ${BUFR_FILE}) + message(FATAL_ERROR "BUFR file ${BUFR_FILE} not found") + endif() + file(COPY ${BUFR_FILE} DESTINATION ${BUFR_INPUT_DIR}) + + # stage the config file + set(CONFIG_FILE_NAME ${TEST}_${DATE}.${CONFIG_TYPE}) + set(CONFIG_FILE_IN "${CONFIG_DIR}/${CONFIG_FILE_NAME}.in") + set(CONFIG_FILE "${TEST_WORKING_DIR}/${CONFIG_FILE_NAME}") + set(IODA_OUTPUT_DIR ${TEST_WORKING_DIR}) + CREATE_CONFIG_FILE( + ${CONFIG_FILE_IN} + ${CONFIG_FILE} + ${BUFR_INPUT_DIR} + ${IODA_OUTPUT_DIR} + ${OCEAN_BASIN_FILE} + ) + + add_test( + NAME test_${TEST} + COMMAND ${MARINE_BUFR2IODA_DIR}/${TEST}.py -c ${CONFIG_FILE} -t ${TESTREF_DIR}/${TESTREF_FILE} + WORKING_DIRECTORY ${TEST_WORKING_DIR} + ) +endfunction() + + +if (GENERATE_BUFR2IODA_TESTS) + ADD_INSITU_TEST("profile_argo" "subpfl") + ADD_INSITU_TEST("profile_bathy" "bathy") + ADD_INSITU_TEST("profile_glider" "subpfl") + ADD_INSITU_TEST("profile_tesac" "tesac") + ADD_INSITU_TEST("profile_xbtctd" "xbtctd") + ADD_INSITU_TEST("surface_trkob" "trkob") +endif() diff --git a/test/marine/testinput/bufr2ioda_insitu_profile_argo_2021063006.yaml.in b/test/marine/testinput/bufr2ioda_insitu_profile_argo_2021063006.yaml.in new file mode 100644 index 000000000..54ef72855 --- /dev/null +++ b/test/marine/testinput/bufr2ioda_insitu_profile_argo_2021063006.yaml.in @@ -0,0 +1,13 @@ +--- +data_format: subpfl +subsets: SUBPFL +source: NCEP data tank +data_type: argo +cycle_type: gdas +cycle_datetime: '2021063006' +dump_directory: __BUFRINPUTDIR__ +ioda_directory: __IODAOUTPUTDIR__ +ocean_basin: __OCEANBASIN__ +data_description: 6-hrly in situ ARGO profiles +data_provider: U.S. NOAA + diff --git a/test/marine/testinput/bufr2ioda_insitu_profile_bathy_2021063006.yaml.in b/test/marine/testinput/bufr2ioda_insitu_profile_bathy_2021063006.yaml.in new file mode 100644 index 000000000..442d0a05a --- /dev/null +++ b/test/marine/testinput/bufr2ioda_insitu_profile_bathy_2021063006.yaml.in @@ -0,0 +1,13 @@ +--- +data_format: bathy +subsets: BATHY +source: NCEP data tank +data_type: bathy +cycle_type: gdas +cycle_datetime: '2021063006' +dump_directory: __BUFRINPUTDIR__ +ioda_directory: __IODAOUTPUTDIR__ +ocean_basin: __OCEANBASIN__ +data_description: 6-hrly in situ Bathythermal profiles +data_provider: U.S. NOAA + diff --git a/test/marine/testinput/bufr2ioda_insitu_profile_glider_2021063006.yaml.in b/test/marine/testinput/bufr2ioda_insitu_profile_glider_2021063006.yaml.in new file mode 100644 index 000000000..8deeff58d --- /dev/null +++ b/test/marine/testinput/bufr2ioda_insitu_profile_glider_2021063006.yaml.in @@ -0,0 +1,13 @@ +--- +data_format: subpfl +subsets: SUBPFL +source: NCEP data tank +data_type: glider +cycle_type: gdas +cycle_datetime: '2021063006' +dump_directory: __BUFRINPUTDIR__ +ioda_directory: __IODAOUTPUTDIR__ +ocean_basin: __OCEANBASIN__ +data_description: 6-hrly in situ GLIDER profiles +data_provider: U.S. NOAA + diff --git a/test/marine/testinput/bufr2ioda_insitu_profile_tesac_2021063006.yaml.in b/test/marine/testinput/bufr2ioda_insitu_profile_tesac_2021063006.yaml.in new file mode 100644 index 000000000..36e5e5ab1 --- /dev/null +++ b/test/marine/testinput/bufr2ioda_insitu_profile_tesac_2021063006.yaml.in @@ -0,0 +1,13 @@ +--- +data_format: tesac +subsets: TESAC +source: NCEP data tank +data_type: tesac +cycle_type: gdas +cycle_datetime: '2021063006' +dump_directory: __BUFRINPUTDIR__ +ioda_directory: __IODAOUTPUTDIR__ +ocean_basin: __OCEANBASIN__ +data_description: 6-hrly in situ TESAC profiles +data_provider: U.S. NOAA + diff --git a/test/marine/testinput/bufr2ioda_insitu_profile_xbtctd_2021063006.yaml.in b/test/marine/testinput/bufr2ioda_insitu_profile_xbtctd_2021063006.yaml.in new file mode 100644 index 000000000..3d0df09e2 --- /dev/null +++ b/test/marine/testinput/bufr2ioda_insitu_profile_xbtctd_2021063006.yaml.in @@ -0,0 +1,13 @@ +--- +data_format: xbtctd +subsets: XBTCTD +source: NCEP data tank +data_type: xbtctd +cycle_type: gdas +cycle_datetime: '2021063006' +dump_directory: __BUFRINPUTDIR__ +ioda_directory: __IODAOUTPUTDIR__ +ocean_basin: __OCEANBASIN__ +data_description: 6-hrly in situ XBT/XCTD profiles +data_provider: U.S. NOAA + diff --git a/test/marine/testinput/bufr2ioda_insitu_surface_trkob_2021063006.yaml.in b/test/marine/testinput/bufr2ioda_insitu_surface_trkob_2021063006.yaml.in new file mode 100644 index 000000000..977f2bcc0 --- /dev/null +++ b/test/marine/testinput/bufr2ioda_insitu_surface_trkob_2021063006.yaml.in @@ -0,0 +1,13 @@ +--- +data_format: trkob +subsets: TRACKOB +source: NCEP data tank +data_type: trackob +cycle_type: gdas +cycle_datetime: '2021063006' +dump_directory: __BUFRINPUTDIR__ +ioda_directory: __IODAOUTPUTDIR__ +ocean_basin: __OCEANBASIN__ +data_description: 6-hrly in situ TRACKOB surface +data_provider: U.S. NOAA + diff --git a/test/marine/testref/bufr2ioda_insitu_profile_argo_2021063006.ref b/test/marine/testref/bufr2ioda_insitu_profile_argo_2021063006.ref new file mode 100644 index 000000000..aeaa3a049 --- /dev/null +++ b/test/marine/testref/bufr2ioda_insitu_profile_argo_2021063006.ref @@ -0,0 +1,23 @@ +dateTime: 97150, int64 min, max = 1625022120, 1625047140 +dateTime hash = efcd92a143d8589288563e635318f2705900b0970e68a193eeffca3526cb4824 +rcptdateTime: 97150, int64 min, max = 1625029260, 1626386700 +rcptdateTime hash = da1886e087ccca11b77474c6cbc2434b4bdce86ea34f76c7a36b7d3f06602103 +lon: 97150, float32 min, max = -174.8090057373047, 170.48915100097656 +lon hash = e79a198bd3a8d0e72c3749ee84db8f5d6b6499a6c231540117c1270abff515ea +lat: 97150, float32 min, max = -59.77299880981445, 76.25086212158203 +lat hash = be568a690e0f32204326993b16d4f7863b4f77db6cbe5febf30523098c9186b3 +depth: 97150, float32 min, max = 0.0, 5700.7001953125 +depth hash = 805664f9daa328e4dce41c33fdadfdf026908d7f6f044b3079bde03f343fe79a +stationID: 97150, self.T_min) & (self.temp <= self.T_max) + + def SalinityFilter(self): + return (self.saln >= self.S_min) & (self.saln <= self.S_max) + + def filter(self): + pass + +########################################################################### + + def write_to_ioda_file(self, obsspace): + self.metadata.write_to_ioda_file(obsspace) + self.additional_vars.write_to_ioda_file(obsspace) + self.write_obs_value_t(obsspace) + self.write_obs_value_s(obsspace) + +########################################################################## + + def set_obs_from_query_result(self, r): + self.temp = r.get('temp', group_by='depth') + self.temp -= 273.15 + self.saln = r.get('saln', group_by='depth') + +########################################################################### + + def write_obs_value_t(self, obsspace): + obsspace.create_var('ObsValue/' + self.T_name, dtype=self.temp.dtype, fillval=self.temp.fill_value) \ + .write_attr('units', 'degC') \ + .write_attr('valid_range', np.array([self.T_min, self.T_max], dtype=np.float32)) \ + .write_attr('long_name', self.T_name) \ + .write_data(self.temp) + + def write_obs_value_s(self, obsspace): + obsspace.create_var( + 'ObsValue/' + self.S_name, + dtype=self.saln.dtype, + fillval=self.saln.fill_value + ) \ + .write_attr('units', 'psu') \ + .write_attr('valid_range', np.array([self.S_min, self.S_max], dtype=np.float32)) \ + .write_attr('long_name', self.S_name) \ + .write_data(self.saln) + +############################################################################## + + def log(self, logger): + self.metadata.log(logger) + self.log_obs(logger) + self.additional_vars.log(logger) + + def log_obs(self, logger): + self.log_temperature(logger) + self.log_salinity(logger) + + def log_temperature(self, logger): + log_variable(logger, "temp", self.temp) + logger.debug(f"temp hash = {compute_hash(self.temp)}") + + def log_salinity(self, logger): + log_variable(logger, "saln", self.saln) + logger.debug(f"saln hash = {compute_hash(self.saln)}") diff --git a/ush/ioda/bufr2ioda/marine/b2i/b2iconverter/ocean.py b/ush/ioda/bufr2ioda/marine/b2i/b2iconverter/ocean.py new file mode 100755 index 000000000..9e848c5f5 --- /dev/null +++ b/ush/ioda/bufr2ioda/marine/b2i/b2iconverter/ocean.py @@ -0,0 +1,131 @@ +#!/usr/bin/env python3 + +import os +import sys +import numpy as np +import numpy.ma as ma +import math +import matplotlib.pyplot as plt +import cartopy.crs as ccrs +import cartopy.feature as cfeature +import netCDF4 as nc +import xarray as xr + +# OceanBasin class provides a facility to add an OceanBasin +# metadata variable using lon and lat +# basic definition of ocean basins is read from an nc file, +# We search for the filename, depending on the system +# The path to the ocean basin nc file can be supplied +# in the implementation of the converter + +# the main method is get_station_basin which returns the ocean basin +# for a list of station coordinates +# there are methods for plotting and printing the ocean basin data +# as well as printing and plotting station basin data + + +class OceanBasin: + def __init__(self): + pass + + def set_ocean_basin_nc_file(self, filename): + self.ocean_basin_nc_file_path = filename + + def read_nc_file(self): + try: + with nc.Dataset(self.ocean_basin_nc_file_path, 'r') as nc_file: + variable_name = 'open_ocean' + if variable_name in nc_file.variables: + lat_dim = nc_file.dimensions['lat'].size + lon_dim = nc_file.dimensions['lon'].size + self.__latitudes = nc_file.variables['lat'][:] + self.__longitudes = nc_file.variables['lon'][:] + + variable = nc_file.variables[variable_name] + # Read the variable data into a numpy array + variable_data = variable[:] + # Convert to 2D numpy array + self.__basin_array = np.reshape(variable_data, (lat_dim, lon_dim)) + except FileNotFoundError: + print(f"The file {file_path} does not exist.") + sys.exit(1) + except IOError as e: + # Handle other I/O errors, such as permission errors + print(f"An IOError occurred: {e}") + sys.exit(1) + + def print_basin(self): + for i in range(n1): + for j in range(n2): + print(i, j, self.__basin_array[i][j]) + + def plot_basin(self): + # Create a figure and axes with Cartopy projection + fig = plt.figure(figsize=(10, 6)) + ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree()) + + # Plot the ocean basins using a colormap with 6 colors + # cmap = plt.cm.get_cmap('rainbow', 6) # Choose a colormap with 6 colors + cmap = plt.get_cmap('viridis', 6) # Create a colormap with 6 discrete colors + im = ax.pcolormesh(self.__longitudes, self.__latitudes, self.__basin_array, cmap='viridis', shading='auto', transform=ccrs.PlateCarree()) + + # Add colorbar + cbar = fig.colorbar(im, ax=ax, orientation='vertical', pad=0.05, ticks=np.arange(0, 6)) + cbar.set_label('Ocean Basin', fontsize=12) + # Add title and gridlines + ax.set_title('Ocean Basin Map', fontsize=16) + ax.coastlines() + ax.gridlines(draw_labels=True) + # Show the plot + plt.show() + plt.savefig('ocean_basin.png', dpi=300) + + # input: 2 vectors of station coordinates + # output: a vector of station ocean basin values + def get_station_basin(self, lat, lon): + n = len(lon) + # print("number of stations = ", n) + + lat0 = self.__latitudes[0] + dlat = self.__latitudes[1] - self.__latitudes[0] + lon0 = self.__longitudes[0] + dlon = self.__longitudes[1] - self.__longitudes[0] + + # the data may be a masked array + ocean_basin = [] + for i in range(n): + if not ma.is_masked(lat[i]): + i1 = round((lat[i] - lat0) / dlat) + i2 = round((lon[i] - lon0) / dlon) + ocean_basin.append(self.__basin_array[i1][i2]) + return ocean_basin + + def print_station_basin(self, lon, lat, file_path): + ocean_basin = self.get_station_basin(lat, lon) + with open(file_path, 'w') as file: + # Iterate over lon, lat, and ocean_basin arrays simultaneously + for lat_val, lon_val, basin_val in zip(lat, lon, ocean_basin): + file.write(f"{lat_val} {lon_val} {basin_val}\n") + + def plot_stations(self, lon, lat, png_file): + ocean_basin = self.get_station_basin(lon, lat) + + # Initialize the plot + plt.figure(figsize=(12, 8)) + # Create a Cartopy map with PlateCarree projection (latitude/longitude) + ax = plt.axes(projection=ccrs.PlateCarree()) + # Add coastlines and borders + ax.coastlines() + ax.add_feature(cartopy.feature.BORDERS, linestyle=':', linewidth=0.5) + + # Scatter plot with colored dots for each basin type + colors = ['blue', 'green', 'red', 'cyan', 'magenta', 'yellow'] + for basin_type in range(6): + indices = np.where(ocean_basin == basin_type)[0] + ax.scatter(lon[indices], lat[indices], color=colors[basin_type], label=f'Basin {basin_type}', alpha=0.7) + + # Add a legend + plt.legend(loc='lower left') + # Add title and show plot + plt.title('Ocean Basins Plot using Cartopy') + plt.savefig(png_file, dpi=300) diff --git a/ush/ioda/bufr2ioda/marine/b2i/b2iconverter/util.py b/ush/ioda/bufr2ioda/marine/b2i/b2iconverter/util.py new file mode 100644 index 000000000..5d3e039b3 --- /dev/null +++ b/ush/ioda/bufr2ioda/marine/b2i/b2iconverter/util.py @@ -0,0 +1,176 @@ +import os +import sys +import argparse +import subprocess +import numpy as np +import tempfile +import hashlib + + +def parse_arguments(): + parser = argparse.ArgumentParser() + parser.add_argument( + '-c', '--config', + type=str, + help='Input JSON or YAML configuration', required=True + ) + parser.add_argument( + '-l', '--log_file', + type=str, + help='Output file for testing ioda variables' + ) + parser.add_argument( + '-t', '--test', + type=str, + help='Input test reference file' + ) + args = parser.parse_args() + config_file = args.config + log_file = args.log_file + test_file = args.test + script_name = sys.argv[0] + return script_name, config_file, log_file, test_file + + +def log_variable(logger, v_name, v): + logger.debug(f"{v_name}: {len(v)}, {v.dtype} min, max = {v.min()}, {v.max()}") + + +def run_diff(file1, file2, logger): + try: + # Run the diff command + result = subprocess.run( + ['diff', file1, file2], + capture_output=True, text=True, check=False + ) + + # Check if diff command succeeded (return code 0) + if result.returncode == 0: + pass + elif result.returncode == 1: + logger.error("diff on files:") + logger.error(f"{file1}") + logger.error(f"{file2}") + logger.error("Files are different:") + logger.error(f"{result.stdout}") + else: + logger.error("Error occurred while running diff command.") + logger.error(f"{result.stdout}") + + except subprocess.CalledProcessError as e: + logger.error(f"Error occurred: {e}") + + return result.returncode + + +# use hash for testing; +def compute_hash(sequence, algorithm='sha256'): + """ + Compute a hash of the given sequence using the specified algorithm. + + :param sequence: A sequence of numbers (e.g., list of integers). + :param algorithm: The hash algorithm to use (e.g., 'sha256'). + :return: The hexadecimal digest of the hash. + """ + # Convert the sequence to a byte string + sequence_bytes = bytes(sequence) + # Create a hash object + hash_obj = hashlib.new(algorithm) + # Update the hash object with the byte string + hash_obj.update(sequence_bytes) + # Return the hexadecimal digest of the hash + return hash_obj.hexdigest() + + +##################################################################### + +def write_date_time(obsspace, dateTime): + obsspace.create_var( + 'MetaData/dateTime', + dtype=dateTime.dtype, fillval=dateTime.fill_value + ) \ + .write_attr('units', 'seconds since 1970-01-01T00:00:00Z') \ + .write_attr('long_name', 'Datetime') \ + .write_data(dateTime) + + +def write_rcpt_date_time(obsspace, rcptdateTime): + obsspace.create_var( + 'MetaData/rcptdateTime', + dtype=rcptdateTime.dtype, fillval=rcptdateTime.fill_value + ) \ + .write_attr('units', 'seconds since 1970-01-01T00:00:00Z') \ + .write_attr('long_name', 'receipt Datetime') \ + .write_data(rcptdateTime) + + +def write_longitude(obsspace, lon): + obsspace.create_var( + 'MetaData/longitude', + dtype=lon.dtype, fillval=lon.fill_value + ) \ + .write_attr('units', 'degrees_east') \ + .write_attr('valid_range', np.array([-180, 180], dtype=np.float32)) \ + .write_attr('long_name', 'Longitude') \ + .write_data(lon) + + +def write_latitude(obsspace, lat): + obsspace.create_var( + 'MetaData/latitude', + dtype=lat.dtype, fillval=lat.fill_value + ) \ + .write_attr('units', 'degrees_north') \ + .write_attr('valid_range', np.array([-90, 90], dtype=np.float32)) \ + .write_attr('long_name', 'Latitude') \ + .write_data(lat) + + +def write_station_id(obsspace, stationID): + obsspace.create_var( + 'MetaData/stationID', + dtype=stationID.dtype, fillval=stationID.fill_value + ) \ + .write_attr('long_name', 'Station Identification') \ + .write_data(stationID) + + +def write_depth(obsspace, depth): + obsspace.create_var( + 'MetaData/depth', + dtype=depth.dtype, + fillval=depth.fill_value + ) \ + .write_attr('units', 'm') \ + .write_attr('long_name', 'Water depth') \ + .write_data(depth) + + +def write_seq_num(obsspace, seqNum, datatype, fillvalue): + obsspace.create_var( + 'MetaData/sequenceNumber', + dtype=datatype, + fillval=fillvalue + ) \ + .write_attr('long_name', 'Sequence Number') \ + .write_data(seqNum) + + +def write_obs_error(obsspace, v_name, units, v): + obsspace.create_var( + v_name, + dtype=v.dtype, fillval=v.fill_value + ) \ + .write_attr('units', units) \ + .write_attr('long_name', 'ObsError') \ + .write_data(v) + + +def write_ocean_basin(obsspace, ocean_basin, datatype, fillvalue): + obsspace.create_var( + 'MetaData/oceanBasin', + dtype=datatype, + fillval=fillvalue + ) \ + .write_attr('long_name', 'Ocean basin') \ + .write_data(ocean_basin) diff --git a/ush/ioda/bufr2ioda/marine/b2i/bathy_ioda_variables.py b/ush/ioda/bufr2ioda/marine/b2i/bathy_ioda_variables.py new file mode 100644 index 000000000..ec49e39c0 --- /dev/null +++ b/ush/ioda/bufr2ioda/marine/b2i/bathy_ioda_variables.py @@ -0,0 +1,59 @@ +import numpy as np +from pyiodaconv import bufr +from b2iconverter.ioda_variables import IODAVariables +from b2iconverter.ioda_addl_vars import IODAAdditionalVariables, compute_seq_num + + +class BathyIODAVariables(IODAVariables): + def __init__(self): + self.construct() + self.additional_vars = BathyAdditionalVariables(self) + + def build_query(self): + q = super().build_query() + q.add('stationID', '*/RPID') + q.add('latitude', '*/CLAT') + q.add('longitude', '*/CLON') + q.add('depth', '*/BTOCN/DBSS') + q.add('temp', '*/BTOCN/STMP') + return q + + def set_obs_from_query_result(self, r): + self.temp = r.get('temp', group_by='depth') + self.temp -= 273.15 + + def filter(self): + mask = self.TemperatureFilter() + self.metadata.filter(mask) + self.temp = self.temp[mask] + + def write_to_ioda_file(self, obsspace): + self.metadata.write_to_ioda_file(obsspace) + self.additional_vars.write_to_ioda_file(obsspace) + self.write_obs_value_t(obsspace) + + def log_obs(self, logger): + self.log_temperature(logger) + + +class BathyAdditionalVariables(IODAAdditionalVariables): + + def construct(self): + self.seqNum = compute_seq_num(self.ioda_vars.metadata.lon, self.ioda_vars.metadata.lat) + n = len(self.seqNum) + self.PreQC = (np.ma.masked_array(np.full(n, 0))).astype(np.int32) + self.ObsError_temp = \ + np.float32(np.ma.masked_array(np.full(n, self.ioda_vars.T_error))) + self.compute_ocean_basin() + + def write_to_ioda_file(self, obsspace): + self.write_seq_num(obsspace) + self.write_preqc(obsspace, self.ioda_vars.T_name) + self.write_obs_errorT(obsspace) + self.write_ocean_basin(obsspace) + + def log(self, logger): + self.log_seq_num(logger) + self.log_preqc(logger) + self.log_obs_error_temp(logger) + self.log_ocean_basin(logger) diff --git a/ush/ioda/bufr2ioda/marine/b2i/bufr2ioda_insitu_profile_argo.py b/ush/ioda/bufr2ioda/marine/b2i/bufr2ioda_insitu_profile_argo.py new file mode 100755 index 000000000..bed8b6692 --- /dev/null +++ b/ush/ioda/bufr2ioda/marine/b2i/bufr2ioda_insitu_profile_argo.py @@ -0,0 +1,39 @@ +#!/usr/bin/env python3 + +import sys +from b2iconverter.util import parse_arguments +from b2iconverter.bufr2ioda_config import Bufr2iodaConfig +from b2iconverter.bufr2ioda_converter import Bufr2ioda_Converter +from argo_ioda_variables import ArgoIODAVariables + + +platform_description = 'ARGO profiles from subpfl: temperature and salinity' + + +class ArgoConfig(Bufr2iodaConfig): + def ioda_filename(self): + return f"{self.cycle_type}.t{self.hh}z.insitu_profile_argo.{self.cycle_datetime}.nc4" + + +if __name__ == '__main__': + + script_name, config_file, log_file, test_file = parse_arguments() + + bufr2ioda_config = ArgoConfig( + script_name, + config_file, + platform_description) + + ioda_vars = ArgoIODAVariables() + ioda_vars.set_temperature_var_name("waterTemperature") + ioda_vars.set_temperature_error(0.02) + ioda_vars.set_salinity_var_name("salinity") + ioda_vars.set_salinity_error(0.01) + + argo = Bufr2ioda_Converter(bufr2ioda_config, ioda_vars, log_file) + + argo.run() + + if test_file: + result = argo.test(test_file) + sys.exit(result) diff --git a/ush/ioda/bufr2ioda/marine/b2i/bufr2ioda_insitu_profile_bathy.py b/ush/ioda/bufr2ioda/marine/b2i/bufr2ioda_insitu_profile_bathy.py new file mode 100755 index 000000000..a492337a4 --- /dev/null +++ b/ush/ioda/bufr2ioda/marine/b2i/bufr2ioda_insitu_profile_bathy.py @@ -0,0 +1,32 @@ +#!/usr/bin/env python3 + +import sys +from b2iconverter.util import parse_arguments +from b2iconverter.bufr2ioda_config import Bufr2iodaConfig +from b2iconverter.bufr2ioda_converter import Bufr2ioda_Converter +from bathy_ioda_variables import BathyIODAVariables + + +platform_description = 'Profiles from BATHYthermal: temperature' + + +if __name__ == '__main__': + + script_name, config_file, log_file, test_file = parse_arguments() + + bufr2ioda_config = Bufr2iodaConfig( + script_name, + config_file, + platform_description) + + ioda_vars = BathyIODAVariables() + ioda_vars.set_temperature_error(0.24) + ioda_vars.set_temperature_var_name("waterTemperature") + + bathy = Bufr2ioda_Converter(bufr2ioda_config, ioda_vars, log_file) + + bathy.run() + + if test_file: + result = bathy.test(test_file) + sys.exit(result) diff --git a/ush/ioda/bufr2ioda/marine/b2i/bufr2ioda_insitu_profile_glider.py b/ush/ioda/bufr2ioda/marine/b2i/bufr2ioda_insitu_profile_glider.py new file mode 100755 index 000000000..627bf0292 --- /dev/null +++ b/ush/ioda/bufr2ioda/marine/b2i/bufr2ioda_insitu_profile_glider.py @@ -0,0 +1,33 @@ +#!/usr/bin/env python3 + +import sys +from b2iconverter.util import parse_arguments +from b2iconverter.bufr2ioda_config import Bufr2iodaConfig +from b2iconverter.bufr2ioda_converter import Bufr2ioda_Converter +from glider_ioda_variables import GliderIODAVariables + + +platform_description = 'GLIDER profiles from subpfl: temperature and salinity' + + +if __name__ == '__main__': + + script_name, config_file, log_file, test_file = parse_arguments() + + bufr2ioda_config = Bufr2iodaConfig( + script_name, + config_file, + platform_description) + + ioda_vars = GliderIODAVariables() + ioda_vars.set_temperature_var_name("waterTemperature") + ioda_vars.set_temperature_error(0.02) + ioda_vars.set_salinity_var_name("salinity") + ioda_vars.set_salinity_error(0.01) + + glider = Bufr2ioda_Converter(bufr2ioda_config, ioda_vars, log_file) + glider.run() + + if test_file: + result = glider.test(test_file) + sys.exit(result) diff --git a/ush/ioda/bufr2ioda/marine/b2i/bufr2ioda_insitu_profile_tesac.py b/ush/ioda/bufr2ioda/marine/b2i/bufr2ioda_insitu_profile_tesac.py new file mode 100755 index 000000000..b1d1d93e6 --- /dev/null +++ b/ush/ioda/bufr2ioda/marine/b2i/bufr2ioda_insitu_profile_tesac.py @@ -0,0 +1,34 @@ +#!/usr/bin/env python3 + +import sys +from b2iconverter.util import parse_arguments +from b2iconverter.bufr2ioda_config import Bufr2iodaConfig +from b2iconverter.bufr2ioda_converter import Bufr2ioda_Converter +from tesac_ioda_variables import TesacIODAVariables + + +platform_description = 'Profiles from TESAC: temperature and salinity' + + +if __name__ == '__main__': + + script_name, config_file, log_file, test_file = parse_arguments() + + bufr2ioda_config = Bufr2iodaConfig( + script_name, + config_file, + platform_description) + + ioda_vars = TesacIODAVariables() + ioda_vars.set_temperature_error(0.02) + ioda_vars.set_temperature_var_name("waterTemperature") + ioda_vars.set_salinity_error(0.01) + ioda_vars.set_salinity_var_name("salinity") + + tesac = Bufr2ioda_Converter(bufr2ioda_config, ioda_vars, log_file) + + tesac.run() + + if test_file: + result = tesac.test(test_file) + sys.exit(result) diff --git a/ush/ioda/bufr2ioda/marine/b2i/bufr2ioda_insitu_profile_xbtctd.py b/ush/ioda/bufr2ioda/marine/b2i/bufr2ioda_insitu_profile_xbtctd.py new file mode 100755 index 000000000..20fbead4e --- /dev/null +++ b/ush/ioda/bufr2ioda/marine/b2i/bufr2ioda_insitu_profile_xbtctd.py @@ -0,0 +1,33 @@ +#!/usr/bin/env python3 + +import sys +from b2iconverter.util import parse_arguments +from b2iconverter.bufr2ioda_config import Bufr2iodaConfig +from b2iconverter.bufr2ioda_converter import Bufr2ioda_Converter +from xbtctd_ioda_variables import XbtctdIODAVariables + + +platform_description = 'Profiles from XBT/CTD: temperature and salinity' + + +if __name__ == '__main__': + + script_name, config_file, log_file, test_file = parse_arguments() + + bufr2ioda_config = Bufr2iodaConfig( + script_name, + config_file, + platform_description) + + ioda_vars = XbtctdIODAVariables() + ioda_vars.set_temperature_var_name("waterTemperature") + ioda_vars.set_temperature_error(0.12) + ioda_vars.set_salinity_var_name("salinity") + ioda_vars.set_salinity_error(1.0) + + xbtctd = Bufr2ioda_Converter(bufr2ioda_config, ioda_vars, log_file) + xbtctd.run() + + if test_file: + result = xbtctd.test(test_file) + sys.exit(result) diff --git a/ush/ioda/bufr2ioda/marine/b2i/bufr2ioda_insitu_surface_trkob.py b/ush/ioda/bufr2ioda/marine/b2i/bufr2ioda_insitu_surface_trkob.py new file mode 100755 index 000000000..122c5cebe --- /dev/null +++ b/ush/ioda/bufr2ioda/marine/b2i/bufr2ioda_insitu_surface_trkob.py @@ -0,0 +1,38 @@ +#!/usr/bin/env python3 + +import sys +from b2iconverter.util import parse_arguments +from b2iconverter.bufr2ioda_config import Bufr2iodaConfig +from b2iconverter.bufr2ioda_converter import Bufr2ioda_Converter +from trkob_ioda_variables import TrkobIODAVariables + + +platform_description = 'Surface obs from TRACKOB: temperature and salinity' + + +class TrkobConfig(Bufr2iodaConfig): + def ioda_filename(self): + return f"{self.cycle_type}.t{self.hh}z.insitu_surface_{self.data_format}.{self.cycle_datetime}.nc4" + + +if __name__ == '__main__': + + script_name, config_file, log_file, test_file = parse_arguments() + + bufr2ioda_config = TrkobConfig( + script_name, + config_file, + platform_description) + + ioda_vars = TrkobIODAVariables() + ioda_vars.set_temperature_var_name("seaSurfaceTemperature") + ioda_vars.set_temperature_error(0.3) + ioda_vars.set_salinity_var_name("seaSurfaceSalinity") + ioda_vars.set_salinity_error(1.0) + + trkob = Bufr2ioda_Converter(bufr2ioda_config, ioda_vars, log_file) + trkob.run() + + if test_file: + result = trkob.test(test_file) + sys.exit(result) diff --git a/ush/ioda/bufr2ioda/marine/b2i/glider_ioda_variables.py b/ush/ioda/bufr2ioda/marine/b2i/glider_ioda_variables.py new file mode 100644 index 000000000..19f24ddc8 --- /dev/null +++ b/ush/ioda/bufr2ioda/marine/b2i/glider_ioda_variables.py @@ -0,0 +1,41 @@ +import numpy as np +from pyiodaconv import bufr +from b2iconverter.ioda_variables import IODAVariables + + +class GliderIODAVariables(IODAVariables): + def __init__(self): + super().__init__() + + def build_query(self): + q = super().build_query() + q.add('stationID', '*/WMOP') + q.add('latitude', '*/CLATH') + q.add('longitude', '*/CLONH') + q.add('depth', '*/GLPFDATA/WPRES') + q.add('temp', '*/GLPFDATA/SSTH') + q.add('saln', '*/GLPFDATA/SALNH') + return q + + def set_from_query_result(self, r): + super().set_from_query_result(r) + # convert depth in pressure units to meters (rho * g * h) + self.metadata.depth = np.float32(self.metadata.depth.astype(float) * 0.0001) + + def filter(self): + # Separate GLIDER profiles from subpfl tank + id = self.metadata.stationID + id_mask = (id >= 68900) & (id <= 68999) | \ + (id >= 1800000) & (id <= 1809999) | \ + (id >= 2800000) & (id <= 2809999) | \ + (id >= 3800000) & (id <= 3809999) | \ + (id >= 4800000) & (id <= 4809999) | \ + (id >= 5800000) & (id <= 5809999) | \ + (id >= 6800000) & (id <= 6809999) | \ + (id >= 7800000) & (id <= 7809999) + mask = self.TemperatureFilter() \ + & self.SalinityFilter() \ + & id_mask + self.metadata.filter(mask) + self.temp = self.temp[mask] + self.saln = self.saln[mask] diff --git a/ush/ioda/bufr2ioda/marine/b2i/tesac_ioda_variables.py b/ush/ioda/bufr2ioda/marine/b2i/tesac_ioda_variables.py new file mode 100644 index 000000000..c950b0c10 --- /dev/null +++ b/ush/ioda/bufr2ioda/marine/b2i/tesac_ioda_variables.py @@ -0,0 +1,27 @@ +from pyiodaconv import bufr +from b2iconverter.ioda_variables import IODAVariables + + +class TesacIODAVariables(IODAVariables): + def __init__(self): + super().__init__() + + def build_query(self): + q = super().build_query() + q.add('stationID', '*/RPID') + q.add('latitude', '*/CLAT') + q.add('longitude', '*/CLON') + q.add('depth', '*/BTOCN/DBSS') + q.add('temp', '*/BTOCN/STMP') + q.add('saln', '*/BTOCN/SALN') + return q + + def filter(self): + # Separate TESAC profiles tesac tank + id_mask = [True if id.isdigit() and id != 0 else False for id in self.metadata.stationID] + mask = id_mask \ + & self.TemperatureFilter() \ + & self.SalinityFilter() + self.metadata.filter(mask) + self.temp = self.temp[mask] + self.saln = self.saln[mask] diff --git a/ush/ioda/bufr2ioda/marine/b2i/trkob_ioda_variables.py b/ush/ioda/bufr2ioda/marine/b2i/trkob_ioda_variables.py new file mode 100644 index 000000000..7eeb3047c --- /dev/null +++ b/ush/ioda/bufr2ioda/marine/b2i/trkob_ioda_variables.py @@ -0,0 +1,87 @@ +import numpy as np +from pyiodaconv import bufr +from b2iconverter.ioda_variables import IODAVariables +from b2iconverter.ioda_addl_vars import IODAAdditionalVariables +from b2iconverter.ioda_metadata import IODAMetadata +from b2iconverter.util import write_date_time, write_rcpt_date_time, write_longitude, write_latitude, write_station_id + + +class TrkobIODAVariables(IODAVariables): + def __init__(self): + self.construct() + self.metadata = TrkobMetadata() + self.additional_vars = TrkobAdditionalVariables(self) + + def build_query(self): + q = super().build_query() + q.add('stationID', '*/RPID') + q.add('latitude', '*/CLAT') + q.add('longitude', '*/CLON') + q.add('depth', '*/BTOCN/DBSS') + q.add('temp', '*/BTOCN/STMP') + q.add('saln', '*/BTOCN/SALN') + return q + + def filter(self): + mask = self.TemperatureFilter() \ + & self.SalinityFilter() + self.temp = self.temp[mask] + self.saln = self.saln[mask] + self.metadata.filter(mask) + + +class TrkobMetadata(IODAMetadata): + + def set_from_query_result(self, r): + self.set_date_time_from_query_result(r) + self.set_rcpt_date_time_from_query_result(r) + self.set_lon_from_query_result(r) + self.set_lat_from_query_result(r) + self.set_station_id_from_query_result(r) + + def filter(self, mask): + self.dateTime = self.dateTime[mask] + self.rcptdateTime = self.rcptdateTime[mask] + self.lat = self.lat[mask] + self.lon = self.lon[mask] + self.stationID = self.stationID[mask] + + def write_to_ioda_file(self, obsspace): + write_date_time(obsspace, self.dateTime) + write_rcpt_date_time(obsspace, self.rcptdateTime) + write_longitude(obsspace, self.lon) + write_latitude(obsspace, self.lat) + write_station_id(obsspace, self.stationID) + + def log(self, logger): + self.log_date_time(logger) + self.log_rcpt_date_time(logger) + # self.logLonLat(logger) + self.log_longitude(logger) + self.log_latitude(logger) + self.log_station_id(logger) + + +class TrkobAdditionalVariables(IODAAdditionalVariables): + + def construct(self): + n = len(self.ioda_vars.temp) + self.PreQC = (np.ma.masked_array(np.full(n, 0))).astype(np.int32) + self.ObsError_temp = \ + np.float32(np.ma.masked_array(np.full(n, self.ioda_vars.T_error))) + self.ObsError_saln = \ + np.float32(np.ma.masked_array(np.full(n, self.ioda_vars.S_error))) + self.compute_ocean_basin() + + def write_to_ioda_file(self, obsspace): + self.write_preqc(obsspace, self.ioda_vars.T_name) + self.write_preqc(obsspace, self.ioda_vars.S_name) + self.write_obs_errorT(obsspace) + self.write_obs_errorS(obsspace) + self.write_ocean_basin(obsspace) + + def log(self, logger): + self.log_preqc(logger) + self.log_obs_error_temp(logger) + self.log_obs_error_saln(logger) + self.log_ocean_basin(logger) diff --git a/ush/ioda/bufr2ioda/marine/b2i/xbtctd_ioda_variables.py b/ush/ioda/bufr2ioda/marine/b2i/xbtctd_ioda_variables.py new file mode 100644 index 000000000..fc85ace38 --- /dev/null +++ b/ush/ioda/bufr2ioda/marine/b2i/xbtctd_ioda_variables.py @@ -0,0 +1,25 @@ +import numpy as np +from pyiodaconv import bufr +from b2iconverter.ioda_variables import IODAVariables + + +class XbtctdIODAVariables(IODAVariables): + def __init__(self): + super().__init__() + + def build_query(self): + q = super().build_query() + q.add('stationID', '*/WMOP') + q.add('latitude', '*/CLATH') + q.add('longitude', '*/CLONH') + q.add('depth', '*/TMSLPFSQ/DBSS') + q.add('temp', '*/TMSLPFSQ/SST1') + q.add('saln', '*/TMSLPFSQ/SALNH') + return q + + def filter(self): + mask = self.TemperatureFilter() \ + & self.SalinityFilter() + self.temp = self.temp[mask] + self.saln = self.saln[mask] + self.metadata.filter(mask)