diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..e7739b1 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +.DS_Store +.idea +.nc diff --git a/bufr_adpupa_prepbufr.py b/bufr_adpupa_prepbufr.py deleted file mode 100644 index 089ab4e..0000000 --- a/bufr_adpupa_prepbufr.py +++ /dev/null @@ -1,297 +0,0 @@ -#!/usr/bin/env python3 -import os -import sys -import bufr -import argparse -import copy -import numpy as np -import numpy.ma as ma -import math -import calendar -import time -from datetime import datetime -from pyioda.ioda.Engines.Bufr import Encoder as iodaEncoder -from bufr.encoders.netcdf import Encoder as netcdfEncoder -from wxflow import Logger - -# Initialize Logger -# Get log level from the environment variable, default to 'INFO it not set -log_level = os.getenv('LOG_LEVEL', 'INFO') -logger = Logger('bufr_adpupa_prepbufr.py', level=log_level, colored_log=False) - - -def logging(comm, level, message): - """ - Logs a message to the console or log file, based on the specified logging level. - - This function ensures that logging is only performed by the root process (`rank 0`) - in a distributed computing environment. The function maps the logging level to - appropriate logger methods and defaults to the 'INFO' level if an invalid level is provided. - - Parameters: - comm: object - The communicator object, typically from a distributed computing framework - (e.g., MPI). It must have a `rank()` method to determine the process rank. - level: str - The logging level as a string. Supported levels are: - - 'DEBUG' - - 'INFO' - - 'WARNING' - - 'ERROR' - - 'CRITICAL' - If an invalid level is provided, a warning will be logged, and the level - will default to 'INFO'. - message: str - The message to be logged. - - Behavior: - - Logs messages only on the root process (`comm.rank() == 0`). - - Maps the provided logging level to a method of the logger object. - - Defaults to 'INFO' and logs a warning if an invalid logging level is given. - - Supports standard logging levels for granular control over log verbosity. - - Example: - >>> logging(comm, 'DEBUG', 'This is a debug message.') - >>> logging(comm, 'ERROR', 'An error occurred!') - - Notes: - - Ensure that a global `logger` object is configured before using this function. - - The `comm` object should conform to MPI-like conventions (e.g., `rank()` method). - """ - - if comm.rank() == 0: - # Define a dictionary to map levels to logger methods - log_methods = { - 'DEBUG': logger.debug, - 'INFO': logger.info, - 'WARNING': logger.warning, - 'ERROR': logger.error, - 'CRITICAL': logger.critical, - } - - # Get the appropriate logging method, default to 'INFO' - log_method = log_methods.get(level.upper(), logger.info) - - if log_method == logger.info and level.upper() not in log_methods: - # Log a warning if the level is invalid - logger.warning(f'log level = {level}: not a valid level --> set to INFO') - - # Call the logging method - log_method(message) - - -def _compute_datetime(cycleTimeSinceEpoch, hrdr): - """ - Compute dateTime using the cycleTimeSinceEpoch and Observation Time Minus Cycle Time - - Parameters: - cycleTimeSinceEpoch: Time of cycle in Epoch Time - hrdr: Observation Time Minus Cycle Time - - Returns: - Masked array of dateTime values - """ - - int64_fill_value = np.int64(0) - dateTime = np.zeros(hrdr.shape, dtype=np.int64) - for i in range(len(dateTime)): - if ma.is_masked(hrdr[i]): - continue - else: - dateTime[i] = np.int64(hrdr[i]*3600) + cycleTimeSinceEpoch - - dateTime = ma.array(dateTime) - dateTime = ma.masked_values(dateTime, int64_fill_value) - - return dateTime - - -def _make_description(mapping_path, cycle_time, update=False): - description = bufr.encoders.Description(mapping_path) - - reference_time = np.int64(calendar.timegm(time.strptime(str(int(cycle_time)), '%Y%m%d%H'))) - - if update: - # Define the variables to be added in a list of dictionaries - variables = [ - { - 'name': 'MetaData/sequenceNumber', - 'source': 'variables/sequenceNumber', - 'units': '1', - 'longName': 'Sequence Number (Obs Subtype)', - }, - ] - - # Loop through each variable and add it to the description - for var in variables: - description.add_variable( - name=var['name'], - source=var['source'], - units=var['units'], - longName=var['longName'] - ) - - #description.add_global(name='datetimeReference', value=str(reference_time)) - - return description - - -def _make_obs(comm, input_path, mapping_path, cycle_time): - """ - Create the ioda adpupa prepbufr observations: - - reads values - - adds sequenceNum - - Parameters - ---------- - comm: object - The communicator object (e.g., MPI) - input_path: str - The input bufr file - mapping_path: str - The input bufr2ioda mapping file - cycle_time: str - The cycle in YYYYMMDDHH format - """ - - # Get container from mapping file first - logging(comm, 'INFO', 'Get container from bufr') - container = bufr.Parser(input_path, mapping_path).parse(comm) - - logging(comm, 'DEBUG', f'container list (original): {container.list()}') - logging(comm, 'DEBUG', f'prepbufrDataLevelCategory') - cat = container.get('variables/prepbufrDataLevelCategory') - - logging(comm, 'DEBUG', f'Do DateTime calculation') - hrdr = container.get('variables/obsTimeMinusCycleTime') - hrdr2 = np.array(hrdr) - cycleTimeSinceEpoch = np.int64(calendar.timegm(time.strptime(str(int(cycle_time)), '%Y%m%d%H'))) - dateTime = _compute_datetime(cycleTimeSinceEpoch, hrdr2) - logging(comm, 'DEBUG', f'dateTime min/max = {dateTime.min()} {dateTime.max()}') - - logging(comm, 'DEBUG', f'Change longitude range from [0,360] to [-180,180]') - lon = container.get('variables/longitude') - lon_paths = container.get_paths('variables/longitude') - lon[lon>180] -= 360 - - logging(comm, 'DEBUG', f'Make an array of 0s for MetaData/sequenceNumber') - sequenceNum = np.zeros(lon.shape, dtype=np.int32) - logging(comm, 'DEBUG', f' sequenceNummin/max = {sequenceNum.min()} {sequenceNum.max()}') - - logging(comm, 'DEBUG', f'Do ps calculation') - pob = container.get('variables/pressure') - ps = np.full(pob.shape[0], pob.fill_value) - ps = np.where(cat == 0, pob, ps) - - logging(comm, 'DEBUG', f'Do tsen and tv calculation') - tpc = container.get('variables/temperatureEventProgramCode') - tob = container.get('variables/airTemperature') - tsen = np.full(tob.shape[0], tob.fill_value) - tsen = np.where(((tpc >= 1) & (tpc < 8)), tob, tsen) - tvo = np.full(tob.shape[0], tob.fill_value) - tvo = np.where((tpc == 8), tob, tvo) - - logging(comm, 'DEBUG', f'Do ps QM calculations') - pqm = container.get('variables/pressureQualityMarker') - psqm = np.full(pqm.shape[0], pqm.fill_value) - psqm = np.where(cat == 0, pqm, psqm) - - logging(comm, 'DEBUG', f'Do tsen and tv QM calculations') - tobqm = container.get('variables/airTemperatureQualityMarker') - tsenqm = np.full(tobqm.shape[0], tobqm.fill_value) - tsenqm = np.where(((tpc >= 1) & (tpc < 8)), tobqm, tsenqm) - tvoqm = np.full(tobqm.shape[0], tobqm.fill_value) - tvoqm = np.where((tpc == 8), tobqm, tvoqm) - - logging(comm, 'DEBUG', f'Do ps ObsError calculations') - poe = container.get('variables/pressureError') - psoe = np.full(poe.shape[0], poe.fill_value) - psoe = np.where(cat == 0, poe, psoe) - - logging(comm, 'DEBUG', f'Do tsen and tv ObsError calculations') - toboe = container.get('variables/airTemperatureError') - tsenoe = np.full(toboe.shape[0], toboe.fill_value) - tsenoe = np.where(((tpc >= 1) & (tpc < 8)), toboe, tsenoe) - tvooe = np.full(toboe.shape[0], toboe.fill_value) - tvooe = np.where((tpc == 8), toboe, tvooe) - - logging(comm, 'DEBUG', f'Update variables in container') - container.replace('variables/longitude', lon) - container.replace('variables/timestamp', dateTime) - container.replace('variables/airTemperature', tsen) - container.replace('variables/airTemperatureQualityMarker', tsenqm) - container.replace('variables/airTemperatureError', tsenoe) - container.replace('variables/virtualTemperature', tvo) - container.replace('variables/virtualTemperatureQualityMarker', tvoqm) - container.replace('variables/virtualTemperatureError', tvooe) - container.replace('variables/stationPressure', ps) - container.replace('variables/stationPressureQualityMarker', psqm) - container.replace('variables/stationPressureError', psoe) - - logging(comm, 'DEBUG', f'Add variables to container') - container.add('variables/sequenceNumber', sequenceNum, lon_paths) - - # Check - logging(comm, 'DEBUG', f'container list (updated): {container.list()}') - - return container - - -def create_obs_group(input_path, mapping_path, cycle_time, env): - - comm = bufr.mpi.Comm(env["comm_name"]) - - logging(comm, 'INFO', f'Make description and make obs') - description = _make_description(mapping_path, cycle_time, update=True) - container = _make_obs(comm, input_path, mapping_path, cycle_time) - - # Gather data from all tasks into all tasks. Each task will have the complete record - logging(comm, 'INFO', f'Gather data from all tasks into all tasks') - container.all_gather(comm) - - logging(comm, 'INFO', f'Encode the data') - data = next(iter(iodaEncoder(description).encode(container).values())) - - logging(comm, 'INFO', f'Return the encoded data.') - return data - - -def create_obs_file(input_path, mapping_path, output_path, cycle_time): - - comm = bufr.mpi.Comm("world") - container = _make_obs(comm, input_path, mapping_path, cycle_time) - container.gather(comm) - - description = _make_description(mapping_path, cycle_time, update=True) - - # Encode the data - if comm.rank() == 0: - netcdfEncoder(description).encode(container, output_path) - - logging(comm, 'INFO', f'Return the encoded data') - - -if __name__ == '__main__': - start_time = time.time() - - bufr.mpi.App(sys.argv) - comm = bufr.mpi.Comm("world") - - # Required input arguments as positional arguments - parser = argparse.ArgumentParser(description="Convert BUFR to NetCDF using a mapping file.") - parser.add_argument('input', type=str, help='Input BUFR file') - parser.add_argument('mapping', type=str, help='BUFR2IODA Mapping File') - parser.add_argument('output', type=str, help='Output NetCDF file') - parser.add_argument('cycle_time', type=str, help='cycle time in YYYYMMDDHH format') - - args = parser.parse_args() - infile = args.input - mapping = args.mapping - output = args.output - cycle_time = args.cycle_time - - create_obs_file(infile, mapping, output, cycle_time) - - end_time = time.time() - running_time = end_time - start_time - logging(comm, 'INFO', f'Total running time: {running_time}') diff --git a/bufr_adpupa_prepbufr_mapping.yaml b/bufr_adpupa_prepbufr_mapping.yaml deleted file mode 100644 index dced147..0000000 --- a/bufr_adpupa_prepbufr_mapping.yaml +++ /dev/null @@ -1,371 +0,0 @@ -# (C) Copyright 2024 NOAA/NWS/NCEP/EMC - -bufr: - group_by_variable: prepbufrDataLevelCategory - subsets: - - ADPUPA - variables: - # ObsType - observationType: - query: "*/TYP" - - # MetaData - timestamp: - timeoffset: - timeOffset: "*/PRSLEVEL/DRFTINFO/HRDR" - transforms: - - scale: 3600 - referenceTime: "2021-08-01T00:00:00Z" - - prepbufrDataLevelCategory: - query: "*/PRSLEVEL/CAT" - - obsTimeMinusCycleTime: - query: "*/PRSLEVEL/DRFTINFO/HRDR" - - longitude: - query: "*/PRSLEVEL/DRFTINFO/XDR" - - latitude: - query: "*/PRSLEVEL/DRFTINFO/YDR" - - stationIdentification: - query: "*/SID" - - stationElevation: - query: "*/ELV" - - temperatureEventProgramCode: - query: "*/PRSLEVEL/T___INFO/T__EVENT{1}/TPC" - - pressure: - query: "*/PRSLEVEL/P___INFO/P__EVENT{1}/POB" - type: float - transforms: - - scale: 100 - - height: - query: "*/PRSLEVEL/Z___INFO/Z__EVENT{1}/ZOB" - type: float - - # ObsValue - stationPressure: - query: "*/PRSLEVEL/P___INFO/P__EVENT{1}/POB" - type: float - transforms: - - scale: 100 - - airTemperature: - query: "*/PRSLEVEL/T___INFO/T__EVENT{1}/TOB" - transforms: - - offset: 273.15 - - dewPointTemperature: - query: "*/PRSLEVEL/Q___INFO/TDO" - transforms: - - offset: 273.15 - - virtualTemperature: - query: "*/PRSLEVEL/T___INFO/TVO" - transforms: - - offset: 273.15 - - windEastward: - query: "*/PRSLEVEL/W___INFO/W__EVENT{1}/UOB" - - windNorthward: - query: "*/PRSLEVEL/W___INFO/W__EVENT{1}/VOB" - - specificHumidity: - query: "*/PRSLEVEL/Q___INFO/Q__EVENT{1}/QOB" - type: float - transforms: - - scale: 0.000001 - - # QualityMark - pressureQualityMarker: - query: "*/PRSLEVEL/P___INFO/P__EVENT{1}/PQM" - - stationPressureQualityMarker: - query: "*/PRSLEVEL/P___INFO/P__EVENT{1}/PQM" - - heightQualityMarker: - query: "*/PRSLEVEL/Z___INFO/Z__EVENT{1}/ZQM" - - airTemperatureQualityMarker: - query: "*/PRSLEVEL/T___INFO/T__EVENT{1}/TQM" - - virtualTemperatureQualityMarker: - query: "*/PRSLEVEL/T___INFO/T__EVENT{1}/TQM" - - windQualityMarker: - query: "*/PRSLEVEL/W___INFO/W__EVENT{1}/WQM" - - specificHumidityQualityMarker: - query: "*/PRSLEVEL/Q___INFO/Q__EVENT{1}/QQM" - - # ObsError - pressureError: - query: "*/PRSLEVEL/P___INFO/P__BACKG/POE" - transforms: - - scale: 100 - - stationPressureError: - query: "*/PRSLEVEL/P___INFO/P__BACKG/POE" - transforms: - - scale: 100 - - airTemperatureError: - query: "*/PRSLEVEL/T___INFO/T__BACKG/TOE" - transforms: - - offset: 273.15 - - virtualTemperatureError: - query: "*/PRSLEVEL/T___INFO/T__BACKG/TOE" - transforms: - - offset: 273.15 - - windError: - query: "*/PRSLEVEL/W___INFO/W__BACKG/WOE" - - specificHumidityError: - query: "*/PRSLEVEL/Q___INFO/Q__BACKG/QOE" - - -encoder: - type: netcdf - - globals: - - name: "dataOriginalFormatSpec" - type: string - value: "prepbufr" - - - name: "platforms" - type: string - value: "ADPUPA" - - - name: "source" - type: string - value: "prepBUFR" - - - name: "description" - type: string - value: "UPPER-AIR (RAOB, PIBAL, RECCO, DROPS) REPORTS" - - - name: "dataProviderOrigin" - type: string - value: "U.S. NOAA NCEP" - - variables: - # Observation Type: stationPressure - - name: "ObsType/stationPressure" - source: variables/observationType - longName: "Observation Type" - - # Observation Type: airTemperature - - name: "ObsType/airTemperature" - source: variables/observationType - longName: "Observation Type" - - # Observation Type: virtualTemperature - - name: "ObsType/virtualTemperature" - source: variables/observationType - longName: "Observation Type" - - # Observation Type: windEastward - - name: "ObsType/windEastward" - source: variables/observationType - longName: "Observation Type" - - # Observation Type: windNorthward - - name: "ObsType/windNorthward" - source: variables/observationType - longName: "Observation Type" - - # Observation Type: specificHumidity - - name: "ObsType/specificHumidity" - source: variables/observationType - longName: "Observation Type" - - # MetaData: dateTime - - name: "MetaData/dateTime" - source: variables/timestamp - longName: "Datetime" - units: "seconds since 1970-01-01T00:00:00Z" - - # MetaData: stationIdentification - - name: "MetaData/stationIdentification" - source: variables/stationIdentification - longName: "Station ID" - - # MetaData: longitude - - name: "MetaData/longitude" - source: variables/longitude - longName: "Longitude" - units: "degrees_east" - range: [0, 360] - - # MetaData: latitude - - name: "MetaData/latitude" - source: variables/latitude - longName: "Latitude" - units: "degrees_north" - range: [-90, 90] - - # MetaData: stationElevation - - name: "MetaData/stationElevation" - source: variables/stationElevation - longName: "Height of Station" - units: "m" - - # MetaData: temperatureEventProgramCode - - name: "MetaData/temperatureEventProgramCode" - source: variables/temperatureEventProgramCode - longName: "Temperature Event Code" - - # MetaData: prepbufrDataLevelCategory - - name: "MetaData/prepbufrDataLevelCategory" - source: variables/prepbufrDataLevelCategory - longName: "Prepbufr Data Level Category" - - # MetaData: pressure - - name: "MetaData/pressure" - source: variables/pressure - longName: "Pressure" - units: "Pa" - - # MetaData: height - - name: "MetaData/height" - source: variables/height - longName: "Height of Observation" - units: "m" - - # QCFlags: qualityFlags - - name: "QCFlags/qualityFlags" - #- name: "MetaData/temperatureEventProgramCode" - source: variables/temperatureEventProgramCode - longName: "Temperature Event Program Code" - - # ObsValue: stationPressure - - name: "ObsValue/stationPressure" - source: variables/stationPressure - longName: "Station Pressure" - units: "Pa" - - # ObsValue: airTemperature - - name: "ObsValue/airTemperature" - source: variables/airTemperature - longName: "Temperature" - units: "K" - - # ObsValue: dewPointTemperature - - name: "ObsValue/dewPointTemperature" - source: variables/dewPointTemperature - longName: "DewPoint Temperature" - units: "K" - - # ObsValue: virtualTemperature - - name: "ObsValue/virtualTemperature" - source: variables/virtualTemperature - longName: "Virtual Temperature Non-Q Controlled" - units: "K" - - # ObsValue: windEastward - - name: "ObsValue/windEastward" - source: variables/windEastward - longName: "Eastward Wind" - units: "m s-1" - - # ObsValue: windNorthward - - name: "ObsValue/windNorthward" - source: variables/windNorthward - longName: "Northward Wind" - units: "m s-1" - - # ObsValue: specificHumidity - - name: "ObsValue/specificHumidity" - source: variables/specificHumidity - longName: "Specific Humidity" - units: "kg kg-1" - - # QualityMarker: pressure - - name: "QualityMarker/pressure" - source: variables/pressureQualityMarker - longName: "Pressure Quality Marker" - - # QualityMarker: height - - name: "QualityMarker/height" - source: variables/heightQualityMarker - longName: "Height Quality Marker" - - # QualityMarker: stationPressure - - name: "QualityMarker/stationPressure" - source: variables/stationPressureQualityMarker - longName: "Station Pressure Quality Marker" - - # QualityMarker: airTemperature - - name: "QualityMarker/airTemperature" - source: variables/airTemperatureQualityMarker - longName: "Temperature Quality Marker" - - # QualityMarker: virtualTemperature - - name: "QualityMarker/virtualTemperature" - source: variables/virtualTemperatureQualityMarker - longName: "Virtual Temperature Quality Marker" - - # QualityMarker: windEastward - - name: "QualityMarker/windEastward" - source: variables/windQualityMarker - longName: "Eastward Wind Quality Marker" - - # QualityMarker: windNorthward - - name: "QualityMarker/windNorthward" - source: variables/windQualityMarker - longName: "Northward Wind Quality Marker" - - # QualityMarker: specificHumidity - - name: "QualityMarker/specificHumidity" - source: variables/specificHumidityQualityMarker - longName: "Specific Humidity Quality Marker" - - # ObsError: pressure - - name: "ObsError/pressure" - source: variables/pressureError - longName: "Pressure Error" - units: "Pa" - - # ObsError: stationPressure - - name: "ObsError/stationPressure" - source: variables/stationPressureError - longName: "Station Pressure Error" - units: "Pa" - - # ObsError: airTemperature - - name: "ObsError/airTemperature" - source: variables/airTemperatureError - longName: "Temperature Error" - units: "K" - - # ObsError: virtualTemperature - - name: "ObsError/virtualTemperature" - source: variables/virtualTemperatureError - longName: "Virtual Temperature Error" - units: "K" - - # ObsError: windEastward - - name: "ObsError/windEastward" - source: variables/windError - longName: "East Wind Error" - units: "m s-1" - - # ObsError: windNorthward - - name: "ObsError/windNorthward" - source: variables/windError - longName: "North Wind Error" - units: "m s-1" - - # ObsError: specificHumidity - - name: "ObsError/specificHumidity" - source: variables/specificHumidityError - longName: "Specific Humidity Error" - units: "kg kg-1" diff --git a/dump/mapping/__init__.py b/dump/mapping/__init__.py new file mode 100644 index 0000000..5e88191 --- /dev/null +++ b/dump/mapping/__init__.py @@ -0,0 +1 @@ +import builder diff --git a/dump/mapping/bufr_adpsfc_prepbufr.py b/dump/mapping/bufr_adpsfc_prepbufr.py deleted file mode 100755 index 57f7d0f..0000000 --- a/dump/mapping/bufr_adpsfc_prepbufr.py +++ /dev/null @@ -1,256 +0,0 @@ -#!/usr/bin/env python3 -import os -import sys -import bufr -import argparse -import copy -import numpy as np -import numpy.ma as ma -import math -import calendar -import time -from datetime import datetime -from pyioda.ioda.Engines.Bufr import Encoder as iodaEncoder -from bufr.encoders.netcdf import Encoder as netcdfEncoder -from wxflow import Logger - -# Initialize Logger -# Get log level from the environment variable, default to 'INFO it not set -log_level = os.getenv('LOG_LEVEL', 'INFO') -logger = Logger('bufr_adpsfc_prepbufr.py', level=log_level, colored_log=False) - - -def logging(comm, level, message): - """ - Logs a message to the console or log file, based on the specified logging level. - - This function ensures that logging is only performed by the root process (`rank 0`) - in a distributed computing environment. The function maps the logging level to - appropriate logger methods and defaults to the 'INFO' level if an invalid level is provided. - - Parameters: - comm: object - The communicator object, typically from a distributed computing framework - (e.g., MPI). It must have a `rank()` method to determine the process rank. - level: str - The logging level as a string. Supported levels are: - - 'DEBUG' - - 'INFO' - - 'WARNING' - - 'ERROR' - - 'CRITICAL' - If an invalid level is provided, a warning will be logged, and the level - will default to 'INFO'. - message: str - The message to be logged. - - Behavior: - - Logs messages only on the root process (`comm.rank() == 0`). - - Maps the provided logging level to a method of the logger object. - - Defaults to 'INFO' and logs a warning if an invalid logging level is given. - - Supports standard logging levels for granular control over log verbosity. - - Example: - >>> logging(comm, 'DEBUG', 'This is a debug message.') - >>> logging(comm, 'ERROR', 'An error occurred!') - - Notes: - - Ensure that a global `logger` object is configured before using this function. - - The `comm` object should conform to MPI-like conventions (e.g., `rank()` method). - """ - - if comm.rank() == 0: - # Define a dictionary to map levels to logger methods - log_methods = { - 'DEBUG': logger.debug, - 'INFO': logger.info, - 'WARNING': logger.warning, - 'ERROR': logger.error, - 'CRITICAL': logger.critical, - } - - # Get the appropriate logging method, default to 'INFO' - log_method = log_methods.get(level.upper(), logger.info) - - if log_method == logger.info and level.upper() not in log_methods: - # Log a warning if the level is invalid - logger.warning(f'log level = {level}: not a valid level --> set to INFO') - - # Call the logging method - log_method(message) - - -def _compute_datetime(cycleTimeSinceEpoch, dhr): - """ - Compute dateTime using the cycleTimeSinceEpoch and Cycle Time - minus Cycle Time - - Parameters: - cycleTimeSinceEpoch: Time of cycle in Epoch Time - dhr: Observation Time Minus Cycle Time - - Returns: - Masked array of dateTime values - """ - - int64_fill_value = np.int64(0) - - dateTime = np.zeros(dhr.shape, dtype=np.int64) - for i in range(len(dateTime)): - if ma.is_masked(dhr[i]): - continue - else: - dateTime[i] = np.int64(dhr[i]*3600) + cycleTimeSinceEpoch - - dateTime = ma.array(dateTime) - dateTime = ma.masked_values(dateTime, int64_fill_value) - - return dateTime - - -def _make_description(mapping_path, cycle_time, update=False): - description = bufr.encoders.Description(mapping_path) - - ReferenceTime = np.int64(calendar.timegm(time.strptime(str(int(cycle_time)), '%Y%m%d%H'))) - - if update: - # Define the variables to be added in a list of dictionaries - variables = [ - { - 'name': 'MetaData/sequenceNumber', - 'source': 'variables/sequenceNumber', - 'units': '1', - 'longName': 'Sequence Number (Obs Subtype)', - } - ] - - # Loop through each variable and add it to the description - for var in variables: - description.add_variable( - name=var['name'], - source=var['source'], - units=var['units'], - longName=var['longName'] - ) - - # description.add_global(name='datetimeReference', value=str(ReferenceTime)) - - return description - - -def _make_obs(comm, input_path, mapping_path, cycle_time): - """ - Create the ioda adpsfc prepbufr observations: - - reads values - - adds sequenceNum - - Parameters - ---------- - comm: object - The communicator object (e.g., MPI) - input_path: str - The input bufr file - mapping_path: str - The input bufr2ioda mapping file - cycle_time: str - The cycle in YYYYMMDDHH format - """ - - # Get container from mapping file first - logging(comm, 'INFO', 'Get container from bufr') - container = bufr.Parser(input_path, mapping_path).parse(comm) - - logging(comm, 'DEBUG', f'container list (original): {container.list()}') - logging(comm, 'DEBUG', f'Change longitude range from [0,360] to [-180,180]') - lon = container.get('variables/longitude') - lon_paths = container.get_paths('variables/longitude') - lon[lon > 180] -= 360 - lon = ma.round(lon, decimals=2) - logging(comm, 'DEBUG', f'longitude max and min are {lon.max()}, {lon.min()}') - - logging(comm, 'DEBUG', f'Do DateTime calculation') - otmct = container.get('variables/obsTimeMinusCycleTime') - otmct_paths = container.get_paths('variables/obsTimeMinusCycleTime') - otmct2 = np.array(otmct) - cycleTimeSinceEpoch = np.int64(calendar.timegm(time.strptime(str(int(cycle_time)), '%Y%m%d%H'))) - dateTime = _compute_datetime(cycleTimeSinceEpoch, otmct2) - min_dateTime_ge_zero = min(x for x in dateTime if x > -1) - logging(comm, 'DEBUG', f'dateTime min/max = {min_dateTime_ge_zero} {dateTime.max()}') - - logging(comm, 'DEBUG', f'Make an array of 0s for MetaData/sequenceNumber') - sequenceNum = np.zeros(lon.shape, dtype=np.int32) - logging(comm, 'DEBUG', f' sequenceNummin/max = {sequenceNum.min()} {sequenceNum.max()}') - - logging(comm, 'DEBUG', f'Update variables in container') - container.replace('variables/longitude', lon) - container.replace('variables/timestamp', dateTime) - - logging(comm, 'DEBUG', f'Add variables to container') - container.add('variables/sequenceNumber', sequenceNum, lon_paths) - - # Check - logging(comm, 'DEBUG', f'container list (updated): {container.list()}') - - return container - - -def create_obs_group(input_path, mapping_path, cycle_time, env): - - comm = bufr.mpi.Comm(env["comm_name"]) - - logging(comm, 'INFO', f'Make description and make obs') - - container = _make_obs(comm, input_path, mapping_path, cycle_time) - description = _make_description(mapping_path, cycle_time, update=True) - - # Gather data from all tasks into all tasks. Each task will have the complete record - logging(comm, 'INFO', f'Gather data from all tasks into all tasks') - container.all_gather(comm) - - logging(comm, 'INFO', f'Encode the data') - data = next(iter(iodaEncoder(description).encode(container).values())) - - logging(comm, 'INFO', f'Return the encoded data.') - - return data - - -def create_obs_file(input_path, mapping_path, output_path, cycle_time): - - comm = bufr.mpi.Comm("world") - container = _make_obs(comm, input_path, mapping_path, cycle_time) - container.gather(comm) - - description = _make_description(mapping_path, cycle_time, update=True) - - # Encode the data - if comm.rank() == 0: - netcdfEncoder(description).encode(container, output_path) - - logging(comm, 'INFO', f'Return the encoded data') - - -if __name__ == '__main__': - start_time = time.time() - - bufr.mpi.App(sys.argv) - comm = bufr.mpi.Comm("world") - - # Required input arguments as positional arguments - parser = argparse.ArgumentParser(description="Convert BUFR to NetCDF using a mapping file.") - parser.add_argument('input', type=str, help='Input BUFR file') - parser.add_argument('mapping', type=str, help='BUFR2IODA Mapping File') - parser.add_argument('output', type=str, help='Output NetCDF file') - parser.add_argument('cycle_time', type=str, help='cycle time in YYYYMMDDHH format') - - args = parser.parse_args() - infile = args.input - mapping = args.mapping - output = args.output - cycle_time = args.cycle_time - - create_obs_file(infile, mapping, output, cycle_time) - - end_time = time.time() - running_time = end_time - start_time - logging(comm, 'INFO', f'Total running time: {running_time}') diff --git a/dump/mapping/bufr_adpupa_prepbufr.py b/dump/mapping/bufr_adpupa_prepbufr.py deleted file mode 100644 index c859732..0000000 --- a/dump/mapping/bufr_adpupa_prepbufr.py +++ /dev/null @@ -1,296 +0,0 @@ -#!/usr/bin/env python3 -import os -import sys -import bufr -import argparse -import copy -import numpy as np -import numpy.ma as ma -import math -import calendar -import time -from datetime import datetime -from pyioda.ioda.Engines.Bufr import Encoder as iodaEncoder -from bufr.encoders.netcdf import Encoder as netcdfEncoder -from wxflow import Logger - -# Initialize Logger -# Get log level from the environment variable, default to 'INFO it not set -log_level = os.getenv('LOG_LEVEL', 'INFO') -logger = Logger('bufr_adpupa_prepbufr.py', level=log_level, colored_log=False) - - -def logging(comm, level, message): - """ - Logs a message to the console or log file, based on the specified logging level. - - This function ensures that logging is only performed by the root process (`rank 0`) - in a distributed computing environment. The function maps the logging level to - appropriate logger methods and defaults to the 'INFO' level if an invalid level is provided. - - Parameters: - comm: object - The communicator object, typically from a distributed computing framework - (e.g., MPI). It must have a `rank()` method to determine the process rank. - level: str - The logging level as a string. Supported levels are: - - 'DEBUG' - - 'INFO' - - 'WARNING' - - 'ERROR' - - 'CRITICAL' - If an invalid level is provided, a warning will be logged, and the level - will default to 'INFO'. - message: str - The message to be logged. - - Behavior: - - Logs messages only on the root process (`comm.rank() == 0`). - - Maps the provided logging level to a method of the logger object. - - Defaults to 'INFO' and logs a warning if an invalid logging level is given. - - Supports standard logging levels for granular control over log verbosity. - - Example: - >>> logging(comm, 'DEBUG', 'This is a debug message.') - >>> logging(comm, 'ERROR', 'An error occurred!') - - Notes: - - Ensure that a global `logger` object is configured before using this function. - - The `comm` object should conform to MPI-like conventions (e.g., `rank()` method). - """ - - if comm.rank() == 0: - # Define a dictionary to map levels to logger methods - log_methods = { - 'DEBUG': logger.debug, - 'INFO': logger.info, - 'WARNING': logger.warning, - 'ERROR': logger.error, - 'CRITICAL': logger.critical, - } - - # Get the appropriate logging method, default to 'INFO' - log_method = log_methods.get(level.upper(), logger.info) - - if log_method == logger.info and level.upper() not in log_methods: - # Log a warning if the level is invalid - logger.warning(f'log level = {level}: not a valid level --> set to INFO') - - # Call the logging method - log_method(message) - - -def _compute_datetime(cycleTimeSinceEpoch, hrdr): - """ - Compute dateTime using the cycleTimeSinceEpoch and Observation Time Minus Cycle Time - - Parameters: - cycleTimeSinceEpoch: Time of cycle in Epoch Time - hrdr: Observation Time Minus Cycle Time - - Returns: - Masked array of dateTime values - """ - - int64_fill_value = np.int64(0) - dateTime = np.zeros(hrdr.shape, dtype=np.int64) - for i in range(len(dateTime)): - if ma.is_masked(hrdr[i]): - continue - else: - dateTime[i] = np.int64(hrdr[i]*3600) + cycleTimeSinceEpoch - - dateTime = ma.array(dateTime) - dateTime = ma.masked_values(dateTime, int64_fill_value) - - return dateTime - - -def _make_description(mapping_path, cycle_time, update=False): - description = bufr.encoders.Description(mapping_path) - - reference_time = np.int64(calendar.timegm(time.strptime(str(int(cycle_time)), '%Y%m%d%H'))) - - if update: - # Define the variables to be added in a list of dictionaries - variables = [ - { - 'name': 'MetaData/sequenceNumber', - 'source': 'variables/sequenceNumber', - 'units': '1', - 'longName': 'Sequence Number (Obs Subtype)', - }, - ] - - # Loop through each variable and add it to the description - for var in variables: - description.add_variable( - name=var['name'], - source=var['source'], - units=var['units'], - longName=var['longName'] - ) - - description.add_global(name='datetimeReference', value=str(reference_time)) - - return description - - -def _make_obs(comm, input_path, mapping_path, cycle_time): - """ - Create the ioda adpupa prepbufr observations: - - reads values - - adds sequenceNum - - Parameters - ---------- - comm: object - The communicator object (e.g., MPI) - input_path: str - The input bufr file - mapping_path: str - The input bufr2ioda mapping file - cycle_time: str - The cycle in YYYYMMDDHH format - """ - - # Get container from mapping file first - logging(comm, 'INFO', 'Get container from bufr') - container = bufr.Parser(input_path, mapping_path).parse(comm) - - logging(comm, 'DEBUG', f'container list (original): {container.list()}') - logging(comm, 'DEBUG', f'prepbufrDataLevelCategory') - cat = container.get('variables/prepbufrDataLevelCategory') - - logging(comm, 'DEBUG', f'Do DateTime calculation') - hrdr = container.get('variables/obsTimeMinusCycleTime') - hrdr2 = np.array(hrdr) - cycleTimeSinceEpoch = np.int64(calendar.timegm(time.strptime(str(int(cycle_time)), '%Y%m%d%H'))) - dateTime = _compute_datetime(cycleTimeSinceEpoch, hrdr2) - logging(comm, 'DEBUG', f'dateTime min/max = {dateTime.min()} {dateTime.max()}') - - logging(comm, 'DEBUG', f'Change longitude range from [0,360] to [-180,180]') - lon = container.get('variables/longitude') - lon[lon>180] -= 360 - - logging(comm, 'DEBUG', f'Make an array of 0s for MetaData/sequenceNumber') - sequenceNum = np.zeros(lon.shape, dtype=np.int32) - logging(comm, 'DEBUG', f' sequenceNummin/max = {sequenceNum.min()} {sequenceNum.max()}') - - logging(comm, 'DEBUG', f'Do ps calculation') - pob = container.get('variables/pressure') - ps = np.full(pob.shape[0], pob.fill_value) - ps = np.where(cat == 0, pob, ps) - - logging(comm, 'DEBUG', f'Do tsen and tv calculation') - tpc = container.get('variables/temperatureEventProgramCode') - tob = container.get('variables/airTemperature') - tsen = np.full(tob.shape[0], tob.fill_value) - tsen = np.where(((tpc >= 1) & (tpc < 8)), tob, tsen) - tvo = np.full(tob.shape[0], tob.fill_value) - tvo = np.where((tpc == 8), tob, tvo) - - logging(comm, 'DEBUG', f'Do ps QM calculations') - pqm = container.get('variables/pressureQualityMarker') - psqm = np.full(pqm.shape[0], pqm.fill_value) - psqm = np.where(cat == 0, pqm, psqm) - - logging(comm, 'DEBUG', f'Do tsen and tv QM calculations') - tobqm = container.get('variables/airTemperatureQualityMarker') - tsenqm = np.full(tobqm.shape[0], tobqm.fill_value) - tsenqm = np.where(((tpc >= 1) & (tpc < 8)), tobqm, tsenqm) - tvoqm = np.full(tobqm.shape[0], tobqm.fill_value) - tvoqm = np.where((tpc == 8), tobqm, tvoqm) - - logging(comm, 'DEBUG', f'Do ps ObsError calculations') - poe = container.get('variables/pressureError') - psoe = np.full(poe.shape[0], poe.fill_value) - psoe = np.where(cat == 0, poe, psoe) - - logging(comm, 'DEBUG', f'Do tsen and tv ObsError calculations') - toboe = container.get('variables/airTemperatureError') - tsenoe = np.full(toboe.shape[0], toboe.fill_value) - tsenoe = np.where(((tpc >= 1) & (tpc < 8)), toboe, tsenoe) - tvooe = np.full(toboe.shape[0], toboe.fill_value) - tvooe = np.where((tpc == 8), toboe, tvooe) - - logging(comm, 'DEBUG', f'Update variables in container') - container.replace('variables/longitude', lon) - container.replace('variables/timestamp', dateTime) - container.replace('variables/airTemperature', tsen) - container.replace('variables/airTemperatureQualityMarker', tsenqm) - container.replace('variables/airTemperatureError', tsenoe) - container.replace('variables/virtualTemperature', tvo) - container.replace('variables/virtualTemperatureQualityMarker', tvoqm) - container.replace('variables/virtualTemperatureError', tvooe) - container.replace('variables/stationPressure', ps) - container.replace('variables/stationPressureQualityMarker', psqm) - container.replace('variables/stationPressureError', psoe) - - logging(comm, 'DEBUG', f'Add variables to container') - container.add('variables/sequenceNumber', sequenceNum, lon_paths) - - # Check - logging(comm, 'DEBUG', f'container list (updated): {container.list()}') - - return container - - -def create_obs_group(input_path, mapping_path, cycle_time, env): - - comm = bufr.mpi.Comm(env["comm_name"]) - - logging(comm, 'INFO', f'Make description and make obs') - description = _make_description(mapping_path, cycle_time, update=True) - container = _make_obs(comm, input_path, mapping_path, cycle_time) - - # Gather data from all tasks into all tasks. Each task will have the complete record - logging(comm, 'INFO', f'Gather data from all tasks into all tasks') - container.all_gather(comm) - - logging(comm, 'INFO', f'Encode the data') - data = next(iter(iodaEncoder(description).encode(container).values())) - - logging(comm, 'INFO', f'Return the encoded data.') - return data - - -def create_obs_file(input_path, mapping_path, output_path, cycle_time): - - comm = bufr.mpi.Comm("world") - container = _make_obs(comm, input_path, mapping_path, cycle_time) - container.gather(comm) - - description = _make_description(mapping_path, cycle_time, update=True) - - # Encode the data - if comm.rank() == 0: - netcdfEncoder(description).encode(container, output_path) - - logging(comm, 'INFO', f'Return the encoded data') - - -if __name__ == '__main__': - start_time = time.time() - - bufr.mpi.App(sys.argv) - comm = bufr.mpi.Comm("world") - - # Required input arguments as positional arguments - parser = argparse.ArgumentParser(description="Convert BUFR to NetCDF using a mapping file.") - parser.add_argument('input', type=str, help='Input BUFR file') - parser.add_argument('mapping', type=str, help='BUFR2IODA Mapping File') - parser.add_argument('output', type=str, help='Output NetCDF file') - parser.add_argument('cycle_time', type=str, help='cycle time in YYYYMMDDHH format') - - args = parser.parse_args() - infile = args.input - mapping = args.mapping - output = args.output - cycle_time = args.cycle_time - - create_obs_file(infile, mapping, output, cycle_time) - - end_time = time.time() - running_time = end_time - start_time - logging(comm, 'INFO', f'Total running time: {running_time}') diff --git a/dump/mapping/bufr_adpupa_prepbufr_mapping.yaml b/dump/mapping/bufr_adpupa_prepbufr_mapping.yaml deleted file mode 100644 index dced147..0000000 --- a/dump/mapping/bufr_adpupa_prepbufr_mapping.yaml +++ /dev/null @@ -1,371 +0,0 @@ -# (C) Copyright 2024 NOAA/NWS/NCEP/EMC - -bufr: - group_by_variable: prepbufrDataLevelCategory - subsets: - - ADPUPA - variables: - # ObsType - observationType: - query: "*/TYP" - - # MetaData - timestamp: - timeoffset: - timeOffset: "*/PRSLEVEL/DRFTINFO/HRDR" - transforms: - - scale: 3600 - referenceTime: "2021-08-01T00:00:00Z" - - prepbufrDataLevelCategory: - query: "*/PRSLEVEL/CAT" - - obsTimeMinusCycleTime: - query: "*/PRSLEVEL/DRFTINFO/HRDR" - - longitude: - query: "*/PRSLEVEL/DRFTINFO/XDR" - - latitude: - query: "*/PRSLEVEL/DRFTINFO/YDR" - - stationIdentification: - query: "*/SID" - - stationElevation: - query: "*/ELV" - - temperatureEventProgramCode: - query: "*/PRSLEVEL/T___INFO/T__EVENT{1}/TPC" - - pressure: - query: "*/PRSLEVEL/P___INFO/P__EVENT{1}/POB" - type: float - transforms: - - scale: 100 - - height: - query: "*/PRSLEVEL/Z___INFO/Z__EVENT{1}/ZOB" - type: float - - # ObsValue - stationPressure: - query: "*/PRSLEVEL/P___INFO/P__EVENT{1}/POB" - type: float - transforms: - - scale: 100 - - airTemperature: - query: "*/PRSLEVEL/T___INFO/T__EVENT{1}/TOB" - transforms: - - offset: 273.15 - - dewPointTemperature: - query: "*/PRSLEVEL/Q___INFO/TDO" - transforms: - - offset: 273.15 - - virtualTemperature: - query: "*/PRSLEVEL/T___INFO/TVO" - transforms: - - offset: 273.15 - - windEastward: - query: "*/PRSLEVEL/W___INFO/W__EVENT{1}/UOB" - - windNorthward: - query: "*/PRSLEVEL/W___INFO/W__EVENT{1}/VOB" - - specificHumidity: - query: "*/PRSLEVEL/Q___INFO/Q__EVENT{1}/QOB" - type: float - transforms: - - scale: 0.000001 - - # QualityMark - pressureQualityMarker: - query: "*/PRSLEVEL/P___INFO/P__EVENT{1}/PQM" - - stationPressureQualityMarker: - query: "*/PRSLEVEL/P___INFO/P__EVENT{1}/PQM" - - heightQualityMarker: - query: "*/PRSLEVEL/Z___INFO/Z__EVENT{1}/ZQM" - - airTemperatureQualityMarker: - query: "*/PRSLEVEL/T___INFO/T__EVENT{1}/TQM" - - virtualTemperatureQualityMarker: - query: "*/PRSLEVEL/T___INFO/T__EVENT{1}/TQM" - - windQualityMarker: - query: "*/PRSLEVEL/W___INFO/W__EVENT{1}/WQM" - - specificHumidityQualityMarker: - query: "*/PRSLEVEL/Q___INFO/Q__EVENT{1}/QQM" - - # ObsError - pressureError: - query: "*/PRSLEVEL/P___INFO/P__BACKG/POE" - transforms: - - scale: 100 - - stationPressureError: - query: "*/PRSLEVEL/P___INFO/P__BACKG/POE" - transforms: - - scale: 100 - - airTemperatureError: - query: "*/PRSLEVEL/T___INFO/T__BACKG/TOE" - transforms: - - offset: 273.15 - - virtualTemperatureError: - query: "*/PRSLEVEL/T___INFO/T__BACKG/TOE" - transforms: - - offset: 273.15 - - windError: - query: "*/PRSLEVEL/W___INFO/W__BACKG/WOE" - - specificHumidityError: - query: "*/PRSLEVEL/Q___INFO/Q__BACKG/QOE" - - -encoder: - type: netcdf - - globals: - - name: "dataOriginalFormatSpec" - type: string - value: "prepbufr" - - - name: "platforms" - type: string - value: "ADPUPA" - - - name: "source" - type: string - value: "prepBUFR" - - - name: "description" - type: string - value: "UPPER-AIR (RAOB, PIBAL, RECCO, DROPS) REPORTS" - - - name: "dataProviderOrigin" - type: string - value: "U.S. NOAA NCEP" - - variables: - # Observation Type: stationPressure - - name: "ObsType/stationPressure" - source: variables/observationType - longName: "Observation Type" - - # Observation Type: airTemperature - - name: "ObsType/airTemperature" - source: variables/observationType - longName: "Observation Type" - - # Observation Type: virtualTemperature - - name: "ObsType/virtualTemperature" - source: variables/observationType - longName: "Observation Type" - - # Observation Type: windEastward - - name: "ObsType/windEastward" - source: variables/observationType - longName: "Observation Type" - - # Observation Type: windNorthward - - name: "ObsType/windNorthward" - source: variables/observationType - longName: "Observation Type" - - # Observation Type: specificHumidity - - name: "ObsType/specificHumidity" - source: variables/observationType - longName: "Observation Type" - - # MetaData: dateTime - - name: "MetaData/dateTime" - source: variables/timestamp - longName: "Datetime" - units: "seconds since 1970-01-01T00:00:00Z" - - # MetaData: stationIdentification - - name: "MetaData/stationIdentification" - source: variables/stationIdentification - longName: "Station ID" - - # MetaData: longitude - - name: "MetaData/longitude" - source: variables/longitude - longName: "Longitude" - units: "degrees_east" - range: [0, 360] - - # MetaData: latitude - - name: "MetaData/latitude" - source: variables/latitude - longName: "Latitude" - units: "degrees_north" - range: [-90, 90] - - # MetaData: stationElevation - - name: "MetaData/stationElevation" - source: variables/stationElevation - longName: "Height of Station" - units: "m" - - # MetaData: temperatureEventProgramCode - - name: "MetaData/temperatureEventProgramCode" - source: variables/temperatureEventProgramCode - longName: "Temperature Event Code" - - # MetaData: prepbufrDataLevelCategory - - name: "MetaData/prepbufrDataLevelCategory" - source: variables/prepbufrDataLevelCategory - longName: "Prepbufr Data Level Category" - - # MetaData: pressure - - name: "MetaData/pressure" - source: variables/pressure - longName: "Pressure" - units: "Pa" - - # MetaData: height - - name: "MetaData/height" - source: variables/height - longName: "Height of Observation" - units: "m" - - # QCFlags: qualityFlags - - name: "QCFlags/qualityFlags" - #- name: "MetaData/temperatureEventProgramCode" - source: variables/temperatureEventProgramCode - longName: "Temperature Event Program Code" - - # ObsValue: stationPressure - - name: "ObsValue/stationPressure" - source: variables/stationPressure - longName: "Station Pressure" - units: "Pa" - - # ObsValue: airTemperature - - name: "ObsValue/airTemperature" - source: variables/airTemperature - longName: "Temperature" - units: "K" - - # ObsValue: dewPointTemperature - - name: "ObsValue/dewPointTemperature" - source: variables/dewPointTemperature - longName: "DewPoint Temperature" - units: "K" - - # ObsValue: virtualTemperature - - name: "ObsValue/virtualTemperature" - source: variables/virtualTemperature - longName: "Virtual Temperature Non-Q Controlled" - units: "K" - - # ObsValue: windEastward - - name: "ObsValue/windEastward" - source: variables/windEastward - longName: "Eastward Wind" - units: "m s-1" - - # ObsValue: windNorthward - - name: "ObsValue/windNorthward" - source: variables/windNorthward - longName: "Northward Wind" - units: "m s-1" - - # ObsValue: specificHumidity - - name: "ObsValue/specificHumidity" - source: variables/specificHumidity - longName: "Specific Humidity" - units: "kg kg-1" - - # QualityMarker: pressure - - name: "QualityMarker/pressure" - source: variables/pressureQualityMarker - longName: "Pressure Quality Marker" - - # QualityMarker: height - - name: "QualityMarker/height" - source: variables/heightQualityMarker - longName: "Height Quality Marker" - - # QualityMarker: stationPressure - - name: "QualityMarker/stationPressure" - source: variables/stationPressureQualityMarker - longName: "Station Pressure Quality Marker" - - # QualityMarker: airTemperature - - name: "QualityMarker/airTemperature" - source: variables/airTemperatureQualityMarker - longName: "Temperature Quality Marker" - - # QualityMarker: virtualTemperature - - name: "QualityMarker/virtualTemperature" - source: variables/virtualTemperatureQualityMarker - longName: "Virtual Temperature Quality Marker" - - # QualityMarker: windEastward - - name: "QualityMarker/windEastward" - source: variables/windQualityMarker - longName: "Eastward Wind Quality Marker" - - # QualityMarker: windNorthward - - name: "QualityMarker/windNorthward" - source: variables/windQualityMarker - longName: "Northward Wind Quality Marker" - - # QualityMarker: specificHumidity - - name: "QualityMarker/specificHumidity" - source: variables/specificHumidityQualityMarker - longName: "Specific Humidity Quality Marker" - - # ObsError: pressure - - name: "ObsError/pressure" - source: variables/pressureError - longName: "Pressure Error" - units: "Pa" - - # ObsError: stationPressure - - name: "ObsError/stationPressure" - source: variables/stationPressureError - longName: "Station Pressure Error" - units: "Pa" - - # ObsError: airTemperature - - name: "ObsError/airTemperature" - source: variables/airTemperatureError - longName: "Temperature Error" - units: "K" - - # ObsError: virtualTemperature - - name: "ObsError/virtualTemperature" - source: variables/virtualTemperatureError - longName: "Virtual Temperature Error" - units: "K" - - # ObsError: windEastward - - name: "ObsError/windEastward" - source: variables/windError - longName: "East Wind Error" - units: "m s-1" - - # ObsError: windNorthward - - name: "ObsError/windNorthward" - source: variables/windError - longName: "North Wind Error" - units: "m s-1" - - # ObsError: specificHumidity - - name: "ObsError/specificHumidity" - source: variables/specificHumidityError - longName: "Specific Humidity Error" - units: "kg kg-1" diff --git a/dump/mapping/bufr_ahicsr.py b/dump/mapping/bufr_ahicsr.py new file mode 100644 index 0000000..a240278 --- /dev/null +++ b/dump/mapping/bufr_ahicsr.py @@ -0,0 +1,76 @@ +#!/usr/bin/env python3 +import os +import numpy as np + +import bufr +from bufr.obs_builder import ObsBuilder, add_main_functions, map_path + +MAPPING_PATH = map_path('bufr_ahicsr.yaml') + + +class BufrAhicsrObsBuilder(ObsBuilder): + + def __init__(self): + super().__init__(MAPPING_PATH, log_name=os.path.basename(__file__)) + + def compute_sensor_channel_number(self, nlocs, nchannels): + """ + Creates a 2D array of sensor channel numbers with shape (nlocs, nchannels), + where each row is [1, 2, ..., nchannels]. + + Args: + nlocs (int): Number of observation locations. + nchannels (int): Number of channels. + + Returns: + np.ndarray: 2D array with shape (nlocs, nchannels). + """ + return np.tile(np.arange(7, nchannels + 1, dtype=np.int32), (nlocs, 1)) + + def make_obs(self, comm, input_path): + # Get container from mapping file + self.log.info('Get container from bufr') + container = super().make_obs(comm, input_path) + + self.log.debug(f'Container list (original): {container.list()}') + self.log.debug(f'All_sub_categories = {container.all_sub_categories()}') + self.log.debug(f'Category map = {container.get_category_map()}') + + # Add a new derived variable, sensorChannelNumber into container + for cat in container.all_sub_categories(): + self.log.debug(f'category: {cat}') + + satId = container.get('satelliteId', cat) + sccf_paths = container.get_paths('sensorCentralFrequency', cat) + nlocs = satId.shape[0] # Number of locations + nchannels = 16 # Number of channels + # Add Ten Channels from 7 to 16 + sensor_channel_number = self.compute_sensor_channel_number(nlocs, nchannels) + self.log.debug(f'Adding derived variable: sensorChannelNumber for {nlocs} locations') + container.add('sensorChannelNumber', sensor_channel_number, sccf_paths, cat) + + if not np.any(satId): + self.log.warning(f'Category {cat[0]} does not exist in input file') + continue # Skip invalid category + + # Final container state + self.log.debug(f'Container list (updated): {container.list()}') + + return container + + def _make_description(self): + description = super()._make_description() + + description.add_variables([ + { + 'name': 'MetaData/sensorChannelNumber', + 'source': 'sensorChannelNumber', + 'units': '', + 'longName': 'Sensor Channel Number', + }, + ]) + + return description + + +add_main_functions(BufrAhicsrObsBuilder) diff --git a/dump/mapping/bufr_ahicsr.yaml b/dump/mapping/bufr_ahicsr.yaml new file mode 100755 index 0000000..d5f4583 --- /dev/null +++ b/dump/mapping/bufr_ahicsr.yaml @@ -0,0 +1,192 @@ +# (C) Copyright 2025 NOAA/NWS/NCEP/EMC +# +# This software is licensed under the terms of the Apache License Version 2.0 +# http://www.apache.org/licenses/LICENSE-2.0. + +bufr: + subsets: + - NC021044 + + variables: + # MetaData + timestamp: + datetime: + year: "*/YEAR" + month: "*/MNTH" + day: "*/DAYS" + hour: "*/HOUR" + minute: "*/MINU" + second: "*/SECO" + + dataReceiptTime: + datetime: + year: "*/RCYR" + month: "*/RCMO" + day: "*/RCDY" + hour: "*/RCHR" + minute: "*/RCMI" + + latitude: + query: "*/CLATH" + + longitude: + query: "*/CLONH" + + satelliteId: + query: "*/SAID" + + satelliteZenithAngle: + query: "*/SAZA" + + solarZenithAngle: + query: "*/SOZA" + + landOrSeaQualifier: + query: "*/LSQL" + + geopotentialHeight: + query: "*/HITE" + + sensorCentralFrequency: + query: "*/RPSEQ11{1-10}/SCCF" + + # ObsValue + cloudType: + query: "*/RPSEQ11{1-10}/CLTP" + + cloudAmount: + query: "*/RPSEQ11{1-10}/CLDMNT" + + cloudFree: + query: "*/RPSEQ11{1-10}/NCLDMNT" + + brightnessTemperature: + query: "*/RPSEQ11{1-10}/TMBRST" + + clearSkyStandardDeviation: + query: "*/RPSEQ11{1-10}/SDTB" + + splits: + satId: + category: + variable: satelliteId + map: + _173: h8 + _174: h9 + +encoder: + dimensions: + - name: Channel + path: "*/RPSEQ11{1-10}" + + globals: + - name: "platformCommonName" + type: string + value: "Himawari" + + - name: "platformLongDescription" + type: string + value: "JMA HIMAWARI-8 Clear Sky Radiances (CSR)" + + - name: "sensor" + type: string + value: "Advanced Himawari Imager (AHI)" + + - name: "providerFullName" + type: string + value: "Japan Meteorological Agency (JMA)" + + - name: "source" + type: string + value: "U.S. National Weather Service, National Centers for Environmental Prediction (NCEP)" + + - name: "sourceFiles" + type: string + value: "NCEP BUFR Dump containing NC021044 subset for JMA HIMAWARI-8 Clear Sky Radiances (CSR)" + + - name: "processingLevel" + type: string + value: "Level-2" + + - name: "converter" + type: string + value: "bufr-query" + + variables: + # MetaData + - name: "MetaData/dateTime" + source: variables/timestamp + longName: "Datetime" + units: "seconds since 1970-01-01T00:00:00Z" + + - name: "MetaData/dataReceiptTime" + source: variables/dataReceiptTime + longName: "Observation Receipt Time" + units: "seconds since 1970-01-01T00:00:00Z" + + - name: "MetaData/latitude" + source: variables/latitude + longName: "Latitude" + units: "degrees_north" + range: [ -90, 90 ] + + - name: "MetaData/longitude" + source: variables/longitude + longName: "Longitude" + units: "degrees_east" + range: [ -180, 180 ] + + - name: "MetaData/satelliteIdentifier" + source: variables/satelliteId + longName: "Satellite Identifier" + + - name: "MetaData/satelliteZenithAngle" + source: variables/satelliteZenithAngle + longName: "Satellite Zenith Angle" + units: "degree" + range: [ 0, 90 ] + + - name: "MetaData/solarZenithAngle" + source: variables/solarZenithAngle + longName: "Solar Zenith Angle" + units: "degree" + range: [0, 180] + + - name: "MetaData/landOrSeaQualifier" + source: variables/landOrSeaQualifier + longName: "Land/Sea Qualifier" + + - name: "MetaData/geopotentialHeight" + source: variables/geopotentialHeight + longName: "Geopotential Height" + units: "m" + + - name: "MetaData/sensorCentralFrequency" + source: variables/sensorCentralFrequency + longName: "Satellite Channel Center Frequency" + units: "Hz" + + # ObsValue + - name: "ObsValue/cloudType" + source: variables/cloudType + longName: "Cloud Type" + + - name: "ObsValue/cloudAmount" + source: variables/cloudAmount + longName: "Cloud Amount In Segment" + units: "1" + + - name: "ObsValue/cloudFree" + source: variables/cloudFree + longName: "Segment Amount Cloud Free" + units: "1" + + - name: "ObsValue/brightnessTemperature" + source: variables/brightnessTemperature + longName: "Brightness Temperature" + units: "K" + + - name: "ObsValue/clearSkyStandardDeviation" + source: variables/clearSkyStandardDeviation + longName: "Brightness Temperature Standard Deviation" + units: "K" diff --git a/dump/mapping/bufr_amsua.py b/dump/mapping/bufr_amsua.py index 6c88b95..12e746c 100644 --- a/dump/mapping/bufr_amsua.py +++ b/dump/mapping/bufr_amsua.py @@ -1,90 +1,13 @@ #!/usr/bin/env python3 -import argparse -import bufr -import netCDF4 as nc import os -import sys -import time -from pyioda.ioda.Engines.Bufr import Encoder -from bufr.encoders.netcdf import Encoder as netcdfEncoder -from wxflow import Logger - -INFO = 'INFO' -DEBUG = 'DEBUG' - -YAML_NORMAL = False # current as normal need remap for path2/1bama -INVALID = 1000.0 - -# Cosmic background temperature. Taken from Mather,J.C. et. al., 1999, "Calibrator Design for the COBE -# Far-Infrared Absolute Spectrophotometer (FIRAS)"Astrophysical Journal, vol 512, pp 511-520 -COSMIC_BACKGROUND_TEMP = 2.7253 - -nc_dir = './aux' - -# Initialize Logger -# Get log level from the environment variable, default to 'INFO it not set -log_level = os.getenv('LOG_LEVEL', INFO) -logger = Logger(os.path.basename(__file__), level=log_level, colored_log=True) - - -def logging(comm, message, level=INFO): - """ - Logs a message to the console or log file, based on the specified logging level. - - This function ensures that logging is only performed by the root process (`rank 0`) - in a distributed computing environment. The function maps the logging level to - appropriate logger methods and defaults to the 'INFO' level if an invalid level is provided. - - Parameters: - comm: object - The communicator object, typically from a distributed computing framework - (e.g., MPI). It must have a `rank()` method to determine the process rank. - level: str - The logging level as a string. Supported levels are: - - 'DEBUG' - - 'INFO' - - 'WARNING' - - 'ERROR' - - 'CRITICAL' - If an invalid level is provided, a warning will be logged, and the level - will default to 'INFO'. - message: str - The message to be logged. - - Behavior: - - Logs messages only on the root process (`comm.rank() == 0`). - - Maps the provided logging level to a method of the logger object. - - Defaults to 'INFO' and logs a warning if an invalid logging level is given. - - Supports standard logging levels for granular control over log verbosity. - - Example: - >>> logging(comm, 'DEBUG', 'This is a debug message.') - >>> logging(comm, 'ERROR', 'An error occurred!') - - Notes: - - Ensure that a global `logger` object is configured before using this function. - - The `comm` object should conform to MPI-like conventions (e.g., `rank()` method). - """ - - if comm.rank() == 0: - # Define a dictionary to map levels to logger methods - log_methods = { - 'DEBUG': logger.debug, - 'INFO': logger.info, - 'WARNING': logger.warning, - 'ERROR': logger.error, - 'CRITICAL': logger.critical, - } +import numpy as np +import numpy.ma as ma - # Get the appropriate logging method, default to 'INFO' - log_method = log_methods.get(level.upper(), logger.info) - - if log_method == logger.info and level.upper() not in log_methods: - # Log a warning if the level is invalid - logger.warning(f'log level = {level}: not a valid level --> set to INFO') +import bufr +from bufr.obs_builder import ObsBuilder, add_main_functions - # Call the logging method - log_method(message) +AMUA_1B_MAPPING = map_path('bufr_amsua_1bamua.yaml') +AMUA_ES_MAPPING = map_path('bufr_amsua_esamua.yaml') class ACCoeff: @@ -100,80 +23,68 @@ def __init__(self, ac_dir, sat_id='n19'): self.a_sp = self.a_space * COSMIC_BACKGROUND_TEMP -def _remove_ant_corr(i, ac, ifov, t): - # AC: Structure containing the antenna correction coefficients for the sensor of interest. - # iFOV: The FOV index for a scanline of the sensor of interest. - # T: On input, this argument contains the brightness - - t = ac.a_ep[i, ifov] * t + ac.a_sp[i, ifov] - t[(ifov < 1) | (ifov > ac.n_fovs)] = [INVALID] - return t - - -def _apply_ant_corr(i, ac, ifov, t): - # t: on input, this argument contains the antenna temperatures for the sensor channels. - t = (t - ac.a_sp[i, ifov]) / ac.a_ep[i, ifov] - t[(ifov < 1) | (ifov > ac.n_fovs)] = [INVALID] - return t - - -def _make_description(yaml_path): - description = bufr.encoders.Description(yaml_path) - return description - - -def _apply_corr(comm, sat_id, ta, ifov): - ac = ACCoeff(nc_dir, sat_id=sat_id) - if sat_id not in ['n15', 'n16']: - # Convert antenna temperature to brightness temperature - ifov = ifov.astype(int) - 1 - for i in range(ta.shape[1]): - logging(comm, f'inside loop for allpy ta to tb: i = {i}', level=DEBUG) - x = ta[:, i] - if YAML_NORMAL: - x = _apply_ant_corr(i, ac, ifov, x) - else: - x = _remove_ant_corr(i, ac, ifov, x) - x[x >= INVALID] = INVALID - ta[:, i] = x - return ta - - -def _re_map_variable(comm, container): - # read_bufrtovs.f90 - # antcorr_application.f90 - # search the keyword “ta2tb” for details - sat_ids = container.all_sub_categories() - for sat_id in sat_ids: - logging(comm, f'Converting for {sat_id[0]}, ...') - ta = container.get('variables/brightnessTemperature', sat_id) - if ta.shape[0]: - ifov = container.get('variables/fieldOfViewNumber', sat_id) - tb = _apply_corr(comm, sat_id[0], ta, ifov) - container.replace('variables/brightnessTemperature', tb, sat_id) - - -def _make_obs(comm, input_path, yaml_path): - cache = bufr.DataCache.has(input_path, yaml_path) - if cache: - logging(comm, f'The cache existed get data container from it') - container = bufr.DataCache.get(input_path, yaml_path) - else: - # If cacache does not exist, get data into cache - # Get data info container first - logging(comm, f'The cache is not existed') - container = bufr.Parser(input_path, yaml_path).parse() - return cache, container - - -def _mark_one_data(comm, cache, input_path, yaml_path, category, container=None): - if cache: - logging(comm, f'The cache existed get data container from it') - bufr.DataCache.mark_finished(input_path, yaml_path, [category]) - else: - logging(comm, f'add original container list into a cache = {container.list()}') - bufr.DataCache.add(input_path, yaml_path, container.all_sub_categories(), container) - bufr.DataCache.mark_finished(input_path, yaml_path, [category]) +class PressureObsBuilder(ObsBuilder): + def __init__(self): + map_dict = {'1bamua': AMUA_1B_MAPPING, + 'esamua': AMUA_ES_MAPPING} + + super().__init__(map_dict, log_name=os.path.basename(__file__)) + + def _remove_ant_corr(self, i, ac, ifov, t): + # AC: Structure containing the antenna correction coefficients for the sensor of interest. + # iFOV: The FOV index for a scanline of the sensor of interest. + # T: On input, this argument contains the brightness + + t = ac.a_ep[i, ifov] * t + ac.a_sp[i, ifov] + t[(ifov < 1) | (ifov > ac.n_fovs)] = [INVALID] + return t + + def _apply_ant_corr(self, i, ac, ifov, t): + # t: on input, this argument contains the antenna temperatures for the sensor channels. + t = (t - ac.a_sp[i, ifov]) / ac.a_ep[i, ifov] + t[(ifov < 1) | (ifov > ac.n_fovs)] = [INVALID] + return t + + def _apply_corr(self, sat_id, ta, ifov): + ac = ACCoeff(nc_dir, sat_id=sat_id) + if sat_id not in ['n15', 'n16']: + # Convert antenna temperature to brightness temperature + ifov = ifov.astype(int) - 1 + for i in range(ta.shape[1]): + self.log.debug(f'inside loop for allpy ta to tb: i = {i}') + x = ta[:, i] + if YAML_NORMAL: + x = _apply_ant_corr(i, ac, ifov, x) + else: + x = _remove_ant_corr(i, ac, ifov, x) + x[x >= INVALID] = INVALID + ta[:, i] = x + return ta + + def _re_map_variable(self, comm, container): + # read_bufrtovs.f90 + # antcorr_application.f90 + # search the keyword “ta2tb” for details + sat_ids = container.all_sub_categories() + for sat_id in sat_ids: + self.log.info(f'Converting for {sat_id[0]}, ...') + ta = container.get('variables/brightnessTemperature', sat_id) + if ta.shape[0]: + ifov = container.get('variables/fieldOfViewNumber', sat_id) + tb = _apply_corr(comm, sat_id[0], ta, ifov) + container.replace('variables/brightnessTemperature', tb, sat_id) + + def make_obs(self, comm, input_path): + cache = bufr.DataCache.has(input_path, yaml_path) + if cache: + logging(comm, f'The cache existed get data container from it') + container = bufr.DataCache.get(input_path, yaml_path) + else: + # If cacache does not exist, get data into cache + # Get data info container first + logging(comm, f'The cache is not existed') + container = bufr.Parser(input_path, yaml_path).parse() + return cache, container def create_obs_group(input_path1, input_path2, yaml_1b, yaml_es, category, env): diff --git a/dump/mapping/bufr_1bamua_mapping.yaml b/dump/mapping/bufr_amsua_1bamua.yaml similarity index 99% rename from dump/mapping/bufr_1bamua_mapping.yaml rename to dump/mapping/bufr_amsua_1bamua.yaml index 5ecc943..ed93594 100644 --- a/dump/mapping/bufr_1bamua_mapping.yaml +++ b/dump/mapping/bufr_amsua_1bamua.yaml @@ -16,6 +16,7 @@ bufr: fieldOfViewNumber: query: '*/FOVN' + heightOfStation: query: '*/HMSL' diff --git a/dump/mapping/bufr_esamua_mapping.yaml b/dump/mapping/bufr_amsua_esamua.yaml similarity index 100% rename from dump/mapping/bufr_esamua_mapping.yaml rename to dump/mapping/bufr_amsua_esamua.yaml diff --git a/dump/mapping/bufr_atms.py b/dump/mapping/bufr_atms.py index 747ac25..f74724b 100755 --- a/dump/mapping/bufr_atms.py +++ b/dump/mapping/bufr_atms.py @@ -1,187 +1,42 @@ #!/usr/bin/env python3 -import sys import os -import argparse -import time import numpy as np -import bufr -from pyioda.ioda.Engines.Bufr import Encoder as iodaEncoder -from bufr.encoders.netcdf import Encoder as netcdfEncoder -from wxflow import Logger - -# Initialize Logger -# Get log level from the environment variable, default to 'INFO it not set -log_level = os.getenv('LOG_LEVEL', 'INFO') -logger = Logger('BUFR_atms.py', level=log_level, colored_log=False) - - -def logging(comm, level, message): - """ - Logs a message to the console or log file, based on the specified logging level. - - This function ensures that logging is only performed by the root process (`rank 0`) - in a distributed computing environment. The function maps the logging level to - appropriate logger methods and defaults to the 'INFO' level if an invalid level is provided. - - Parameters: - comm: object - The communicator object, typically from a distributed computing framework - (e.g., MPI). It must have a `rank()` method to determine the process rank. - level: str - The logging level as a string. Supported levels are: - - 'DEBUG' - - 'INFO' - - 'WARNING' - - 'ERROR' - - 'CRITICAL' - If an invalid level is provided, a warning will be logged, and the level - will default to 'INFO'. - message: str - The message to be logged. - - Behavior: - - Logs messages only on the root process (`comm.rank() == 0`). - - Maps the provided logging level to a method of the logger object. - - Defaults to 'INFO' and logs a warning if an invalid logging level is given. - - Supports standard logging levels for granular control over log verbosity. - - Example: - >>> logging(comm, 'DEBUG', 'This is a debug message.') - >>> logging(comm, 'ERROR', 'An error occurred!') - - Notes: - - Ensure that a global `logger` object is configured before using this function. - - The `comm` object should conform to MPI-like conventions (e.g., `rank()` method). - """ - - if comm.rank() == 0: - # Define a dictionary to map levels to logger methods - log_methods = { - 'DEBUG': logger.debug, - 'INFO': logger.info, - 'WARNING': logger.warning, - 'ERROR': logger.error, - 'CRITICAL': logger.critical, - } - - # Get the appropriate logging method, default to 'INFO' - log_method = log_methods.get(level.upper(), logger.info) - - if log_method == logger.info and level.upper() not in log_methods: - # Log a warning if the level is invalid - logger.warning(f'log level = {level}: not a valid level --> set to INFO') - - # Call the logging method - log_method(message) - - -def _make_description(mapping_path, update=False): - - description = bufr.encoders.Description(mapping_path) - - return description - - -def _make_obs(comm, input_path, mapping_path): - - # Get container from mapping file first - logging(comm, 'INFO', 'Get container from bufr') - container = bufr.Parser(input_path, mapping_path).parse(comm) - - logging(comm, 'DEBUG', f'container list (original): {container.list()}') - logging(comm, 'DEBUG', f'all_sub_categories = {container.all_sub_categories()}') - logging(comm, 'DEBUG', f'category map = {container.get_category_map()}') - - # Add new/derived data into container - for cat in container.all_sub_categories(): - - logging(comm, 'DEBUG', f'category = {cat}') - - satid = container.get('variables/satelliteId', cat) - if satid.size == 0: - logging(comm, 'WARNING', f'category {cat[0]} does not exist in input file') - - # Check - logging(comm, 'DEBUG', f'container list (updated): {container.list()}') - logging(comm, 'DEBUG', f'all_sub_categories {container.all_sub_categories()}') - - return container - - -def create_obs_group(input_path, mapping_path, category, env): - - comm = bufr.mpi.Comm(env["comm_name"]) - - description = _make_description(mapping_path, update=False) - - # Check the cache for the data and return it if it exists - logging(comm, 'DEBUG', f'Check if bufr.DataCache exists? {bufr.DataCache.has(input_path, mapping_path)}') - if bufr.DataCache.has(input_path, mapping_path): - container = bufr.DataCache.get(input_path, mapping_path) - logging(comm, 'INFO', f'Encode {category} from cache') - data = iodaEncoder(description).encode(container)[(category,)] - logging(comm, 'INFO', f'Mark {category} as finished in the cache') - bufr.DataCache.mark_finished(input_path, mapping_path, [category]) - logging(comm, 'INFO', f'Return the encoded data for {category}') - return data - - container = _make_obs(comm, input_path, mapping_path) - - # Gather data from all tasks into all tasks. Each task will have the complete record - logging(comm, 'INFO', f'Gather data from all tasks into all tasks') - container.all_gather(comm) - - logging(comm, 'INFO', f'Add container to cache') - # Add the container to the cache - bufr.DataCache.add(input_path, mapping_path, container.all_sub_categories(), container) - - # Encode the data - logging(comm, 'INFO', f'Encode {category}') - data = iodaEncoder(description).encode(container)[(category,)] - - logging(comm, 'INFO', f'Mark {category} as finished in the cache') - # Mark the data as finished in the cache - bufr.DataCache.mark_finished(input_path, mapping_path, [category]) - - logging(comm, 'INFO', f'Return the encoded data for {category}') - return data +import numpy.ma as ma +import bufr +from bufr.obs_builder import ObsBuilder, add_main_functions, map_path -def create_obs_file(input_path, mapping_path, output_path): - comm = bufr.mpi.Comm("world") - container = _make_obs(comm, input_path, mapping_path) - container.gather(comm) +MAPPING_PATH = map_path('bufr_atms.yaml') - description = _make_description(mapping_path, update=False) - # Encode the data - if comm.rank() == 0: - netcdfEncoder(description).encode(container, output_path) +class BufrAtmsObsBuilder(ObsBuilder): + def __init__(self): + super().__init__(MAPPING_PATH, log_name=os.path.basename(__file__)) - logging(comm, 'INFO', f'Return the encoded data') + def make_obs(self, comm, input_path): + # Get container from mapping file first + self.log.info('Get container from bufr') + container = super().make_obs(comm, input_path) + self.log.debug(f'container list (original): {container.list()}') + self.log.debug(f'all_sub_categories = {container.all_sub_categories()}') + self.log.debug(f'category map = {container.get_category_map()}') -if __name__ == '__main__': + # Add new/derived data into container + for cat in container.all_sub_categories(): - start_time = time.time() + self.log.debug(f'category = {cat}') - bufr.mpi.App(sys.argv) - comm = bufr.mpi.Comm("world") + satId = container.get('satelliteId', cat) + if not np.any(satId): + self.log.warning(f'category {cat[0]} does not exist in input file') - # Required input arguments as positional arguments - parser = argparse.ArgumentParser(description="Convert BUFR to NetCDF using a mapping file.") - parser.add_argument('input', type=str, help='Input BUFR file') - parser.add_argument('mapping', type=str, help='BUFR2IODA Mapping File') - parser.add_argument('output', type=str, help='Output NetCDF file') + # Check + self.log.debug(f'container list (updated): {container.list()}') + self.log.debug(f'all_sub_categories {container.all_sub_categories()}') - args = parser.parse_args() - mapping = args.mapping - infile = args.input - output = args.output + return container - create_obs_file(infile, mapping, output) - end_time = time.time() - running_time = end_time - start_time - logging(comm, 'INFO', f'Total running time: {running_time}') +add_main_functions(BufrAtmsObsBuilder) diff --git a/dump/mapping/bufr_atms_mapping.yaml b/dump/mapping/bufr_atms.yaml similarity index 85% rename from dump/mapping/bufr_atms_mapping.yaml rename to dump/mapping/bufr_atms.yaml index 4570445..e51cb3a 100755 --- a/dump/mapping/bufr_atms_mapping.yaml +++ b/dump/mapping/bufr_atms.yaml @@ -19,14 +19,12 @@ bufr: satelliteId: query: "*/SAID" - satelliteInstrument: - query: "*/SIID" - fieldOfViewNumber: query: "*/FOVN" heightOfStation: query: "*/HMSL" + type: float solarZenithAngle: query: "*/SOZA" @@ -75,14 +73,16 @@ bufr: _226: n21 encoder: -# type: netcdf - dimensions: - name: Channel source: variables/sensorChannelNumber path: "*/ATMSCH" globals: + - name: "description" + type: string + value: "ATMS Level-1B data processed using FFT-based beam convolution to standardize all channels to a common 3°×3° spatial resolution" + - name: "platformCommonName" type: string value: "JPSS" @@ -91,17 +91,21 @@ encoder: type: string value: "Joint Polar Satellite System" - - name: "sensor" + - name: "sensorCommonName" type: string - value: "Advanced Baseline Imager (ATMS)" + value: "ATMS" + + - name: "sensorLongDescription" + type: string + value: "Advanced Baseline Imager (ATMS) - 22 channels including the 54 and 183 GHz bands; cross-track scanning" - name: "source" type: string - value: "NCEP BUFR Dump for ATMS brightness temperature data (atms normal feed)" + value: "U.S. National Weather Service, National Centres for Environmental Prediction (NCEP))" - - name: "providerFullName" + - name: "sourceFiles" type: string - value: "National Environmental Satellite, Data, and Information Service" + value: "NCEP BUFR Dump - atms NC021203 Normal Feed" - name: "processingLevel" type: string @@ -109,7 +113,7 @@ encoder: - name: "converter" type: string - value: "BUFR" + value: "bufr-query" variables: # MetaData @@ -133,10 +137,6 @@ encoder: source: variables/satelliteId longName: "Satellite Identifier" - - name: "MetaData/satelliteInstrument" - source: variables/satelliteInstrument - longName: "Satellite Instrument" - - name: "MetaData/sensorScanPosition" source: variables/fieldOfViewNumber longName: "Field of View Number" diff --git a/dump/mapping/bufr_cris-fsr.py b/dump/mapping/bufr_cris-fsr.py new file mode 100755 index 0000000..93cecc7 --- /dev/null +++ b/dump/mapping/bufr_cris-fsr.py @@ -0,0 +1,41 @@ +#!/usr/bin/env python3 +import os +import numpy as np +import numpy.ma as ma + +import bufr +from bufr.obs_builder import ObsBuilder, add_main_functions, map_path + + +MAPPING_PATH = map_path('bufr_cris-fsr.yaml') + + +class BufrCrisObsBuilder(ObsBuilder): + def __init__(self): + super().__init__(MAPPING_PATH, log_name=os.path.basename(__file__)) + + def make_obs(self, comm, input_path): + # Get container from mapping file first + self.log.info('Get container from bufr') + container = super().make_obs(comm, input_path) + + self.log.debug(f'container list (original): {container.list()}') + self.log.debug(f'all_sub_categories = {container.all_sub_categories()}') + self.log.debug(f'category map = {container.get_category_map()}') + + # Add new/derived data into container + for cat in container.all_sub_categories(): + self.log.debug(f'category = {cat}') + + satId = container.get('satelliteId', cat) + if not np.any(satId): + self.log.warning(f'category {cat[0]} does not exist in input file') + + # Check + self.log.debug(f'container list (updated): {container.list()}') + self.log.debug(f'all_sub_categories {container.all_sub_categories()}') + + return container + + +add_main_functions(BufrCrisObsBuilder) diff --git a/dump/mapping/bufr_cris-fsr.yaml b/dump/mapping/bufr_cris-fsr.yaml new file mode 100755 index 0000000..a7824bd --- /dev/null +++ b/dump/mapping/bufr_cris-fsr.yaml @@ -0,0 +1,219 @@ +bufr: + variables: + # MetaData + timestamp: + datetime: + year: "*/YEAR" + month: "*/MNTH" + day: "*/DAYS" + hour: "*/HOUR" + minute: "*/MINU" + second: "*/SECO" + + latitude: + query: "*/CLATH" + + longitude: + query: "*/CLONH" + + satelliteId: + query: "*/SAID" + + sensorId: + query: "*/SIID[1]" + + scanLineNumber: + query: "*/SLNM" + + fieldOfViewNumber: + query: "*/FOVN" + + fieldOfRegardNumber: + query: "*/FORN" + + sensorScanPosition: + query: "*/FOVN" + + solarZenithAngle: + query: "*/SOZA" + + solarAzimuthAngle: + query: "*/SOLAZI" + + sensorZenithAngle: + query: "*/SAZA" + + sensorAzimuthAngle: + query: "*/BEARAZ" + + sensorViewAngle: + sensorScanAngle: + fieldOfViewNumber: "*/FOVN" + fieldOfRegardNumber: "*/FORN" + scanStart: -48.330 + scanStep: 3.3331 + sensor: cris + + sensorChannelNumber: + query: "*/CRCHNM/CHNM" + + # ObsValue + cloudCoverTotal: + query: "*/TOCC" + type: float + transforms: + - scale: 0.01 + + cloudTopHeight: + query: "*/HOCT" + type: float + + spectralRadiance: + query: "*/CRCHNM/SRAD" + transforms: + - scale: 0.01 + + splits: + satId: + category: + variable: satelliteId + map: + _224: npp + _225: n20 + _226: n21 + +encoder: + dimensions: + - name: Channel + source: variables/sensorChannelNumber + path: "*/CRCHNM" + + globals: + - name: "description" + type: string + value: "Level-1C CrIS radiance spectra apodized using the standard Hamming window, with full radiometric calibration and geolocation" + + - name: "platformCommonName" + type: string + value: "JPSS" + + - name: "platformLongDescription" + type: string + value: "NOAA Joint Polar Satellite System" + + - name: "sensorCommonName" + type: string + value: "CrIS" + + - name: "sensorLongDescription" + type: string + value: "Cross-track Infrared Sounder (Interferometer) with three IR bands with 2211 channels; cross-track scanning" + + - name: "source" + type: string + value: "U.S. National Weather Service, National Centres for Environmental Prediction (NCEP))" + + - name: "sourceFiles" + type: string + value: "NCEP BUFR Dump - crisf4 NC021206 Normal Feed 431 channel subset" + + - name: "processingLevel" + type: string + value: "Level-1C" + + - name: "converter" + type: string + value: "bufr-query" + + variables: + + # MetaData + - name: "MetaData/dateTime" + source: variables/timestamp + longName: "Datetime" + units: "seconds since 1970-01-01T00:00:00Z" + + - name: "MetaData/latitude" + source: variables/latitude + longName: "Latitude" + units: "degrees_north" + range: [ -90, 90 ] + + - name: "MetaData/longitude" + source: variables/longitude + longName: "Longitude" + units: "degrees_east" + range: [ -180, 180 ] + + - name: "MetaData/satelliteIdentifier" + source: variables/satelliteId + longName: "Satellite Identifier" + +# - name: "MetaData/instrumentIdentifier" +# source: variables/sensorId +# longName: "Satellite Instrument Identifier" + + - name: "MetaData/scanLineNumber" + source: variables/scanLineNumber + longName: "Scan Line Number" + + - name: "MetaData/fieldOfViewNumber" + source: variables/fieldOfViewNumber + longName: "Field of View Number" + + - name: "MetaData/fieldOfRegardNumber" + source: variables/fieldOfRegardNumber + longName: "Field of Regard Number" + + - name: "MetaData/sensorScanPosition" + source: variables/sensorScanPosition + longName: "Sensor Scan Position" + + - name: "MetaData/solarZenithAngle" + source: variables/solarZenithAngle + longName: "Solar Zenith Angle" + units: "degree" + range: [ 0, 180 ] + + - name: "MetaData/solarAzimuthAngle" + source: variables/solarAzimuthAngle + longName: "Solar Azimuth Angle" + units: "degree" + range: [ 0, 360 ] + + - name: "MetaData/sensorZenithAngle" + source: variables/sensorZenithAngle + longName: "Sensor Zenith Angle" + units: "degree" + range: [ 0, 90 ] + + - name: "MetaData/sensorAzimuthAngle" + source: variables/sensorAzimuthAngle + longName: "Sensor Azimuth Angle" + units: "degree" + range: [ 0, 360 ] + + - name: "MetaData/sensorViewAngle" + source: variables/sensorViewAngle + longName: "Sensor View Angle" + units: "degree" + + - name: "MetaData/sensorChannelNumber" + source: variables/sensorChannelNumber + longName: "Sensor Channel Number" + + - name: "ObsValue/cloudCoverTotal" + source: variables/cloudCoverTotal + longName: "Total Cloud Cover" + units: "1" + range: [ 0, 1 ] + + - name: "ObsValue/heightOfTopOfCloud" + source: variables/cloudTopHeight + longName: "Height of Top of Cloud" + units: "m" + + - name: "ObsValue/spectralRadiance" + source: variables/spectralRadiance + longName: "CrIS Spectral Radiance" + units: "W m-2 sr-1 m-1" diff --git a/dump/mapping/bufr_iasi.py b/dump/mapping/bufr_iasi.py new file mode 100755 index 0000000..e04b060 --- /dev/null +++ b/dump/mapping/bufr_iasi.py @@ -0,0 +1,42 @@ +#!/usr/bin/env python3 +import os +import numpy as np +import numpy.ma as ma + +import bufr +from bufr.obs_builder import ObsBuilder, add_main_functions, map_path + + +MAPPING_PATH = map_path('bufr_iasi.yaml') + + +class BufrIasiObsBuilder(ObsBuilder): + def __init__(self): + super().__init__(MAPPING_PATH, log_name=os.path.basename(__file__)) + + def make_obs(self, comm, input_path): + # Get container from mapping file first + self.log.info('Get container from bufr') + container = super().make_obs(comm, input_path) + + self.log.debug(f'container list (original): {container.list()}') + self.log.debug(f'all_sub_categories = {container.all_sub_categories()}') + self.log.debug(f'category map = {container.get_category_map()}') + + # Add new/derived data into container + for cat in container.all_sub_categories(): + + self.log.debug(f'category = {cat}') + + satId = container.get('satelliteId', cat) + if not np.any(satId): + self.log.warning(f'category {cat[0]} does not exist in input file') + + # Check + self.log.debug(f'container list (updated): {container.list()}') + self.log.debug(f'all_sub_categories {container.all_sub_categories()}') + + return container + + +add_main_functions(BufrIasiObsBuilder) diff --git a/dump/mapping/bufr_mtiasi_mapping.yaml b/dump/mapping/bufr_iasi.yaml similarity index 85% rename from dump/mapping/bufr_mtiasi_mapping.yaml rename to dump/mapping/bufr_iasi.yaml index f28e98f..bb78c3b 100755 --- a/dump/mapping/bufr_mtiasi_mapping.yaml +++ b/dump/mapping/bufr_iasi.yaml @@ -60,13 +60,16 @@ bufr: scanStepAdjust: 0.625 sensor: iasi - fractionOfClearPixelsInFOV: - query: "*/IASIL1CS{1}/FCPH" - sensorChannelNumber: query: "*/IASICHN/CHNM" # ObsValue + fractionOfClearPixelsInFOV: + query: "*/IASIL1CS{1}/FCPH" + type: float + transforms: + - scale: 0.01 + spectralRadiance: spectralRadiance: sensorChannelNumber: "*/IASICHN/CHNM" @@ -85,43 +88,39 @@ bufr: _5: metop-c encoder: - type: netcdf - dimensions: - name: Channel source: variables/sensorChannelNumber path: "*/IASICHN" - - name: Cluster - path: "*/IASIL1CS" - - - name: Band - path: "*/IASIL1CB" - globals: + - name: "description" + type: string + value: "IASI Level-1C radiance spectra, calibrated, geolocated, and apodized using the standard Hamming window to suppress spectral sidelobes" + - name: "platformCommonName" type: string - value: "Meteorological Operational Satellite" + value: "MetOp" - name: "platformLongDescription" type: string - value: "EUMETSAT Polar System in sunsynchronous orbit" + value: "Meteorological Operational Satellite" - - name: "source" + - name: "sensorCommonName" type: string - value: "MTYP 021-241 IASI 1C RADIANCES (VARIABLE CHNS) (METOP)" + value: "IASI" - - name: "sourceFiles" + - name: "sensorLongDescription" type: string - value: "gdas.t00z.mtiasi.tm00.bufr_d" + value: "Infrared Atmospheric Sounding Interferometer with 8461 channels with one embedded IR imaging channel; corss-track scanning" - - name: "datetimeReference" + - name: "source" type: string - value: "2021-08-01T00:00:00Z" + value: "U.S. National Weather Service, National Centres for Environmental Prediction (NCEP))" - - name: "sensor" + - name: "sourceFiles" type: string - value: "IASI" + value: "NCEP BUFR Dump - mtiasi NC021241 Normal Feed 616 channel subset" - name: "processingLevel" type: string @@ -129,7 +128,7 @@ encoder: - name: "converter" type: string - value: "BUFR" + value: "bufr-query" variables: @@ -213,17 +212,16 @@ encoder: source: variables/sensorChannelNumber longName: "Sensor Channel Number" - - name: "MetaData/fractionOfClearPixelsInFOV" + - name: "ObsValue/cloudFree" source: variables/fractionOfClearPixelsInFOV longName: "Fraction of Clear Pixels in a Field of View" units: "1" - range: [0, 100] + range: [ 0, 1 ] -# The unit from BUFR is W m-2 sr-1 m -- this is radiance per wavenumber - - name: "ObsValue/radiance" + - name: "ObsValue/spectralRadiance" source: variables/spectralRadiance longName: "IASI Spectral Radiance" - units: "W m-2 sr-1" # units: "W m-2 sr-1 m-1" + units: "W m-2 sr-1 m-1" # chunks: [1000, 154] # chunks: [1000, 616] # compressionLevel: 4 diff --git a/dump/mapping/bufr_mtiasi.py b/dump/mapping/bufr_mtiasi.py deleted file mode 100755 index 67d59a6..0000000 --- a/dump/mapping/bufr_mtiasi.py +++ /dev/null @@ -1,187 +0,0 @@ -#!/usr/bin/env python3 -import sys -import os -import argparse -import time -import numpy as np -import bufr -from pyioda.ioda.Engines.Bufr import Encoder as iodaEncoder -from bufr.encoders.netcdf import Encoder as netcdfEncoder -from wxflow import Logger - -# Initialize Logger -# Get log level from the environment variable, default to 'INFO it not set -log_level = os.getenv('LOG_LEVEL', 'INFO') -logger = Logger('BUFR2IODA_iasi.py', level=log_level, colored_log=False) - - -def logging(comm, level, message): - """ - Logs a message to the console or log file, based on the specified logging level. - - This function ensures that logging is only performed by the root process (`rank 0`) - in a distributed computing environment. The function maps the logging level to - appropriate logger methods and defaults to the 'INFO' level if an invalid level is provided. - - Parameters: - comm: object - The communicator object, typically from a distributed computing framework - (e.g., MPI). It must have a `rank()` method to determine the process rank. - level: str - The logging level as a string. Supported levels are: - - 'DEBUG' - - 'INFO' - - 'WARNING' - - 'ERROR' - - 'CRITICAL' - If an invalid level is provided, a warning will be logged, and the level - will default to 'INFO'. - message: str - The message to be logged. - - Behavior: - - Logs messages only on the root process (`comm.rank() == 0`). - - Maps the provided logging level to a method of the logger object. - - Defaults to 'INFO' and logs a warning if an invalid logging level is given. - - Supports standard logging levels for granular control over log verbosity. - - Example: - >>> logging(comm, 'DEBUG', 'This is a debug message.') - >>> logging(comm, 'ERROR', 'An error occurred!') - - Notes: - - Ensure that a global `logger` object is configured before using this function. - - The `comm` object should conform to MPI-like conventions (e.g., `rank()` method). - """ - - if comm.rank() == 0: - # Define a dictionary to map levels to logger methods - log_methods = { - 'DEBUG': logger.debug, - 'INFO': logger.info, - 'WARNING': logger.warning, - 'ERROR': logger.error, - 'CRITICAL': logger.critical, - } - - # Get the appropriate logging method, default to 'INFO' - log_method = log_methods.get(level.upper(), logger.info) - - if log_method == logger.info and level.upper() not in log_methods: - # Log a warning if the level is invalid - logger.warning(f'log level = {level}: not a valid level --> set to INFO') - - # Call the logging method - log_method(message) - - -def _make_description(mapping_path, update=False): - - description = bufr.encoders.Description(mapping_path) - - return description - - -def _make_obs(comm, input_path, mapping_path): - - # Get container from mapping file first - logging(comm, 'INFO', 'Get container from bufr') - container = bufr.Parser(input_path, mapping_path).parse(comm) - - logging(comm, 'DEBUG', f'container list (original): {container.list()}') - logging(comm, 'DEBUG', f'all_sub_categories = {container.all_sub_categories()}') - logging(comm, 'DEBUG', f'category map = {container.get_category_map()}') - - # Add new/derived data into container - for cat in container.all_sub_categories(): - - logging(comm, 'DEBUG', f'category = {cat}') - - satid = container.get('variables/satelliteId', cat) - if satid.size == 0: - logging(comm, 'WARNING', f'category {cat[0]} does not exist in input file') - - # Check - logging(comm, 'DEBUG', f'container list (updated): {container.list()}') - logging(comm, 'DEBUG', f'all_sub_categories {container.all_sub_categories()}') - - return container - - -def create_obs_group(input_path, mapping_path, category, env): - - comm = bufr.mpi.Comm(env["comm_name"]) - - description = _make_description(mapping_path, update=False) - - # Check the cache for the data and return it if it exists - logging(comm, 'DEBUG', f'Check if bufr.DataCache exists? {bufr.DataCache.has(input_path, mapping_path)}') - if bufr.DataCache.has(input_path, mapping_path): - container = bufr.DataCache.get(input_path, mapping_path) - logging(comm, 'INFO', f'Encode {category} from cache') - data = iodaEncoder(description).encode(container)[(category,)] - logging(comm, 'INFO', f'Mark {category} as finished in the cache') - bufr.DataCache.mark_finished(input_path, mapping_path, [category]) - logging(comm, 'INFO', f'Return the encoded data for {category}') - return data - - container = _make_obs(comm, input_path, mapping_path) - - # Gather data from all tasks into all tasks. Each task will have the complete record - logging(comm, 'INFO', f'Gather data from all tasks into all tasks') - container.all_gather(comm) - - logging(comm, 'INFO', f'Add container to cache') - # Add the container to the cache - bufr.DataCache.add(input_path, mapping_path, container.all_sub_categories(), container) - - # Encode the data - logging(comm, 'INFO', f'Encode {category}') - data = iodaEncoder(description).encode(container)[(category,)] - - logging(comm, 'INFO', f'Mark {category} as finished in the cache') - # Mark the data as finished in the cache - bufr.DataCache.mark_finished(input_path, mapping_path, [category]) - - logging(comm, 'INFO', f'Return the encoded data for {category}') - return data - - -def create_obs_file(input_path, mapping_path, output_path): - - comm = bufr.mpi.Comm("world") - container = _make_obs(comm, input_path, mapping_path) - container.gather(comm) - - description = _make_description(mapping_path, update=False) - - # Encode the data - if comm.rank() == 0: - netcdfEncoder(description).encode(container, output_path) - - logging(comm, 'INFO', f'Return the encoded data') - - -if __name__ == '__main__': - - start_time = time.time() - - bufr.mpi.App(sys.argv) - comm = bufr.mpi.Comm("world") - - # Required input arguments as positional arguments - parser = argparse.ArgumentParser(description="Convert BUFR to NetCDF using a mapping file.") - parser.add_argument('input', type=str, help='Input BUFR file') - parser.add_argument('mapping', type=str, help='BUFR2IODA Mapping File') - parser.add_argument('output', type=str, help='Output NetCDF file') - - args = parser.parse_args() - mapping = args.mapping - infile = args.input - output = args.output - - create_obs_file(infile, mapping, output) - - end_time = time.time() - running_time = end_time - start_time - logging(comm, 'INFO', f'Total running time: {running_time}') diff --git a/dump/mapping/bufr_satwnd_amv_abi.py b/dump/mapping/bufr_satwnd_amv_abi.py index 41ba77a..06fea6e 100755 --- a/dump/mapping/bufr_satwnd_amv_abi.py +++ b/dump/mapping/bufr_satwnd_amv_abi.py @@ -1,316 +1,37 @@ #!/usr/bin/env python3 -import sys + import os -import argparse -import time import numpy as np -import bufr -from pyioda.ioda.Engines.Bufr import Encoder as iodaEncoder -from bufr.encoders.netcdf import Encoder as netcdfEncoder -from wxflow import Logger - -# Initialize Logger -# Get log level from the environment variable, default to 'INFO it not set -log_level = os.getenv('LOG_LEVEL', 'INFO') -logger = Logger('BUFR2IODA_satwnd_amv_goes.py', level=log_level, colored_log=False) - - -def logging(comm, level, message): - """ - Logs a message to the console or log file, based on the specified logging level. - - This function ensures that logging is only performed by the root process (`rank 0`) - in a distributed computing environment. The function maps the logging level to - appropriate logger methods and defaults to the 'INFO' level if an invalid level is provided. - - Parameters: - comm: object - The communicator object, typically from a distributed computing framework - (e.g., MPI). It must have a `rank()` method to determine the process rank. - level: str - The logging level as a string. Supported levels are: - - 'DEBUG' - - 'INFO' - - 'WARNING' - - 'ERROR' - - 'CRITICAL' - If an invalid level is provided, a warning will be logged, and the level - will default to 'INFO'. - message: str - The message to be logged. - - Behavior: - - Logs messages only on the root process (`comm.rank() == 0`). - - Maps the provided logging level to a method of the logger object. - - Defaults to 'INFO' and logs a warning if an invalid logging level is given. - - Supports standard logging levels for granular control over log verbosity. - - Example: - >>> logging(comm, 'DEBUG', 'This is a debug message.') - >>> logging(comm, 'ERROR', 'An error occurred!') - - Notes: - - Ensure that a global `logger` object is configured before using this function. - - The `comm` object should conform to MPI-like conventions (e.g., `rank()` method). - """ - - if comm.rank() == 0: - # Define a dictionary to map levels to logger methods - log_methods = { - 'DEBUG': logger.debug, - 'INFO': logger.info, - 'WARNING': logger.warning, - 'ERROR': logger.error, - 'CRITICAL': logger.critical, - } - - # Get the appropriate logging method, default to 'INFO' - log_method = log_methods.get(level.upper(), logger.info) - - if log_method == logger.info and level.upper() not in log_methods: - # Log a warning if the level is invalid - logger.warning(f'log level = {level}: not a valid level --> set to INFO') - - # Call the logging method - log_method(message) - - -def _make_description(mapping_path, update=False): - description = bufr.encoders.Description(mapping_path) - - if update: - # Define the variables to be added in a list of dictionaries - variables = [ - { - 'name': 'ObsType/windEastward', - 'source': 'variables/obstype_uwind', - 'units': '1', - 'longName': 'Observation Type based on Satellite-derived Wind Computation Method and Spectral Band', - }, - { - 'name': 'ObsType/windNorthward', - 'source': 'variables/obstype_vwind', - 'units': '1', - 'longName': 'Observation Type based on Satellite-derived Wind Computation Method and Spectral Band', - }, - { - 'name': 'ObsValue/windEastward', - 'source': 'variables/windEastward', - 'units': 'm s-1', - 'longName': 'Eastward Wind Component', - }, - { - 'name': 'ObsValue/windNorthward', - 'source': 'variables/windNorthward', - 'units': 'm s-1', - 'longName': 'Northward Wind Component', - }, - ] - - # Loop through each variable and add it to the description - for var in variables: - description.add_variable( - name=var['name'], - source=var['source'], - units=var['units'], - longName=var['longName'] - ) - - return description - - -def compute_wind_components(wdir, wspd): - """ - Compute the U and V wind components from wind direction and wind speed. - - Parameters: - wdir (array-like): Wind direction in degrees (meteorological convention: 0° = North, 90° = East). - wspd (array-like): Wind speed (m/s). - - Returns: - tuple: U and V wind components as numpy arrays with dtype float32. - """ - wdir_rad = np.radians(wdir) # Convert degrees to radians - u = -wspd * np.sin(wdir_rad) - v = -wspd * np.cos(wdir_rad) - - return u.astype(np.float32), v.astype(np.float32) - - -def _get_obs_type(swcm, chanfreq): - """ - Determine the observation type based on `swcm` and `chanfreq`. - - Parameters: - swcm (array-like): Satellite derived wind calculation method. - chanfreq (array-like): Satellite channel center frequency (Hz). - - Returns: - numpy.ndarray: Observation type array. - - Raises: - ValueError: If any `obstype` is unassigned. - """ - - obstype = swcm.copy() - - # Use numpy vectorized operations - obstype = np.where(swcm == 5, 247, obstype) # WVCA/DL - obstype = np.where(swcm == 3, 246, obstype) # WVCT - obstype = np.where(swcm == 2, 251, obstype) # VIS - obstype = np.where(swcm == 1, 245, obstype) # IRLW - - condition = np.logical_and(swcm == 1, chanfreq >= 5e13) # IRSW - obstype = np.where(condition, 240, obstype) - - if not np.any(np.isin(obstype, [247, 246, 251, 245, 240])): - raise ValueError("Error: Unassigned ObsType found ... ") - - return obstype +import bufr +from bufr.obs_builder import add_main_functions +from bufr_satwnd_amv_obs_builder import SatWndAmvObsBuilder, map_path -def _make_obs(comm, input_path, mapping_path): - - # Get container from mapping file first - logging(comm, 'INFO', 'Get container from bufr') - container = bufr.Parser(input_path, mapping_path).parse(comm) - - logging(comm, 'DEBUG', f'container list (original): {container.list()}') - logging(comm, 'DEBUG', f'all_sub_categories = {container.all_sub_categories()}') - logging(comm, 'DEBUG', f'category map = {container.get_category_map()}') - - # Add new/derived data into container - for cat in container.all_sub_categories(): - - logging(comm, 'DEBUG', f'category = {cat}') - - satid = container.get('variables/satelliteId', cat) - if satid.size == 0: - logging(comm, 'WARNING', f'category {cat[0]} does not exist in input file') - paths = container.get_paths('variables/windComputationMethod', cat) - obstype = container.get('variables/windComputationMethod', cat) - container.add('variables/obstype_uwind', obstype, paths, cat) - container.add('variables/obstype_vwind', obstype, paths, cat) - - paths = container.get_paths('variables/windSpeed', cat) - wob = container.get('variables/windSpeed', cat) - container.add('variables/windEastward', wob, paths, cat) - container.add('variables/windNorthward', wob, paths, cat) - - else: - # Add new variables: ObsType/windEastward & ObsType/windNorthward - swcm = container.get('variables/windComputationMethod', cat) - chanfreq = container.get('variables/sensorCentralFrequency', cat) - - logging(comm, 'DEBUG', f'swcm min/max = {swcm.min()} {swcm.max()}') - logging(comm, 'DEBUG', f'chanfreq min/max = {chanfreq.min()} {chanfreq.max()}') - - obstype = _get_obs_type(swcm, chanfreq) - - logging(comm, 'DEBUG', f'obstype = {obstype}') - logging(comm, 'DEBUG', f'obstype min/max = {obstype.min()} {obstype.max()}') - - paths = container.get_paths('variables/windComputationMethod', cat) - container.add('variables/obstype_uwind', obstype, paths, cat) - container.add('variables/obstype_vwind', obstype, paths, cat) - - # Add new variables: ObsValue/windEastward & ObsValue/windNorthward - wdir = container.get('variables/windDirection', cat) - wspd = container.get('variables/windSpeed', cat) - - logging(comm, 'DEBUG', f'wdir min/max = {wdir.min()} {wdir.max()}') - logging(comm, 'DEBUG', f'wspd min/max = {wspd.min()} {wspd.max()}') - - uob, vob = compute_wind_components(wdir, wspd) - - logging(comm, 'DEBUG', f'uob min/max = {uob.min()} {uob.max()}') - logging(comm, 'DEBUG', f'vob min/max = {vob.min()} {vob.max()}') - - paths = container.get_paths('variables/windSpeed', cat) - container.add('variables/windEastward', uob, paths, cat) - container.add('variables/windNorthward', vob, paths, cat) - - # Check - logging(comm, 'DEBUG', f'container list (updated): {container.list()}') - logging(comm, 'DEBUG', f'all_sub_categories {container.all_sub_categories()}') - - return container - - -def create_obs_group(input_path, mapping_path, category, env): - - comm = bufr.mpi.Comm(env["comm_name"]) - - description = _make_description(mapping_path, update=True) - - # Check the cache for the data and return it if it exists - logging(comm, 'DEBUG', f'Check if bufr.DataCache exists? {bufr.DataCache.has(input_path, mapping_path)}') - if bufr.DataCache.has(input_path, mapping_path): - container = bufr.DataCache.get(input_path, mapping_path) - logging(comm, 'INFO', f'Encode {category} from cache') - data = iodaEncoder(description).encode(container)[(category,)] - logging(comm, 'INFO', f'Mark {category} as finished in the cache') - bufr.DataCache.mark_finished(input_path, mapping_path, [category]) - logging(comm, 'INFO', f'Return the encoded data for {category}') - return data - - container = _make_obs(comm, input_path, mapping_path) - - # Gather data from all tasks into all tasks. Each task will have the complete record - logging(comm, 'INFO', f'Gather data from all tasks into all tasks') - container.all_gather(comm) - - logging(comm, 'INFO', f'Add container to cache') - # Add the container to the cache - bufr.DataCache.add(input_path, mapping_path, container.all_sub_categories(), container) - - # Encode the data - logging(comm, 'INFO', f'Encode {category}') - data = iodaEncoder(description).encode(container)[(category,)] - - logging(comm, 'INFO', f'Mark {category} as finished in the cache') - # Mark the data as finished in the cache - bufr.DataCache.mark_finished(input_path, mapping_path, [category]) - - logging(comm, 'INFO', f'Return the encoded data for {category}') - return data - - -def create_obs_file(input_path, mapping_path, output_path): - - comm = bufr.mpi.Comm("world") - container = _make_obs(comm, input_path, mapping_path) - container.gather(comm) - - description = _make_description(mapping_path, update=True) - # Encode the data - if comm.rank() == 0: - netcdfEncoder(description).encode(container, output_path) +MAPPING_PATH = map_path('bufr_satwnd_amv_abi.yaml') - logging(comm, 'INFO', f'Return the encoded data') +class SatWndAmvAbiObsBuilder(SatWndAmvObsBuilder): + def __init__(self): + super().__init__(MAPPING_PATH, log_name=os.path.basename(__file__)) -if __name__ == '__main__': + def _get_obs_type(self, swcm, chanfreq): + obstype = swcm.copy() - start_time = time.time() + # Use numpy vectorized operations + obstype = np.where(swcm == 5, 247, obstype) # WVCA/DL + obstype = np.where(swcm == 3, 246, obstype) # WVCT + obstype = np.where(swcm == 2, 251, obstype) # VIS + obstype = np.where(swcm == 1, 245, obstype) # IRLW - bufr.mpi.App(sys.argv) - comm = bufr.mpi.Comm("world") + condition = np.logical_and(swcm == 1, chanfreq >= 5e13) # IRSW + obstype = np.where(condition, 240, obstype) - # Required input arguments as positional arguments - parser = argparse.ArgumentParser(description="Convert BUFR to NetCDF using a mapping file.") - parser.add_argument('input', type=str, help='Input BUFR file') - parser.add_argument('mapping', type=str, help='BUFR2IODA Mapping File') - parser.add_argument('output', type=str, help='Output NetCDF file') + if not np.any(np.isin(obstype, [247, 246, 251, 245, 240])): + raise ValueError("Error: Unassigned ObsType found ... ") - args = parser.parse_args() - infile = args.input - mapping = args.mapping - output = args.output + return obstype - create_obs_file(infile, mapping, output) - end_time = time.time() - running_time = end_time - start_time - logging(comm, 'INFO', f'Total running time: {running_time}') +# Add main functions create_obs_file and create_obs_group +add_main_functions(SatWndAmvAbiObsBuilder) diff --git a/dump/mapping/bufr_satwnd_amv_abi_mapping.yaml b/dump/mapping/bufr_satwnd_amv_abi.yaml similarity index 98% rename from dump/mapping/bufr_satwnd_amv_abi_mapping.yaml rename to dump/mapping/bufr_satwnd_amv_abi.yaml index 410d15f..2b133f8 100755 --- a/dump/mapping/bufr_satwnd_amv_abi_mapping.yaml +++ b/dump/mapping/bufr_satwnd_amv_abi.yaml @@ -161,7 +161,7 @@ encoder: source: variables/dataProviderOrigin longName: "Identification of Originating/Generating Center" - - name: "MetaData/qualityInformationWithoutForecast" + - name: "MetaData/qiWithoutForecast" source: variables/qualityInformationWithoutForecast longName: "Quality Information Without Forecast" diff --git a/dump/mapping/bufr_satwnd_amv_ahi.py b/dump/mapping/bufr_satwnd_amv_ahi.py index ba6f5d7..a659abc 100755 --- a/dump/mapping/bufr_satwnd_amv_ahi.py +++ b/dump/mapping/bufr_satwnd_amv_ahi.py @@ -1,342 +1,50 @@ #!/usr/bin/env python3 -import sys + import os -import argparse -import time import numpy as np -import bufr -from pyioda.ioda.Engines.Bufr import Encoder as iodaEncoder -from bufr.encoders.netcdf import Encoder as netcdfEncoder -from wxflow import Logger - -# Initialize Logger -# Get log level from the environment variable, default to 'INFO it not set -log_level = os.getenv('LOG_LEVEL', 'INFO') -logger = Logger('BUFR2IODA_satwnd_amv_ahi.py', level=log_level, colored_log=False) - - -def logging(comm, level, message): - - if comm.rank() == 0: - # Define a dictionary to map levels to logger methods - log_methods = { - 'DEBUG': logger.debug, - 'INFO': logger.info, - 'WARNING': logger.warning, - 'ERROR': logger.error, - 'CRITICAL': logger.critical, - } - - # Get the appropriate logging method, default to 'INFO' - log_method = log_methods.get(level.upper(), logger.info) - - if log_method == logger.info and level.upper() not in log_methods: - # Log a warning if the level is invalid - logger.warning(f'log level = {level}: not a valid level --> set to INFO') - - # Call the logging method - log_method(message) - - -def _make_description(mapping_path, update=False): - description = bufr.encoders.Description(mapping_path) - - if update: - # Define the variables to be added in a list of dictionaries - variables = [ - { - 'name': 'ObsType/windEastward', - 'source': 'variables/obstype_uwind', - 'units': '1', - 'longName': 'Observation Type based on Satellite-derived Wind Computation Method and Spectral Band', - }, - { - 'name': 'ObsType/windNorthward', - 'source': 'variables/obstype_vwind', - 'units': '1', - 'longName': 'Observation Type based on Satellite-derived Wind Computation Method and Spectral Band', - }, - { - 'name': 'ObsValue/windEastward', - 'source': 'variables/windEastward', - 'units': 'm s-1', - 'longName': 'Eastward Wind Component', - }, - { - 'name': 'ObsValue/windNorthward', - 'source': 'variables/windNorthward', - 'units': 'm s-1', - 'longName': 'Northward Wind Component', - }, - # MetaData/windGeneratingApplication will be inferred from variables/generatingApplication - # following a search for the proper variables/generatingApplication column - { - 'name': 'MetaData/windGeneratingApplication', - 'source': 'variables/windGeneratingApplication', - 'units': '1', - 'longName': 'Wind Generating Application', - }, - # MetaData/qualityInformationWithoutForecast will be inferred from variables/qualityInformation - # following a search for the proper variables/generatingApplication column - { - 'name': 'MetaData/qualityInformationWithoutForecast', - 'source': 'variables/qualityInformationWithoutForecast', - 'units': 'percent', - 'longName': 'Quality Information Without Forecast', - } - ] - - # Loop through each variable and add it to the description - for var in variables: - description.add_variable( - name=var['name'], - source=var['source'], - units=var['units'], - longName=var['longName'] - ) - - return description - - -def compute_wind_components(wdir, wspd): - """ - Compute the U and V wind components from wind direction and wind speed. - - Parameters: - wdir (array-like): Wind direction in degrees (meteorological convention: 0° = North, 90° = East). - wspd (array-like): Wind speed. - - Returns: - tuple: U and V wind components as numpy arrays with dtype float32. - """ - wdir_rad = np.radians(wdir) # Convert degrees to radians - u = -wspd * np.sin(wdir_rad) - v = -wspd * np.cos(wdir_rad) - - return u.astype(np.float32), v.astype(np.float32) - - -def _get_quality_information_and_generating_application(comm, gnap2D, pccf2D): - # For Himawari AHI data, qi w/o forecast (qifn) is packaged in same - # vector of qi with ga = 102 (EUMETSAT QI without forecast). Must - # conduct a search and extract the correct vector for gnap and qi - findQI = 102 - # 1. Find dimension-sizes of ga and qi (should be the same!) - gDim1, gDim2 = np.shape(gnap2D) - qDim1, qDim2 = np.shape(pccf2D) - logging(comm, 'INFO', 'Generating Application and Quality Information SEARCH') - logging(comm, 'DEBUG', f'Dimension size of GNAP ({gDim1},{gDim2})') - logging(comm, 'DEBUG', f'Dimension size of PCCF ({qDim1},{qDim2})') - # 2. Initialize gnap and qifn as None, and search for dimension of - # ga with values of 5. If the same column exists for qi, assign - # gnap to ga[:,i] and qifn to qi[:,i], else raise warning that no - # appropriate GNAP/PCCF combination was found - gnap = None - qifn = None - for i in range(gDim2): - if np.unique(gnap2D[:, i].squeeze()) == findQI: - if i <= qDim2: - logging(comm, 'INFO', f'GNAP/PCCF found for column {i}') - gnap = gnap2D[:, i].squeeze() - qifn = pccf2D[:, i].squeeze() - else: - logging(comm, 'INFO', f'ERROR: GNAP column {i} outside of PCCF dimension {qDim2}') - if (gnap is None) & (qifn is None): - raise ValueError(f'GNAP == {findQI} NOT FOUND OR OUT OF PCCF DIMENSION-RANGE, WILL FAIL!') - # If EE is needed, key search on np.unique(gnap2D[:,i].squeeze()) == 7 instead - # NOTE: Make sure to return np.float32 or np.int32 types as appropriate!!! - return gnap.astype(np.int32), qifn.astype(np.int32) - - -def _get_obs_type(swcm): - """ - Determine the observation type based on `swcm` and `chanfreq`. - - Parameters: - swcm (array-like): Switch mode values. - chanfreq (array-like): Channel frequency values (Hz). - - Returns: - numpy.ndarray: Observation type array. - - Raises: - ValueError: If any `obstype` is unassigned. - """ - - obstype = swcm.copy() - - # Use numpy vectorized operations - obstype = np.where(swcm == 5, 250, obstype) # WVCA/DL - obstype = np.where(swcm == 3, 250, obstype) # WVCT - obstype = np.where(swcm == 2, 242, obstype) # VIS - obstype = np.where(swcm == 1, 252, obstype) # IRLW - - if not np.any(np.isin(obstype, [242, 250, 252])): - raise ValueError("Error: Unassigned ObsType found ... ") - - return obstype - - -def _make_obs(comm, input_path, mapping_path): - - # Get container from mapping file first - logging(comm, 'INFO', 'Get container from bufr') - container = bufr.Parser(input_path, mapping_path).parse(comm) - logging(comm, 'DEBUG', f'container list (original): {container.list()}') - logging(comm, 'DEBUG', f'all_sub_categories = {container.all_sub_categories()}') - logging(comm, 'DEBUG', f'category map = {container.get_category_map()}') - - # Add new/derived data into container - for cat in container.all_sub_categories(): - - logging(comm, 'DEBUG', f'category = {cat}') - - satid = container.get('variables/satelliteId', cat) - if satid.size == 0: - logging(comm, 'WARNING', f'category {cat[0]} does not exist in input file') - paths = container.get_paths('variables/windComputationMethod', cat) - obstype = container.get('variables/windComputationMethod', cat) - container.add('variables/obstype_uwind', obstype, paths, cat) - container.add('variables/obstype_vwind', obstype, paths, cat) - - paths = container.get_paths('variables/windSpeed', cat) - wob = container.get('variables/windSpeed', cat) - container.add('variables/windEastward', wob, paths, cat) - container.add('variables/windNorthward', wob, paths, cat) - - paths = container.get_paths('variables/windComputationMethod', cat) - dummy = container.get('variables/windSpeed', cat) - container.add('variables/windGeneratingApplication', dummy, paths, cat) - container.add('variables/qualityInformationWithoutForecast', dummy, paths, cat) - - else: - # Add new variables: ObsType/windEastward & ObsType/windNorthward - swcm = container.get('variables/windComputationMethod', cat) - chanfreq = container.get('variables/sensorCentralFrequency', cat) - - logging(comm, 'DEBUG', f'swcm min/max = {swcm.min()} {swcm.max()}') - logging(comm, 'DEBUG', f'chanfreq min/max = {chanfreq.min()} {chanfreq.max()}') - - obstype = _get_obs_type(swcm) - - logging(comm, 'DEBUG', f'obstype = {obstype}') - logging(comm, 'DEBUG', f'obstype min/max = {obstype.min()} {obstype.max()}') - - paths = container.get_paths('variables/windComputationMethod', cat) - container.add('variables/obstype_uwind', obstype, paths, cat) - container.add('variables/obstype_vwind', obstype, paths, cat) - - # Add new variables: ObsValue/windEastward & ObsValue/windNorthward - wdir = container.get('variables/windDirection', cat) - wspd = container.get('variables/windSpeed', cat) - - logging(comm, 'DEBUG', f'wdir min/max = {wdir.min()} {wdir.max()}') - logging(comm, 'DEBUG', f'wspd min/max = {wspd.min()} {wspd.max()}') - - uob, vob = compute_wind_components(wdir, wspd) - - logging(comm, 'DEBUG', f'uob min/max = {uob.min()} {uob.max()}') - logging(comm, 'DEBUG', f'vob min/max = {vob.min()} {vob.max()}') - - paths = container.get_paths('variables/windSpeed', cat) - container.add('variables/windEastward', uob, paths, cat) - container.add('variables/windNorthward', vob, paths, cat) - - # Add new variables: MetaData/windGeneratingApplication and qualityInformationWithoutForecast - gnap2D = container.get('variables/generatingApplication', cat) - pccf2D = container.get('variables/qualityInformation', cat) - - gnap, qifn = _get_quality_information_and_generating_application(comm, gnap2D, pccf2D) - - logging(comm, 'DEBUG', f'gnap min/max = {gnap.min()} {gnap.max()}') - logging(comm, 'DEBUG', f'qifn min/max = {qifn.min()} {qifn.max()}') - - paths = container.get_paths('variables/windComputationMethod', cat) - container.add('variables/windGeneratingApplication', gnap, paths, cat) - container.add('variables/qualityInformationWithoutForecast', qifn, paths, cat) - - # Check - logging(comm, 'DEBUG', f'container list (updated): {container.list()}') - logging(comm, 'DEBUG', f'all_sub_categories {container.all_sub_categories()}') - - return container - - -def create_obs_group(input_path, mapping_path, category, env): - - comm = bufr.mpi.Comm(env["comm_name"]) - - description = _make_description(mapping_path, update=True) - - # Check the cache for the data and return it if it exists - logging(comm, 'DEBUG', f'Check if bufr.DataCache exists? {bufr.DataCache.has(input_path, mapping_path)}') - if bufr.DataCache.has(input_path, mapping_path): - container = bufr.DataCache.get(input_path, mapping_path) - logging(comm, 'INFO', f'Encode {category} from cache') - data = iodaEncoder(description).encode(container)[(category,)] - logging(comm, 'INFO', f'Mark {category} as finished in the cache') - bufr.DataCache.mark_finished(input_path, mapping_path, [category]) - logging(comm, 'INFO', f'Return the encoded data for {category}') - return data - - container = _make_obs(comm, input_path, mapping_path) - - # Gather data from all tasks into all tasks. Each task will have the complete record - logging(comm, 'INFO', f'Gather data from all tasks into all tasks') - container.all_gather(comm) - - logging(comm, 'INFO', f'Add container to cache') - # Add the container to the cache - bufr.DataCache.add(input_path, mapping_path, container.all_sub_categories(), container) - - # Encode the data - logging(comm, 'INFO', f'Encode {category}') - data = iodaEncoder(description).encode(container)[(category,)] - - logging(comm, 'INFO', f'Mark {category} as finished in the cache') - # Mark the data as finished in the cache - bufr.DataCache.mark_finished(input_path, mapping_path, [category]) +import bufr +from bufr.obs_builder import add_main_functions +from bufr_satwnd_amv_obs_builder import SatWndAmvObsBuilder, map_path - logging(comm, 'INFO', f'Return the encoded data for {category}') - return data +MAPPING_PATH = map_path('bufr_satwnd_amv_ahi.yaml') +FIND_QI = 102 -def create_obs_file(input_path, mapping_path, output_path): - comm = bufr.mpi.Comm("world") - container = _make_obs(comm, input_path, mapping_path) - container.gather(comm) +class SatWndAmvAhiObsBuilder(SatWndAmvObsBuilder): + def __init__(self): + super().__init__(MAPPING_PATH, log_name=os.path.basename(__file__)) - description = _make_description(mapping_path, update=True) + def make_obs(self, comm, input_path): + # Get container from mapping file first + container = super().make_obs(comm, input_path) - # Encode the data - if comm.rank() == 0: - netcdfEncoder(description).encode(container, output_path) + # Add new/derived data into container + for cat in container.all_sub_categories(): + self._add_quality_info_and_gen_app(FIND_QI, container, cat) - logging(comm, 'INFO', f'Return the encoded data') + return container + def _make_description(self): + description = super()._make_description() + self._add_quality_info_and_gen_app_descriptions(description) -if __name__ == '__main__': + return description - start_time = time.time() + def _get_obs_type(self, swcm, chan_freq): + obstype = swcm.copy() - bufr.mpi.App(sys.argv) - comm = bufr.mpi.Comm("world") + # Use numpy vectorized operations + obstype = np.where(swcm == 5, 250, obstype) # WVCA/DL + obstype = np.where(swcm == 3, 250, obstype) # WVCT + obstype = np.where(swcm == 2, 242, obstype) # VIS + obstype = np.where(swcm == 1, 252, obstype) # IRLW - # Required input arguments as positional arguments - parser = argparse.ArgumentParser(description="Convert BUFR to NetCDF using a mapping file.") - parser.add_argument('input', type=str, help='Input BUFR file') - parser.add_argument('mapping', type=str, help='BUFR2IODA Mapping File') - parser.add_argument('output', type=str, help='Output NetCDF file') + if not np.any(np.isin(obstype, [242, 250, 252])): + raise ValueError("Error: Unassigned ObsType found ... ") - args = parser.parse_args() - infile = args.input - mapping = args.mapping - output = args.output + return obstype - create_obs_file(infile, mapping, output) - end_time = time.time() - running_time = end_time - start_time - logging(comm, 'INFO', f'Total running time: {running_time}') +add_main_functions(SatWndAmvAhiObsBuilder) diff --git a/dump/mapping/bufr_satwnd_amv_ahi_mapping.yaml b/dump/mapping/bufr_satwnd_amv_ahi.yaml similarity index 100% rename from dump/mapping/bufr_satwnd_amv_ahi_mapping.yaml rename to dump/mapping/bufr_satwnd_amv_ahi.yaml diff --git a/dump/mapping/bufr_satwnd_amv_avhrr.py b/dump/mapping/bufr_satwnd_amv_avhrr.py index dbd9983..5436388 100755 --- a/dump/mapping/bufr_satwnd_amv_avhrr.py +++ b/dump/mapping/bufr_satwnd_amv_avhrr.py @@ -1,369 +1,130 @@ #!/usr/bin/env python3 -import sys + import os -import argparse -import time import numpy as np + import bufr -from pyioda.ioda.Engines.Bufr import Encoder as iodaEncoder -from bufr.encoders.netcdf import Encoder as netcdfEncoder -from wxflow import Logger +from bufr.obs_builder import add_main_functions +from bufr_satwnd_amv_obs_builder import SatWndAmvObsBuilder, map_path -# Initialize Logger -# Get log level from the environment variable, default to 'INFO it not set -log_level = os.getenv('LOG_LEVEL', 'INFO') -logger = Logger('BUFR2IODA_satwnd_amv_avhrr.py', level=log_level, colored_log=False) +MAPPING_PATH = map_path('bufr_satwnd_amv_avhrr.yaml') -def logging(comm, level, message): - if comm.rank() == 0: - # Define a dictionary to map levels to logger methods - log_methods = { - 'DEBUG': logger.debug, - 'INFO': logger.info, - 'WARNING': logger.warning, - 'ERROR': logger.error, - 'CRITICAL': logger.critical, - } +class SatWndAmvAvhrrObsBuilder(SatWndAmvObsBuilder): + def __init__(self): + super().__init__(MAPPING_PATH, log_name=os.path.basename(__file__)) - # Get the appropriate logging method, default to 'INFO' - log_method = log_methods.get(level.upper(), logger.info) + def make_obs(self, comm, input_path): + # Get container from mapping file first + container = super().make_obs(comm, input_path) - if log_method == logger.info and level.upper() not in log_methods: - # Log a warning if the level is invalid - logger.warning(f'log level = {level}: not a valid level --> set to INFO') + # Add new/derived data into container + for cat in container.all_sub_categories(): + self._add_avhrr_quality_info_and_gen_app(container, cat) - # Call the logging method - log_method(message) + return container + def _make_description(self): + description = super()._make_description() + self._add_quality_info_and_gen_app_descriptions(description) -def _make_description(mapping_path, update=False): - description = bufr.encoders.Description(mapping_path) + return description - if update: - # Define the variables to be added in a list of dictionaries - variables = [ - { - 'name': 'ObsType/windEastward', - 'source': 'variables/obstype_uwind', - 'units': '1', - 'longName': 'Observation Type based on Satellite-derived Wind Computation Method and Spectral Band', - }, - { - 'name': 'ObsType/windNorthward', - 'source': 'variables/obstype_vwind', - 'units': '1', - 'longName': 'Observation Type based on Satellite-derived Wind Computation Method and Spectral Band', - }, - { - 'name': 'ObsValue/windEastward', - 'source': 'variables/windEastward', - 'units': 'm s-1', - 'longName': 'Eastward Wind Component', - }, - { - 'name': 'ObsValue/windNorthward', - 'source': 'variables/windNorthward', - 'units': 'm s-1', - 'longName': 'Northward Wind Component', - }, - # MetaData/windGeneratingApplication will be inferred from variables/generatingApplication - # following a search for the proper variables/generatingApplication column - { - 'name': 'MetaData/windGeneratingApplication', - 'source': 'variables/windGeneratingApplication', - 'units': '1', - 'longName': 'Wind Generating Application', - }, - # MetaData/qualityInformationWithoutForecast will be inferred from variables/qualityInformation - # following a search for the proper variables/generatingApplication column - { - 'name': 'MetaData/qualityInformationWithoutForecast', - 'source': 'variables/qualityInformationWithoutForecast', - 'units': 'percent', - 'longName': 'Quality Information Without Forecast', - } - ] + def _get_obs_type(self, swcm, chan_freq): + obstype = swcm.copy() - # Loop through each variable and add it to the description - for var in variables: - description.add_variable( - name=var['name'], - source=var['source'], - units=var['units'], - longName=var['longName'] - ) + # Use numpy vectorized operations + obstype = np.where(swcm == 1, 244, obstype) # IRLW - return description + if not np.any(np.isin(obstype, [244])): + raise ValueError("Error: Unassigned ObsType found ... ") + return obstype.astype(np.int32) -def compute_wind_components(wdir, wspd): - """ - Compute the U and V wind components from wind direction and wind speed. + def _add_avhrr_quality_info_and_gen_app(self, container, cat): + # Add new variables: MetaData/windGeneratingApplication and qualityInformationWithoutForecast + gnap2D = container.get('generatingApplication', cat) + pccf2D = container.get('qualityInformation', cat) + satId = container.get('satelliteId', cat) - Parameters: - wdir (array-like): Wind direction in degrees (meteorological convention: 0° = North, 90° = East). - wspd (array-like): Wind speed. + if not satId.size: + paths = container.get_paths('windComputationMethod', cat) + dummy = container.get('windSpeed', cat) + container.add('windGeneratingApplication', dummy, paths, cat) + container.add('qualityInformationWithoutForecast', dummy, paths, cat) + return - Returns: - tuple: U and V wind components as numpy arrays with dtype float32. - """ - wdir_rad = np.radians(wdir) # Convert degrees to radians - u = -wspd * np.sin(wdir_rad) - v = -wspd * np.cos(wdir_rad) + gnap, qifn = self._get_avhrr_quality_info_and_gen_app(gnap2D, pccf2D, satId) - return u.astype(np.float32), v.astype(np.float32) + self.log.debug(f'gnap min/max = {gnap.min()} {gnap.max()}') + self.log.debug(f'qifn min/max = {qifn.min()} {qifn.max()}') + paths = container.get_paths('windComputationMethod', cat) + container.add('windGeneratingApplication', gnap, paths, cat) + container.add('qualityInformationWithoutForecast', qifn, paths, cat) -def _get_quality_information_and_generating_application(comm, gnap2D, pccf2D, satID): - # For METOP-A/B/C AVHRR data (satID 3,4,5), qi w/o forecast (qifn) is - # packaged in same vector of qi with ga = 5 (QI without forecast), and EE - # is packaged in same vector of qi with ga=7 (Estimated Error (EE) in m/s - # converted to a percent confidence) shape (4,nobs). - # - # For NOAA-15/18/19 AVHRR data (satID 206,209,223), qi w/o forecast - # (qifn) is packaged in same vector of qi with ga = 1 (EUMETSAT QI - # without forecast), and EE is packaged in same vector of qi with ga=4 - # (Estimated Error (EE) in m/s converted to a percent confidence) shape - # (4,nobs). - # - # Must conduct a search and extract the correct vector for gnap and qi - # 0. Define the appropriate QI and EE search values, based on satID - if np.all(np.isin(satID, [206, 209, 223])): # NESDIS AVHRR set - findQI = 1 - findEE = 4 - elif np.all(np.isin(satID, [3, 4, 5])): # EUMETSAT AVHRR set - findQI = 5 - findEE = 7 - # There is a catch: prior to 2023 AVHRR winds from EUMETSAT were formatted in the same - # way as NESDIS AVHRR winds and both were passed through the NC005080 tank as a single - # dataset. In that case, we need to actually set findQI=1 and findEE=4 here. - # Let's do a preliminary check to see if any gnap2D values match findQI. If not, let's - # automatically switch to findQI=1, findEE=4 and presume pre-2023 EUMETSAT AVHRR format - if not (np.any(np.isin(gnap2D, [findQI]))): - logging(comm, 'DEBUG', f'NO GNAP VALUE OF {findQI} EXISTS FOR EUMETSAT AVHRR DATASET, PRESUMING PRE-2023 FORMATTING') + def _get_avhrr_quality_info_and_gen_app(self, gnap2D, pccf2D, satID): + # For METOP-A/B/C AVHRR data (satID 3,4,5), qi w/o forecast (qifn) is + # packaged in same vector of qi with ga = 5 (QI without forecast), and EE + # is packaged in same vector of qi with ga=7 (Estimated Error (EE) in m/s + # converted to a percent confidence) shape (4,nobs). + # + # For NOAA-15/18/19 AVHRR data (satID 206,209,223), qi w/o forecast + # (qifn) is packaged in same vector of qi with ga = 1 (EUMETSAT QI + # without forecast), and EE is packaged in same vector of qi with ga=4 + # (Estimated Error (EE) in m/s converted to a percent confidence) shape + # (4,nobs). + # + # Must conduct a search and extract the correct vector for gnap and qi + # 0. Define the appropriate QI and EE search values, based on satID + if np.all(np.isin(satID, [206, 209, 223])): # NESDIS AVHRR set findQI = 1 findEE = 4 - else: - logging(comm, 'DEBUG', f'satID set not found (all satID values follow):') - for sid in np.unique(satID): - logging(comm, 'DEBUG', f'satID: {sid}') - logging(comm, 'DEBUG', f'BTH: findQI={findQI}') - # 1. Find dimension-sizes of ga and qi (should be the same!) - gDim1, gDim2 = np.shape(gnap2D) - qDim1, qDim2 = np.shape(pccf2D) - logging(comm, 'INFO', f'Generating Application and Quality Information SEARCH:') - logging(comm, 'DEBUG', f'Dimension size of GNAP ({gDim1},{gDim2})') - logging(comm, 'DEBUG', f'Dimension size of PCCF ({qDim1},{qDim2})') - # 2. Initialize gnap and qifn as None, and search for dimension of - # ga with values of findQI. If the same column exists for qi, assign - # gnap to ga[:,i] and qifn to qi[:,i], else raise warning that no - # appropriate GNAP/PCCF combination was found - gnap = None - qifn = None - for i in range(gDim2): - if np.unique(gnap2D[:, i].squeeze()) == findQI: - if i <= qDim2: - logging(comm, 'INFO', f'GNAP/PCCF found for column {i}') - gnap = gnap2D[:, i].squeeze() - qifn = pccf2D[:, i].squeeze() - else: - logging(comm, 'INFO', f'ERROR: GNAP column {i} outside of PCCF dimension {qDim2}') - if (gnap is None) & (qifn is None): - raise ValueError(f'GNAP == {findQI} NOT FOUND OR OUT OF PCCF DIMENSION-RANGE, WILL FAIL!') - # If EE is needed, key search on np.unique(gnap2D[:,i].squeeze()) == findEE instead - # NOTE: Make sure to return np.float32 or np.int32 types as appropriate!!! - return gnap.astype(np.int32), qifn.astype(np.int32) - - -def _get_obs_type(swcm): - """ - Determine the observation type based on `swcm` and `chanfreq`. - - Parameters: - swcm (array-like): Switch mode values. - chanfreq (array-like): Channel frequency values (Hz). - - Returns: - numpy.ndarray: Observation type array. - - Raises: - ValueError: If any `obstype` is unassigned. - """ - - obstype = swcm.copy() - - # Use numpy vectorized operations - obstype = np.where(swcm == 1, 244, obstype) # IRLW - - if not np.any(np.isin(obstype, [244])): - raise ValueError("Error: Unassigned ObsType found ... ") - - return obstype.astype(np.int32) - - -def _make_obs(comm, input_path, mapping_path): - - # Get container from mapping file first - logging(comm, 'INFO', 'Get container from bufr') - container = bufr.Parser(input_path, mapping_path).parse(comm) - - logging(comm, 'DEBUG', f'container list (original): {container.list()}') - logging(comm, 'DEBUG', f'all_sub_categories = {container.all_sub_categories()}') - logging(comm, 'DEBUG', f'category map = {container.get_category_map()}') - - # Add new/derived data into container - for cat in container.all_sub_categories(): - - logging(comm, 'DEBUG', f'category = {cat}') - - satid = container.get('variables/satelliteId', cat) - if satid.size == 0: - logging(comm, 'WARNING', f'category {cat[0]} does not exist in input file') - paths = container.get_paths('variables/windComputationMethod', cat) - obstype = container.get('variables/windComputationMethod', cat) - container.add('variables/obstype_uwind', obstype, paths, cat) - container.add('variables/obstype_vwind', obstype, paths, cat) - - paths = container.get_paths('variables/windSpeed', cat) - wob = container.get('variables/windSpeed', cat) - container.add('variables/windEastward', wob, paths, cat) - container.add('variables/windNorthward', wob, paths, cat) - - paths = container.get_paths('variables/windComputationMethod', cat) - dummy = container.get('variables/windSpeed', cat) - container.add('variables/windGeneratingApplication', dummy, paths, cat) - container.add('variables/qualityInformationWithoutForecast', dummy, paths, cat) - + elif np.all(np.isin(satID, [3, 4, 5])): # EUMETSAT AVHRR set + findQI = 5 + findEE = 7 + # There is a catch: prior to 2023 AVHRR winds from EUMETSAT were formatted in the same + # way as NESDIS AVHRR winds and both were passed through the NC005080 tank as a single + # dataset. In that case, we need to actually set findQI=1 and findEE=4 here. + # Let's do a preliminary check to see if any gnap2D values match findQI. If not, let's + # automatically switch to findQI=1, findEE=4 and presume pre-2023 EUMETSAT AVHRR format + if not np.any(np.isin(gnap2D, [findQI])): + self.log.debug( + f'NO GNAP VALUE OF {findQI} EXISTS FOR EUMETSAT AVHRR DATASET, PRESUMING PRE-2023 FORMATTING') + findQI = 1 + findEE = 4 else: - # Add new variables: ObsType/windEastward & ObsType/windNorthward - swcm = container.get('variables/windComputationMethod', cat) - chanfreq = container.get('variables/sensorCentralFrequency', cat) - - logging(comm, 'DEBUG', f'swcm min/max = {swcm.min()} {swcm.max()}') - logging(comm, 'DEBUG', f'chanfreq min/max = {chanfreq.min()} {chanfreq.max()}') - - obstype = _get_obs_type(swcm) - - logging(comm, 'DEBUG', f'obstype = {obstype}') - logging(comm, 'DEBUG', f'obstype min/max = {obstype.min()} {obstype.max()}') - - paths = container.get_paths('variables/windComputationMethod', cat) - container.add('variables/obstype_uwind', obstype, paths, cat) - container.add('variables/obstype_vwind', obstype, paths, cat) - - # Add new variables: ObsValue/windEastward & ObsValue/windNorthward - wdir = container.get('variables/windDirection', cat) - wspd = container.get('variables/windSpeed', cat) - - logging(comm, 'DEBUG', f'wdir min/max = {wdir.min()} {wdir.max()}') - logging(comm, 'DEBUG', f'wspd min/max = {wspd.min()} {wspd.max()}') - - uob, vob = compute_wind_components(wdir, wspd) - - logging(comm, 'DEBUG', f'uob min/max = {uob.min()} {uob.max()}') - logging(comm, 'DEBUG', f'vob min/max = {vob.min()} {vob.max()}') - - paths = container.get_paths('variables/windSpeed', cat) - container.add('variables/windEastward', uob, paths, cat) - container.add('variables/windNorthward', vob, paths, cat) - - # Add new variables: MetaData/windGeneratingApplication and qualityInformationWithoutForecast - satID = container.get('variables/satelliteId', cat) - gnap2D = container.get('variables/generatingApplication', cat) - pccf2D = container.get('variables/qualityInformation', cat) - - gnap, qifn = _get_quality_information_and_generating_application(comm, gnap2D, pccf2D, satID) - - logging(comm, 'DEBUG', f'gnap min/max = {gnap.min()} {gnap.max()}') - logging(comm, 'DEBUG', f'qifn min/max = {qifn.min()} {qifn.max()}') - - paths = container.get_paths('variables/windComputationMethod', cat) - container.add('variables/windGeneratingApplication', gnap, paths, cat) - container.add('variables/qualityInformationWithoutForecast', qifn, paths, cat) - - # Check - logging(comm, 'DEBUG', f'container list (updated): {container.list()}') - logging(comm, 'DEBUG', f'all_sub_categories {container.all_sub_categories()}') - - return container - - -def create_obs_group(input_path, mapping_path, category, env): - - comm = bufr.mpi.Comm(env["comm_name"]) - - description = _make_description(mapping_path, update=True) - - # Check the cache for the data and return it if it exists - logging(comm, 'DEBUG', f'Check if bufr.DataCache exists? {bufr.DataCache.has(input_path, mapping_path)}') - if bufr.DataCache.has(input_path, mapping_path): - container = bufr.DataCache.get(input_path, mapping_path) - logging(comm, 'INFO', f'Encode {category} from cache') - data = iodaEncoder(description).encode(container)[(category,)] - logging(comm, 'INFO', f'Mark {category} as finished in the cache') - bufr.DataCache.mark_finished(input_path, mapping_path, [category]) - logging(comm, 'INFO', f'Return the encoded data for {category}') - return data - - container = _make_obs(comm, input_path, mapping_path) - - # Gather data from all tasks into all tasks. Each task will have the complete record - logging(comm, 'INFO', f'Gather data from all tasks into all tasks') - container.all_gather(comm) - - logging(comm, 'INFO', f'Add container to cache') - # Add the container to the cache - bufr.DataCache.add(input_path, mapping_path, container.all_sub_categories(), container) - - # Encode the data - logging(comm, 'INFO', f'Encode {category}') - data = iodaEncoder(description).encode(container)[(category,)] - - logging(comm, 'INFO', f'Mark {category} as finished in the cache') - # Mark the data as finished in the cache - bufr.DataCache.mark_finished(input_path, mapping_path, [category]) - - logging(comm, 'INFO', f'Return the encoded data for {category}') - return data - - -def create_obs_file(input_path, mapping_path, output_path): - - comm = bufr.mpi.Comm("world") - container = _make_obs(comm, input_path, mapping_path) - container.gather(comm) - - description = _make_description(mapping_path, update=True) - - # Encode the data - if comm.rank() == 0: - netcdfEncoder(description).encode(container, output_path) - - logging(comm, 'INFO', f'Return the encoded data') - - -if __name__ == '__main__': - - start_time = time.time() - - bufr.mpi.App(sys.argv) - comm = bufr.mpi.Comm("world") - - # Required input arguments as positional arguments - parser = argparse.ArgumentParser(description="Convert BUFR to NetCDF using a mapping file.") - parser.add_argument('input', type=str, help='Input BUFR file') - parser.add_argument('mapping', type=str, help='BUFR2IODA Mapping File') - parser.add_argument('output', type=str, help='Output NetCDF file') - - args = parser.parse_args() - infile = args.input - mapping = args.mapping - output = args.output - - create_obs_file(infile, mapping, output) - - end_time = time.time() - running_time = end_time - start_time - logging(comm, 'INFO', f'Total running time: {running_time}') + self.log.debug(f'satID set not found (all satID values follow):') + for sid in np.unique(satID): + self.log.debug(f'satID: {sid}') + self.log.debug(f'BTH: findQI={findQI}') + # 1. Find dimension-sizes of ga and qi (should be the same!) + gDim1, gDim2 = np.shape(gnap2D) + qDim1, qDim2 = np.shape(pccf2D) + self.log.info(f'Generating Application and Quality Information SEARCH:') + self.log.debug(f'Dimension size of GNAP ({gDim1},{gDim2})') + self.log.debug(f'Dimension size of PCCF ({qDim1},{qDim2})') + # 2. Initialize gnap and qifn as None, and search for dimension of + # ga with values of findQI. If the same column exists for qi, assign + # gnap to ga[:,i] and qifn to qi[:,i], else raise warning that no + # appropriate GNAP/PCCF combination was found + gnap = None + qifn = None + for i in range(gDim2): + if np.unique(gnap2D[:, i].squeeze()) == findQI: + if i <= qDim2: + self.log.info(f'GNAP/PCCF found for column {i}') + gnap = gnap2D[:, i].squeeze() + qifn = pccf2D[:, i].squeeze() + else: + self.log.info(f'ERROR: GNAP column {i} outside of PCCF dimension {qDim2}') + if (gnap is None) and (qifn is None): + raise ValueError(f'GNAP == {findQI} NOT FOUND OR OUT OF PCCF DIMENSION-RANGE, WILL FAIL!') + # If EE is needed, key search on np.unique(gnap2D[:,i].squeeze()) == findEE instead + # NOTE: Make sure to return np.float32 or np.int32 types as appropriate!!! + return gnap.astype(np.int32), qifn.astype(np.int32) + + +# Add main functions create_obs_file and create_obs_group +add_main_functions(SatWndAmvAvhrrObsBuilder) diff --git a/dump/mapping/bufr_satwnd_amv_avhrr_mapping.yaml b/dump/mapping/bufr_satwnd_amv_avhrr.yaml similarity index 96% rename from dump/mapping/bufr_satwnd_amv_avhrr_mapping.yaml rename to dump/mapping/bufr_satwnd_amv_avhrr.yaml index 9723f27..b5554f3 100755 --- a/dump/mapping/bufr_satwnd_amv_avhrr_mapping.yaml +++ b/dump/mapping/bufr_satwnd_amv_avhrr.yaml @@ -67,12 +67,12 @@ bufr: category: variable: satelliteId map: - _3: metop-1b - _4: metop-1a - _5: metop-1c - _206: noaa-15 - _209: noaa-18 - _223: noaa-19 + _3: metop-b + _4: metop-a + _5: metop-c + _206: n15 + _209: n18 + _223: n19 encoder: # type: netcdf diff --git a/dump/mapping/bufr_satwnd_amv_leogeo.py b/dump/mapping/bufr_satwnd_amv_leogeo.py index f68ec23..a2bc08e 100755 --- a/dump/mapping/bufr_satwnd_amv_leogeo.py +++ b/dump/mapping/bufr_satwnd_amv_leogeo.py @@ -1,273 +1,30 @@ #!/usr/bin/env python3 -import sys + import os -import argparse -import time import numpy as np -import bufr -from pyioda.ioda.Engines.Bufr import Encoder as iodaEncoder -from bufr.encoders.netcdf import Encoder as netcdfEncoder -from wxflow import Logger - -# Initialize Logger -# Get log level from the environment variable, default to 'INFO it not set -log_level = os.getenv('LOG_LEVEL', 'INFO') -logger = Logger('BUFR2IODA_satwnd_amv_leogeo.py', level=log_level, colored_log=False) - - -def logging(comm, level, message): - - if comm.rank() == 0: - # Define a dictionary to map levels to logger methods - log_methods = { - 'DEBUG': logger.debug, - 'INFO': logger.info, - 'WARNING': logger.warning, - 'ERROR': logger.error, - 'CRITICAL': logger.critical, - } - - # Get the appropriate logging method, default to 'INFO' - log_method = log_methods.get(level.upper(), logger.info) - - if log_method == logger.info and level.upper() not in log_methods: - # Log a warning if the level is invalid - logger.warning(f'log level = {level}: not a valid level --> set to INFO') - - # Call the logging method - log_method(message) - - -def _make_description(mapping_path, update=False): - description = bufr.encoders.Description(mapping_path) - - if update: - # Define the variables to be added in a list of dictionaries - variables = [ - { - 'name': 'ObsType/windEastward', - 'source': 'variables/obstype_uwind', - 'units': '1', - 'longName': 'Observation Type based on Satellite-derived Wind Computation Method and Spectral Band', - }, - { - 'name': 'ObsType/windNorthward', - 'source': 'variables/obstype_vwind', - 'units': '1', - 'longName': 'Observation Type based on Satellite-derived Wind Computation Method and Spectral Band', - }, - { - 'name': 'ObsValue/windEastward', - 'source': 'variables/windEastward', - 'units': 'm s-1', - 'longName': 'Eastward Wind Component', - }, - { - 'name': 'ObsValue/windNorthward', - 'source': 'variables/windNorthward', - 'units': 'm s-1', - 'longName': 'Northward Wind Component', - } - ] - - # Loop through each variable and add it to the description - for var in variables: - description.add_variable( - name=var['name'], - source=var['source'], - units=var['units'], - longName=var['longName'] - ) - - return description - - -def compute_wind_components(wdir, wspd): - """ - Compute the U and V wind components from wind direction and wind speed. - - Parameters: - wdir (array-like): Wind direction in degrees (meteorological convention: 0° = North, 90° = East). - wspd (array-like): Wind speed. - - Returns: - tuple: U and V wind components as numpy arrays with dtype float32. - """ - wdir_rad = np.radians(wdir) # Convert degrees to radians - u = -wspd * np.sin(wdir_rad) - v = -wspd * np.cos(wdir_rad) - - return u.astype(np.float32), v.astype(np.float32) - - -def _get_obs_type(swcm): - """ - Determine the observation type based on `swcm` and `chanfreq`. - - Parameters: - swcm (array-like): Switch mode values. - chanfreq (array-like): Channel frequency values (Hz). - - Returns: - numpy.ndarray: Observation type array. - - Raises: - ValueError: If any `obstype` is unassigned. - """ - - obstype = swcm.copy() - - # Use numpy vectorized operations - obstype = np.where(swcm == 1, 255, obstype) # IRLW - - if not np.any(np.isin(obstype, [255])): - raise ValueError("Error: Unassigned ObsType found ... ") - - return obstype - - -def _make_obs(comm, input_path, mapping_path): - - # Get container from mapping file first - logging(comm, 'INFO', 'Get container from bufr') - container = bufr.Parser(input_path, mapping_path).parse(comm) - - logging(comm, 'DEBUG', f'container list (original): {container.list()}') - logging(comm, 'DEBUG', f'all_sub_categories = {container.all_sub_categories()}') - logging(comm, 'DEBUG', f'category map = {container.get_category_map()}') - # Add new/derived data into container - for cat in container.all_sub_categories(): - - logging(comm, 'DEBUG', f'category = {cat}') - - satid = container.get('variables/satelliteId', cat) - if satid.size == 0: - logging(comm, 'WARNING', f'category {cat[0]} does not exist in input file') - paths = container.get_paths('variables/windComputationMethod', cat) - obstype = container.get('variables/windComputationMethod', cat) - container.add('variables/obstype_uwind', obstype, paths, cat) - container.add('variables/obstype_vwind', obstype, paths, cat) - - paths = container.get_paths('variables/windSpeed', cat) - wob = container.get('variables/windSpeed', cat) - container.add('variables/windEastward', wob, paths, cat) - container.add('variables/windNorthward', wob, paths, cat) - - else: - # Add new variables: ObsType/windEastward & ObsType/windNorthward - swcm = container.get('variables/windComputationMethod', cat) - chanfreq = container.get('variables/sensorCentralFrequency', cat) - - logging(comm, 'DEBUG', f'swcm min/max = {swcm.min()} {swcm.max()}') - logging(comm, 'DEBUG', f'chanfreq min/max = {chanfreq.min()} {chanfreq.max()}') - - obstype = _get_obs_type(swcm) - - logging(comm, 'DEBUG', f'obstype = {obstype}') - logging(comm, 'DEBUG', f'obstype min/max = {obstype.min()} {obstype.max()}') - - paths = container.get_paths('variables/windComputationMethod', cat) - container.add('variables/obstype_uwind', obstype, paths, cat) - container.add('variables/obstype_vwind', obstype, paths, cat) - - # Add new variables: ObsValue/windEastward & ObsValue/windNorthward - wdir = container.get('variables/windDirection', cat) - wspd = container.get('variables/windSpeed', cat) - - logging(comm, 'DEBUG', f'wdir min/max = {wdir.min()} {wdir.max()}') - logging(comm, 'DEBUG', f'wspd min/max = {wspd.min()} {wspd.max()}') - - uob, vob = compute_wind_components(wdir, wspd) - - logging(comm, 'DEBUG', f'uob min/max = {uob.min()} {uob.max()}') - logging(comm, 'DEBUG', f'vob min/max = {vob.min()} {vob.max()}') - - paths = container.get_paths('variables/windSpeed', cat) - container.add('variables/windEastward', uob, paths, cat) - container.add('variables/windNorthward', vob, paths, cat) - - # Check - logging(comm, 'DEBUG', f'container list (updated): {container.list()}') - logging(comm, 'DEBUG', f'all_sub_categories {container.all_sub_categories()}') - - return container - - -def create_obs_group(input_path, mapping_path, category, env): - - comm = bufr.mpi.Comm(env["comm_name"]) - - description = _make_description(mapping_path, update=True) - - # Check the cache for the data and return it if it exists - logging(comm, 'DEBUG', f'Check if bufr.DataCache exists? {bufr.DataCache.has(input_path, mapping_path)}') - if bufr.DataCache.has(input_path, mapping_path): - container = bufr.DataCache.get(input_path, mapping_path) - logging(comm, 'INFO', f'Encode {category} from cache') - data = iodaEncoder(description).encode(container)[(category,)] - logging(comm, 'INFO', f'Mark {category} as finished in the cache') - bufr.DataCache.mark_finished(input_path, mapping_path, [category]) - logging(comm, 'INFO', f'Return the encoded data for {category}') - return data - - container = _make_obs(comm, input_path, mapping_path) - - # Gather data from all tasks into all tasks. Each task will have the complete record - logging(comm, 'INFO', f'Gather data from all tasks into all tasks') - container.all_gather(comm) - - logging(comm, 'INFO', f'Add container to cache') - # Add the container to the cache - bufr.DataCache.add(input_path, mapping_path, container.all_sub_categories(), container) - - # Encode the data - logging(comm, 'INFO', f'Encode {category}') - data = iodaEncoder(description).encode(container)[(category,)] - - logging(comm, 'INFO', f'Mark {category} as finished in the cache') - # Mark the data as finished in the cache - bufr.DataCache.mark_finished(input_path, mapping_path, [category]) - - logging(comm, 'INFO', f'Return the encoded data for {category}') - return data - - -def create_obs_file(input_path, mapping_path, output_path): - - comm = bufr.mpi.Comm("world") - container = _make_obs(comm, input_path, mapping_path) - container.gather(comm) - - description = _make_description(mapping_path, update=True) +import bufr +from bufr.obs_builder import add_main_functions +from bufr_satwnd_amv_obs_builder import SatWndAmvObsBuilder, map_path - # Encode the data - if comm.rank() == 0: - netcdfEncoder(description).encode(container, output_path) - logging(comm, 'INFO', f'Return the encoded data') +MAPPING_PATH = map_path('bufr_satwnd_amv_leogeo.yaml') -if __name__ == '__main__': +class SatWndAmvLeogeoObsBuilder(SatWndAmvObsBuilder): + def __init__(self): + super().__init__(MAPPING_PATH, log_name=os.path.basename(__file__)) - start_time = time.time() + def _get_obs_type(self, swcm, chan_freq): + obstype = swcm.copy() - bufr.mpi.App(sys.argv) - comm = bufr.mpi.Comm("world") + # Use numpy vectorized operations + obstype = np.where(swcm == 1, 255, obstype) # IRLW - # Required input arguments as positional arguments - parser = argparse.ArgumentParser(description="Convert BUFR to NetCDF using a mapping file.") - parser.add_argument('input', type=str, help='Input BUFR file') - parser.add_argument('mapping', type=str, help='BUFR2IODA Mapping File') - parser.add_argument('output', type=str, help='Output NetCDF file') + if not np.any(np.isin(obstype, [255])): + raise ValueError("Error: Unassigned ObsType found ... ") - args = parser.parse_args() - infile = args.input - mapping = args.mapping - output = args.output + return obstype - create_obs_file(infile, mapping, output) - end_time = time.time() - running_time = end_time - start_time - logging(comm, 'INFO', f'Total running time: {running_time}') +add_main_functions(SatWndAmvLeogeoObsBuilder) diff --git a/dump/mapping/bufr_satwnd_amv_leogeo_mapping.yaml b/dump/mapping/bufr_satwnd_amv_leogeo.yaml similarity index 97% rename from dump/mapping/bufr_satwnd_amv_leogeo_mapping.yaml rename to dump/mapping/bufr_satwnd_amv_leogeo.yaml index 895b475..48b2ca0 100755 --- a/dump/mapping/bufr_satwnd_amv_leogeo_mapping.yaml +++ b/dump/mapping/bufr_satwnd_amv_leogeo.yaml @@ -41,7 +41,7 @@ bufr: # Quality Information Without Forecast qualityInformationWithoutForecast: - query: '*/LGRSQ4[1]/PCCF' + query: '*/LGRSQ4[1]{1}/PCCF' # Wind Retrieval Method Information - Computation windComputationMethod: @@ -143,7 +143,7 @@ encoder: source: variables/windGeneratingApplication longName: "Wind Generating Application" - - name: "MetaData/qualityInformationWithoutForecast" + - name: "MetaData/qiWithoutForecast" source: variables/qualityInformationWithoutForecast longName: "Quality Information Without Forecast" diff --git a/dump/mapping/bufr_satwnd_amv_modis.py b/dump/mapping/bufr_satwnd_amv_modis.py index 11a5f69..d10cda1 100755 --- a/dump/mapping/bufr_satwnd_amv_modis.py +++ b/dump/mapping/bufr_satwnd_amv_modis.py @@ -1,341 +1,49 @@ #!/usr/bin/env python3 -import sys + import os -import argparse -import time import numpy as np -import bufr -from pyioda.ioda.Engines.Bufr import Encoder as iodaEncoder -from bufr.encoders.netcdf import Encoder as netcdfEncoder -from wxflow import Logger - -# Initialize Logger -# Get log level from the environment variable, default to 'INFO it not set -log_level = os.getenv('LOG_LEVEL', 'INFO') -logger = Logger('BUFR2IODA_satwnd_amv_modis.py', level=log_level, colored_log=False) - - -def logging(comm, level, message): - - if comm.rank() == 0: - # Define a dictionary to map levels to logger methods - log_methods = { - 'DEBUG': logger.debug, - 'INFO': logger.info, - 'WARNING': logger.warning, - 'ERROR': logger.error, - 'CRITICAL': logger.critical, - } - - # Get the appropriate logging method, default to 'INFO' - log_method = log_methods.get(level.upper(), logger.info) - - if log_method == logger.info and level.upper() not in log_methods: - # Log a warning if the level is invalid - logger.warning(f'log level = {level}: not a valid level --> set to INFO') - - # Call the logging method - log_method(message) - - -def _make_description(mapping_path, update=False): - description = bufr.encoders.Description(mapping_path) - - if update: - # Define the variables to be added in a list of dictionaries - variables = [ - { - 'name': 'ObsType/windEastward', - 'source': 'variables/obstype_uwind', - 'units': '1', - 'longName': 'Observation Type based on Satellite-derived Wind Computation Method and Spectral Band', - }, - { - 'name': 'ObsType/windNorthward', - 'source': 'variables/obstype_vwind', - 'units': '1', - 'longName': 'Observation Type based on Satellite-derived Wind Computation Method and Spectral Band', - }, - { - 'name': 'ObsValue/windEastward', - 'source': 'variables/windEastward', - 'units': 'm s-1', - 'longName': 'Eastward Wind Component', - }, - { - 'name': 'ObsValue/windNorthward', - 'source': 'variables/windNorthward', - 'units': 'm s-1', - 'longName': 'Northward Wind Component', - }, - # MetaData/windGeneratingApplication will be inferred from variables/generatingApplication - # following a search for the proper variables/generatingApplication column - { - 'name': 'MetaData/windGeneratingApplication', - 'source': 'variables/windGeneratingApplication', - 'units': '1', - 'longName': 'Wind Generating Application', - }, - # MetaData/qualityInformationWithoutForecast will be inferred from variables/qualityInformation - # following a search for the proper variables/generatingApplication column - { - 'name': 'MetaData/qualityInformationWithoutForecast', - 'source': 'variables/qualityInformationWithoutForecast', - 'units': 'percent', - 'longName': 'Quality Information Without Forecast', - } - ] - - # Loop through each variable and add it to the description - for var in variables: - description.add_variable( - name=var['name'], - source=var['source'], - units=var['units'], - longName=var['longName'] - ) - - return description - - -def compute_wind_components(wdir, wspd): - """ - Compute the U and V wind components from wind direction and wind speed. - - Parameters: - wdir (array-like): Wind direction in degrees (meteorological convention: 0° = North, 90° = East). - wspd (array-like): Wind speed. - - Returns: - tuple: U and V wind components as numpy arrays with dtype float32. - """ - wdir_rad = np.radians(wdir) # Convert degrees to radians - u = -wspd * np.sin(wdir_rad) - v = -wspd * np.cos(wdir_rad) - - return u.astype(np.float32), v.astype(np.float32) - - -def _get_quality_information_and_generating_application(comm, gnap2D, pccf2D): - # For MODIS Terra/Aqua data, qi w/o forecast (qifn) is packaged in same - # vector of qi with ga = 1 (EUMETSAT QI without forecast). Must - # conduct a search and extract the correct vector for gnap and qi - findQI = 1 - # 1. Find dimension-sizes of ga and qi (should be the same!) - gDim1, gDim2 = np.shape(gnap2D) - qDim1, qDim2 = np.shape(pccf2D) - logging(comm, 'INFO', 'Generating Application and Quality Information SEARCH') - logging(comm, 'DEBUG', f'Dimension size of GNAP ({gDim1},{gDim2})') - logging(comm, 'DEBUG', f'Dimension size of PCCF ({qDim1},{qDim2})') - # 2. Initialize gnap and qifn as None, and search for dimension of - # ga with values of 5. If the same column exists for qi, assign - # gnap to ga[:,i] and qifn to qi[:,i], else raise warning that no - # appropriate GNAP/PCCF combination was found - gnap = None - qifn = None - for i in range(gDim2): - if np.unique(gnap2D[:, i].squeeze()) == findQI: - if i <= qDim2: - logging(comm, 'INFO', f'GNAP/PCCF found for column {i}') - gnap = gnap2D[:, i].squeeze() - qifn = pccf2D[:, i].squeeze() - else: - logging(comm, 'INFO', f'ERROR: GNAP column {i} outside of PCCF dimension {qDim2}') - if (gnap is None) & (qifn is None): - raise ValueError(f'GNAP == {findQI} NOT FOUND OR OUT OF PCCF DIMENSION-RANGE, WILL FAIL!') - # If EE is needed, key search on np.unique(gnap2D[:,i].squeeze()) == 4 instead - # NOTE: Make sure to return np.float32 or np.int32 types as appropriate!!! - return gnap.astype(np.int32), qifn.astype(np.int32) - - -def _get_obs_type(swcm): - """ - Determine the observation type based on `swcm` and `chanfreq`. - - Parameters: - swcm (array-like): Switch mode values. - chanfreq (array-like): Channel frequency values (Hz). - - Returns: - numpy.ndarray: Observation type array. - - Raises: - ValueError: If any `obstype` is unassigned. - """ - - obstype = swcm.copy() - - # Use numpy vectorized operations - obstype = np.where(swcm >= 4, 259, obstype) # WVCA/DL - obstype = np.where(swcm == 3, 258, obstype) # WVCT - obstype = np.where(swcm == 1, 257, obstype) # IRLW - - if not np.any(np.isin(obstype, [257, 258, 259])): - raise ValueError("Error: Unassigned ObsType found ... ") - - return obstype - - -def _make_obs(comm, input_path, mapping_path): - - # Get container from mapping file first - logging(comm, 'INFO', 'Get container from bufr') - container = bufr.Parser(input_path, mapping_path).parse(comm) - logging(comm, 'DEBUG', f'container list (original): {container.list()}') - logging(comm, 'DEBUG', f'all_sub_categories = {container.all_sub_categories()}') - logging(comm, 'DEBUG', f'category map = {container.get_category_map()}') - - # Add new/derived data into container - for cat in container.all_sub_categories(): - - logging(comm, 'DEBUG', f'category = {cat}') - - satid = container.get('variables/satelliteId', cat) - if satid.size == 0: - logging(comm, 'WARNING', f'category {cat[0]} does not exist in input file') - paths = container.get_paths('variables/windComputationMethod', cat) - obstype = container.get('variables/windComputationMethod', cat) - container.add('variables/obstype_uwind', obstype, paths, cat) - container.add('variables/obstype_vwind', obstype, paths, cat) - - paths = container.get_paths('variables/windSpeed', cat) - wob = container.get('variables/windSpeed', cat) - container.add('variables/windEastward', wob, paths, cat) - container.add('variables/windNorthward', wob, paths, cat) - - paths = container.get_paths('variables/windComputationMethod', cat) - dummy = container.get('variables/windSpeed', cat) - container.add('variables/windGeneratingApplication', dummy, paths, cat) - container.add('variables/qualityInformationWithoutForecast', dummy, paths, cat) - - else: - # Add new variables: ObsType/windEastward & ObsType/windNorthward - swcm = container.get('variables/windComputationMethod', cat) - chanfreq = container.get('variables/sensorCentralFrequency', cat) - - logging(comm, 'DEBUG', f'swcm min/max = {swcm.min()} {swcm.max()}') - logging(comm, 'DEBUG', f'chanfreq min/max = {chanfreq.min()} {chanfreq.max()}') - - obstype = _get_obs_type(swcm) - - logging(comm, 'DEBUG', f'obstype = {obstype}') - logging(comm, 'DEBUG', f'obstype min/max = {obstype.min()} {obstype.max()}') - - paths = container.get_paths('variables/windComputationMethod', cat) - container.add('variables/obstype_uwind', obstype, paths, cat) - container.add('variables/obstype_vwind', obstype, paths, cat) - - # Add new variabls: ObsValue/windEastward & ObsValue/windNorthward - wdir = container.get('variables/windDirection', cat) - wspd = container.get('variables/windSpeed', cat) - - logging(comm, 'DEBUG', f'wdir min/max = {wdir.min()} {wdir.max()}') - logging(comm, 'DEBUG', f'wspd min/max = {wspd.min()} {wspd.max()}') - - uob, vob = compute_wind_components(wdir, wspd) - - logging(comm, 'DEBUG', f'uob min/max = {uob.min()} {uob.max()}') - logging(comm, 'DEBUG', f'vob min/max = {vob.min()} {vob.max()}') - - paths = container.get_paths('variables/windSpeed', cat) - container.add('variables/windEastward', uob, paths, cat) - container.add('variables/windNorthward', vob, paths, cat) - - # Add new variables: MetaData/windGeneratingApplication and qualityInformationWithoutForecast - gnap2D = container.get('variables/generatingApplication', cat) - pccf2D = container.get('variables/qualityInformation', cat) - - gnap, qifn = _get_quality_information_and_generating_application(comm, gnap2D, pccf2D) - - logging(comm, 'DEBUG', f'gnap min/max = {gnap.min()} {gnap.max()}') - logging(comm, 'DEBUG', f'qifn min/max = {qifn.min()} {qifn.max()}') - - paths = container.get_paths('variables/windComputationMethod', cat) - container.add('variables/windGeneratingApplication', gnap, paths, cat) - container.add('variables/qualityInformationWithoutForecast', qifn, paths, cat) - - # Check - logging(comm, 'DEBUG', f'container list (updated): {container.list()}') - logging(comm, 'DEBUG', f'all_sub_categories {container.all_sub_categories()}') - - return container - - -def create_obs_group(input_path, mapping_path, category, env): - - comm = bufr.mpi.Comm(env["comm_name"]) - - description = _make_description(mapping_path, update=True) - - # Check the cache for the data and return it if it exists - logging(comm, 'DEBUG', f'Check if bufr.DataCache exists? {bufr.DataCache.has(input_path, mapping_path)}') - if bufr.DataCache.has(input_path, mapping_path): - container = bufr.DataCache.get(input_path, mapping_path) - logging(comm, 'INFO', f'Encode {category} from cache') - data = iodaEncoder(description).encode(container)[(category,)] - logging(comm, 'INFO', f'Mark {category} as finished in the cache') - bufr.DataCache.mark_finished(input_path, mapping_path, [category]) - logging(comm, 'INFO', f'Return the encoded data for {category}') - return data - - container = _make_obs(comm, input_path, mapping_path) - - # Gather data from all tasks into all tasks. Each task will have the complete record - logging(comm, 'INFO', f'Gather data from all tasks into all tasks') - container.all_gather(comm) - - logging(comm, 'INFO', f'Add container to cache') - # Add the container to the cache - bufr.DataCache.add(input_path, mapping_path, container.all_sub_categories(), container) - - # Encode the data - logging(comm, 'INFO', f'Encode {category}') - data = iodaEncoder(description).encode(container)[(category,)] - - logging(comm, 'INFO', f'Mark {category} as finished in the cache') - # Mark the data as finished in the cache - bufr.DataCache.mark_finished(input_path, mapping_path, [category]) +import bufr +from bufr.obs_builder import add_main_functions +from bufr_satwnd_amv_obs_builder import SatWndAmvObsBuilder, map_path - logging(comm, 'INFO', f'Return the encoded data for {category}') - return data +MAPPING_PATH = map_path('bufr_satwnd_amv_modis.yaml') +FIND_QI = 1 -def create_obs_file(input_path, mapping_path, output_path): - comm = bufr.mpi.Comm("world") - container = _make_obs(comm, input_path, mapping_path) - container.gather(comm) +class SatWndAmvModisObsBuilder(SatWndAmvObsBuilder): + def __init__(self): + super().__init__(MAPPING_PATH, log_name=os.path.basename(__file__)) - description = _make_description(mapping_path, update=True) + def make_obs(self, comm, input_path): + container = super().make_obs(comm, input_path) - # Encode the data - if comm.rank() == 0: - netcdfEncoder(description).encode(container, output_path) + # Add new/derived data into container + for cat in container.all_sub_categories(): + self._add_quality_info_and_gen_app(FIND_QI, container, cat) - logging(comm, 'INFO', f'Return the encoded data') + return container + def _make_description(self): + description = super()._make_description() + self._add_quality_info_and_gen_app_descriptions(description) -if __name__ == '__main__': + return description - start_time = time.time() + def _get_obs_type(self, swcm, chan_freq): + obstype = swcm.copy() - bufr.mpi.App(sys.argv) - comm = bufr.mpi.Comm("world") + # Use numpy vectorized operations + obstype = np.where(swcm >= 4, 259, obstype) # WVCA/DL + obstype = np.where(swcm == 3, 258, obstype) # WVCT + obstype = np.where(swcm == 1, 257, obstype) # IRLW - # Required input arguments as positional arguments - parser = argparse.ArgumentParser(description="Convert BUFR to NetCDF using a mapping file.") - parser.add_argument('input', type=str, help='Input BUFR file') - parser.add_argument('mapping', type=str, help='BUFR2IODA Mapping File') - parser.add_argument('output', type=str, help='Output NetCDF file') + if not np.any(np.isin(obstype, [257, 258, 259])): + raise ValueError("Error: Unassigned ObsType found ... ") - args = parser.parse_args() - infile = args.input - mapping = args.mapping - output = args.output + return obstype - create_obs_file(infile, mapping, output) - end_time = time.time() - running_time = end_time - start_time - logging(comm, 'INFO', f'Total running time: {running_time}') +# Add main functions create_obs_file and create_obs_group +add_main_functions(SatWndAmvModisObsBuilder) diff --git a/dump/mapping/bufr_satwnd_amv_modis_mapping.yaml b/dump/mapping/bufr_satwnd_amv_modis.yaml similarity index 100% rename from dump/mapping/bufr_satwnd_amv_modis_mapping.yaml rename to dump/mapping/bufr_satwnd_amv_modis.yaml diff --git a/dump/mapping/bufr_satwnd_amv_obs_builder.py b/dump/mapping/bufr_satwnd_amv_obs_builder.py new file mode 100644 index 0000000..6b59909 --- /dev/null +++ b/dump/mapping/bufr_satwnd_amv_obs_builder.py @@ -0,0 +1,228 @@ +#!/usr/bin/env python3 + +import os +import numpy as np + +import bufr +from bufr.obs_builder import ObsBuilder + + +def map_path(map_file_name): + script_dir = os.path.dirname(os.path.abspath(__file__)) + return os.path.join(script_dir, map_file_name) + + +class SatWndAmvObsBuilder(ObsBuilder): + def __init__(self, mapping_path, log_name=os.path.basename(__file__)): + super().__init__(mapping_path, log_name=log_name) + + def make_obs(self, comm, input_path): + # Get container from mapping file first + self.log.info('Get container from bufr') + container = super().make_obs(comm, input_path) + + self.log.debug(f'container list (original): {container.list()}') + self.log.debug(f'all_sub_categories = {container.all_sub_categories()}') + self.log.debug(f'category map = {container.get_category_map()}') + + # Add new/derived data into container + for cat in container.all_sub_categories(): + self.log.debug(f'category = {cat}') + + satId = container.get('satelliteId', cat) + if not np.any(satId): + self.log.warning(f'category {cat[0]} does not exist in input file') + + self._add_wind_obs(container, cat) + + # Check + self.log.debug(f'container list (updated): {container.list()}') + self.log.debug('all_sub_categories {container.all_sub_categories()}') + + return container + + # Provide defualt implementations for methods from the ObsBuilder class + def _make_description(self): + description = super()._make_description() + self._add_wind_descriptions(description) + + return description + + # Methods that are used to extend the export description + def _add_wind_descriptions(self, description): + description.add_variables([ + { + 'name': 'ObsType/windEastward', + 'source': 'obstype_uwind', + 'units': '1', + 'longName': 'Observation Type based on Satellite-derived Wind Computation Method and Spectral Band', + }, + { + 'name': 'ObsType/windNorthward', + 'source': 'obstype_vwind', + 'units': '1', + 'longName': 'Observation Type based on Satellite-derived Wind Computation Method and Spectral Band', + }, + { + 'name': 'ObsValue/windEastward', + 'source': 'windEastward', + 'units': 'm s-1', + 'longName': 'Eastward Wind Component', + }, + { + 'name': 'ObsValue/windNorthward', + 'source': 'windNorthward', + 'units': 'm s-1', + 'longName': 'Northward Wind Component', + }]) + + def _add_quality_info_and_gen_app_descriptions(self, description): + description.add_variables([ + { + 'name': 'MetaData/windGeneratingApplication', + 'source': 'windGeneratingApplication', + 'units': '1', + 'longName': 'Wind Generating Application', + }, + { + 'name': 'MetaData/qiWithoutForecast', + 'source': 'qualityInformationWithoutForecast', + 'units': '1', + 'longName': 'Quality Information Without Forecast', + }]) + + # Methods that are used to extend the obs data container + def _add_wind_obs(self, container, cat): + # Add new variables: ObsType/windEastward & ObsType/windNorthward + swcm = container.get('windComputationMethod', cat) + chanfreq = container.get('sensorCentralFrequency', cat) + + if swcm.size == 0: + self.log.warning(f'category {cat[0]} does not exist in input file') + paths = container.get_paths('variables/windComputationMethod', cat) + obstype = container.get('variables/windComputationMethod', cat) + container.add('variables/obstype_uwind', obstype, paths, cat) + container.add('variables/obstype_vwind', obstype, paths, cat) + + paths = container.get_paths('variables/windSpeed', cat) + wob = container.get('variables/windSpeed', cat) + container.add('variables/windEastward', wob, paths, cat) + container.add('variables/windNorthward', wob, paths, cat) + return + + # self.log.debug(f'swcm min/max = {swcm.min()} {swcm.max()}') + self.log.debug('chanfreq min/max = {chanfreq.min()} {chanfreq.max()}') + + obstype = self._get_obs_type(swcm, chanfreq) + + self.log.debug(f'obstype = {obstype}') + self.log.debug(f'obstype min/max = {obstype.min()} {obstype.max()}') + + paths = container.get_paths('windComputationMethod', cat) + container.add('obstype_uwind', obstype, paths, cat) + container.add('obstype_vwind', obstype, paths, cat) + + # Add new variables: ObsValue/windEastward & ObsValue/windNorthward + wdir = container.get('windDirection', cat) + wspd = container.get('windSpeed', cat) + + self.log.debug(f'wdir min/max = {wdir.min()} {wdir.max()}') + self.log.debug(f'wspd min/max = {wspd.min()} {wspd.max()}') + + uob, vob = self._compute_wind_components(wdir, wspd) + + self.log.debug(f'uob min/max = {uob.min()} {uob.max()}') + self.log.debug(f'vob min/max = {vob.min()} {vob.max()}') + + paths = container.get_paths('windSpeed', cat) + container.add('windEastward', uob, paths, cat) + container.add('windNorthward', vob, paths, cat) + + def _add_quality_info_and_gen_app(self, findQi, container, cat): + # Add new variables: MetaData/windGeneratingApplication and qiWithoutForecast + gnap2D = container.get('generatingApplication', cat) + pccf2D = container.get('qualityInformation', cat) + satId = container.get('satelliteId', cat) + + if not satId.size: + paths = container.get_paths('windComputationMethod', cat) + dummy = container.get('windSpeed', cat) + container.add('windGeneratingApplication', dummy, paths, cat) + container.add('qualityInformationWithoutForecast', dummy, paths, cat) + return + + gnap, qifn = self._get_quality_info_and_gen_app(findQi, gnap2D, pccf2D) + + self.log.debug(f'gnap min/max = {gnap.min()} {gnap.max()}') + self.log.debug(f'qifn min/max = {qifn.min()} {qifn.max()}') + + paths = container.get_paths('windComputationMethod', cat) + container.add('windGeneratingApplication', gnap, paths, cat) + container.add('qualityInformationWithoutForecast', qifn, paths, cat) + + def _get_obs_type(self, swcm, chan_freq=0): + """ + Determine the observation type based on `swcm` and `chanfreq`. + + Parameters: + swcm (array-like): Switch mode values. + chanfreq (array-like): Channel frequency values (Hz). + + Returns: + numpy.ndarray: Observation type array. + + Raises: + ValueError: If any `obstype` is unassigned. + """ + + raise NotImplementedError('Method _get_obs_type must be implemented in derived classes') + + # Private methods + def _compute_wind_components(self, wdir, wspd): + """ + Compute the U and V wind components from wind direction and wind speed. + + Parameters: + wdir (array-like): Wind direction in degrees (meteorological convention: 0° = North, 90° = East). + wspd (array-like): Wind speed. + + Returns: + tuple: U and V wind components as numpy arrays with dtype float32. + """ + wdir_rad = np.radians(wdir) # Convert degrees to radians + u = -wspd * np.sin(wdir_rad) + v = -wspd * np.cos(wdir_rad) + + return u.astype(np.float32), v.astype(np.float32) + + def _get_quality_info_and_gen_app(self, findQi, gnap2D, pccf2D): + # For NOAA VIIRS data, qi w/o forecast (qifn) is packaged in same + # vector of qi with ga = 5 (EUMETSAT QI without forecast). Must + # conduct a search and extract the correct vector for gnap and qi + + # 1. Find dimension-sizes of ga and qi (should be the same!) + gDim1, gDim2 = np.shape(gnap2D) + qDim1, qDim2 = np.shape(pccf2D) + self.log.info('Generating Application and Quality Information SEARCH') + self.log.debug(f'Dimension size of GNAP ({gDim1},{gDim2})') + self.log.debug(f'Dimension size of PCCF ({qDim1},{qDim2})') + + # 2. Initialize gnap and qifn as None, and search for dimension of + # ga with values of 5. If the same column exists for qi, assign + # gnap to ga[:,i] and qifn to qi[:,i], else raise warning that no + # appropriate GNAP/PCCF combination was found + gnap = None + qifn = None + for i in range(gDim2): + if np.unique(gnap2D[:, i].squeeze()) == findQi: + if i <= qDim2: + self.log.info(f'GNAP/PCCF found for column {i}') + gnap = gnap2D[:, i].squeeze() + qifn = pccf2D[:, i].squeeze() + else: + self.log.info(f'ERROR: GNAP column {i} outside of PCCF dimension {qDim2}') + if (gnap is None) & (qifn is None): + raise ValueError(f'GNAP == {findQI} NOT FOUND OR OUT OF PCCF DIMENSION-RANGE, WILL FAIL!') + # If EE is needed, key search on np.unique(gnap2D[:,i].squeeze()) == 7 instead + # NOTE: Make sure to return np.float32 or np.int32 types as appropriate!!! + return gnap.astype(np.int32), qifn.astype(np.int32) diff --git a/dump/mapping/bufr_satwnd_amv_seviri.py b/dump/mapping/bufr_satwnd_amv_seviri.py index d673900..422bce3 100755 --- a/dump/mapping/bufr_satwnd_amv_seviri.py +++ b/dump/mapping/bufr_satwnd_amv_seviri.py @@ -1,276 +1,34 @@ #!/usr/bin/env python3 -import sys + import os -import argparse -import time import numpy as np -import bufr -from pyioda.ioda.Engines.Bufr import Encoder as iodaEncoder -from bufr.encoders.netcdf import Encoder as netcdfEncoder -from wxflow import Logger - -# Initialize Logger -# Get log level from the environment variable, default to 'INFO it not set -log_level = os.getenv('LOG_LEVEL', 'INFO') -logger = Logger('BUFR2IODA_satwnd_amv_seviri.py', level=log_level, colored_log=False) - - -def logging(comm, level, message): - - if comm.rank() == 0: - # Define a dictionary to map levels to logger methods - log_methods = { - 'DEBUG': logger.debug, - 'INFO': logger.info, - 'WARNING': logger.warning, - 'ERROR': logger.error, - 'CRITICAL': logger.critical, - } - - # Get the appropriate logging method, default to 'INFO' - log_method = log_methods.get(level.upper(), logger.info) - - if log_method == logger.info and level.upper() not in log_methods: - # Log a warning if the level is invalid - logger.warning(f'log level = {level}: not a valid level --> set to INFO') - - # Call the logging method - log_method(message) - - -def _make_description(mapping_path, update=False): - description = bufr.encoders.Description(mapping_path) - - if update: - # Define the variables to be added in a list of dictionaries - variables = [ - { - 'name': 'ObsType/windEastward', - 'source': 'variables/obstype_uwind', - 'units': '1', - 'longName': 'Observation Type based on Satellite-derived Wind Computation Method and Spectral Band', - }, - { - 'name': 'ObsType/windNorthward', - 'source': 'variables/obstype_vwind', - 'units': '1', - 'longName': 'Observation Type based on Satellite-derived Wind Computation Method and Spectral Band', - }, - { - 'name': 'ObsValue/windEastward', - 'source': 'variables/windEastward', - 'units': 'm s-1', - 'longName': 'Eastward Wind Component', - }, - { - 'name': 'ObsValue/windNorthward', - 'source': 'variables/windNorthward', - 'units': 'm s-1', - 'longName': 'Northward Wind Component', - } - ] - - # Loop through each variable and add it to the description - for var in variables: - description.add_variable( - name=var['name'], - source=var['source'], - units=var['units'], - longName=var['longName'] - ) - - return description - - -def compute_wind_components(wdir, wspd): - """ - Compute the U and V wind components from wind direction and wind speed. - - Parameters: - wdir (array-like): Wind direction in degrees (meteorological convention: 0° = North, 90° = East). - wspd (array-like): Wind speed. - - Returns: - tuple: U and V wind components as numpy arrays with dtype float32. - """ - wdir_rad = np.radians(wdir) # Convert degrees to radians - u = -wspd * np.sin(wdir_rad) - v = -wspd * np.cos(wdir_rad) - - return u.astype(np.float32), v.astype(np.float32) - - -def _get_obs_type(swcm): - """ - Determine the observation type based on `swcm` and `chanfreq`. - - Parameters: - swcm (array-like): Switch mode values. - chanfreq (array-like): Channel frequency values (Hz). - - Returns: - numpy.ndarray: Observation type array. - - Raises: - ValueError: If any `obstype` is unassigned. - """ - - obstype = swcm.copy() - - # Use numpy vectorized operations - obstype = np.where(swcm == 5, 254, obstype) # WVCA/DL - obstype = np.where(swcm == 3, 254, obstype) # WVCT - obstype = np.where(swcm == 2, 243, obstype) # VIS - obstype = np.where(swcm == 1, 253, obstype) # IRLW - - if not np.any(np.isin(obstype, [243, 253, 254])): - raise ValueError("Error: Unassigned ObsType found ... ") - - return obstype - - -def _make_obs(comm, input_path, mapping_path): - - # Get container from mapping file first - logging(comm, 'INFO', 'Get container from bufr') - container = bufr.Parser(input_path, mapping_path).parse(comm) - - logging(comm, 'DEBUG', f'container list (original): {container.list()}') - logging(comm, 'DEBUG', f'all_sub_categories = {container.all_sub_categories()}') - logging(comm, 'DEBUG', f'category map = {container.get_category_map()}') - # Add new/derived data into container - for cat in container.all_sub_categories(): - - logging(comm, 'DEBUG', f'category = {cat}') - - satid = container.get('variables/satelliteId', cat) - if satid.size == 0: - logging(comm, 'WARNING', f'category {cat[0]} does not exist in input file') - paths = container.get_paths('variables/windComputationMethod', cat) - obstype = container.get('variables/windComputationMethod', cat) - container.add('variables/obstype_uwind', obstype, paths, cat) - container.add('variables/obstype_vwind', obstype, paths, cat) - - paths = container.get_paths('variables/windSpeed', cat) - wob = container.get('variables/windSpeed', cat) - container.add('variables/windEastward', wob, paths, cat) - container.add('variables/windNorthward', wob, paths, cat) - - else: - # Add new variables: ObsType/windEastward & ObsType/windNorthward - swcm = container.get('variables/windComputationMethod', cat) - chanfreq = container.get('variables/sensorCentralFrequency', cat) - - logging(comm, 'DEBUG', f'swcm min/max = {swcm.min()} {swcm.max()}') - logging(comm, 'DEBUG', f'chanfreq min/max = {chanfreq.min()} {chanfreq.max()}') - - obstype = _get_obs_type(swcm) - - logging(comm, 'DEBUG', f'obstype = {obstype}') - logging(comm, 'DEBUG', f'obstype min/max = {obstype.min()} {obstype.max()}') - - paths = container.get_paths('variables/windComputationMethod', cat) - container.add('variables/obstype_uwind', obstype, paths, cat) - container.add('variables/obstype_vwind', obstype, paths, cat) - - # Add new variables: ObsValue/windEastward & ObsValue/windNorthward - wdir = container.get('variables/windDirection', cat) - wspd = container.get('variables/windSpeed', cat) - - logging(comm, 'DEBUG', f'wdir min/max = {wdir.min()} {wdir.max()}') - logging(comm, 'DEBUG', f'wspd min/max = {wspd.min()} {wspd.max()}') - - uob, vob = compute_wind_components(wdir, wspd) - - logging(comm, 'DEBUG', f'uob min/max = {uob.min()} {uob.max()}') - logging(comm, 'DEBUG', f'vob min/max = {vob.min()} {vob.max()}') - - paths = container.get_paths('variables/windSpeed', cat) - container.add('variables/windEastward', uob, paths, cat) - container.add('variables/windNorthward', vob, paths, cat) - - # Check - logging(comm, 'DEBUG', f'container list (updated): {container.list()}') - logging(comm, 'DEBUG', f'all_sub_categories {container.all_sub_categories()}') - - return container - - -def create_obs_group(input_path, mapping_path, category, env): - - comm = bufr.mpi.Comm(env["comm_name"]) - - description = _make_description(mapping_path, update=True) - - # Check the cache for the data and return it if it exists - logging(comm, 'DEBUG', f'Check if bufr.DataCache exists? {bufr.DataCache.has(input_path, mapping_path)}') - if bufr.DataCache.has(input_path, mapping_path): - container = bufr.DataCache.get(input_path, mapping_path) - logging(comm, 'INFO', f'Encode {category} from cache') - data = iodaEncoder(description).encode(container)[(category,)] - logging(comm, 'INFO', f'Mark {category} as finished in the cache') - bufr.DataCache.mark_finished(input_path, mapping_path, [category]) - logging(comm, 'INFO', f'Return the encoded data for {category}') - return data - - container = _make_obs(comm, input_path, mapping_path) - - # Gather data from all tasks into all tasks. Each task will have the complete record - logging(comm, 'INFO', f'Gather data from all tasks into all tasks') - container.all_gather(comm) - - logging(comm, 'INFO', f'Add container to cache') - # Add the container to the cache - bufr.DataCache.add(input_path, mapping_path, container.all_sub_categories(), container) - - # Encode the data - logging(comm, 'INFO', f'Encode {category}') - data = iodaEncoder(description).encode(container)[(category,)] - - logging(comm, 'INFO', f'Mark {category} as finished in the cache') - # Mark the data as finished in the cache - bufr.DataCache.mark_finished(input_path, mapping_path, [category]) - - logging(comm, 'INFO', f'Return the encoded data for {category}') - return data - - -def create_obs_file(input_path, mapping_path, output_path): - - comm = bufr.mpi.Comm("world") - container = _make_obs(comm, input_path, mapping_path) - container.gather(comm) - - description = _make_description(mapping_path, update=True) +import bufr +from bufr.obs_builder import add_main_functions +from bufr_satwnd_amv_obs_builder import SatWndAmvObsBuilder, map_path - # Encode the data - if comm.rank() == 0: - netcdfEncoder(description).encode(container, output_path) - logging(comm, 'INFO', f'Return the encoded data') +MAPPING_PATH = map_path('bufr_satwnd_amv_seviri.yaml') -if __name__ == '__main__': +class SatWndAmvSeviriObsBuilder(SatWndAmvObsBuilder): + def __init__(self): + super().__init__(MAPPING_PATH, log_name=os.path.basename(__file__)) - start_time = time.time() + def _get_obs_type(self, swcm, chan_freq): + obstype = swcm.copy() - bufr.mpi.App(sys.argv) - comm = bufr.mpi.Comm("world") + # Use numpy vectorized operations + obstype = np.where(swcm == 5, 254, obstype) # WVCA/DL + obstype = np.where(swcm == 3, 254, obstype) # WVCT + obstype = np.where(swcm == 2, 243, obstype) # VIS + obstype = np.where(swcm == 1, 253, obstype) # IRLW - # Required input arguments as positional arguments - parser = argparse.ArgumentParser(description="Convert BUFR to NetCDF using a mapping file.") - parser.add_argument('input', type=str, help='Input BUFR file') - parser.add_argument('mapping', type=str, help='BUFR2IODA Mapping File') - parser.add_argument('output', type=str, help='Output NetCDF file') + if not np.any(np.isin(obstype, [243, 253, 254])): + raise ValueError("Error: Unassigned ObsType found ... ") - args = parser.parse_args() - infile = args.input - mapping = args.mapping - output = args.output + return obstype - create_obs_file(infile, mapping, output) - end_time = time.time() - running_time = end_time - start_time - logging(comm, 'INFO', f'Total running time: {running_time}') +# Add main functions create_obs_file and create_obs_group +add_main_functions(SatWndAmvSeviriObsBuilder) diff --git a/dump/mapping/bufr_satwnd_amv_seviri_mapping.yaml b/dump/mapping/bufr_satwnd_amv_seviri.yaml similarity index 98% rename from dump/mapping/bufr_satwnd_amv_seviri_mapping.yaml rename to dump/mapping/bufr_satwnd_amv_seviri.yaml index 89d9c86..8c78a67 100755 --- a/dump/mapping/bufr_satwnd_amv_seviri_mapping.yaml +++ b/dump/mapping/bufr_satwnd_amv_seviri.yaml @@ -148,7 +148,7 @@ encoder: source: variables/windGeneratingApplication longName: "Wind Generating Application" - - name: "MetaData/qualityInformationWithoutForecast" + - name: "MetaData/qiWithoutForecast" source: variables/qualityInformationWithoutForecast longName: "Quality Information Without Forecast" diff --git a/dump/mapping/bufr_satwnd_amv_viirs.py b/dump/mapping/bufr_satwnd_amv_viirs.py index 7d0a2c9..1f1d7ea 100755 --- a/dump/mapping/bufr_satwnd_amv_viirs.py +++ b/dump/mapping/bufr_satwnd_amv_viirs.py @@ -1,339 +1,47 @@ #!/usr/bin/env python3 -import sys + import os -import argparse -import time import numpy as np -import bufr -from pyioda.ioda.Engines.Bufr import Encoder as iodaEncoder -from bufr.encoders.netcdf import Encoder as netcdfEncoder -from wxflow import Logger - -# Initialize Logger -# Get log level from the environment variable, default to 'INFO it not set -log_level = os.getenv('LOG_LEVEL', 'INFO') -logger = Logger('BUFR2IODA_satwnd_amv_viirs.py', level=log_level, colored_log=False) - - -def logging(comm, level, message): - - if comm.rank() == 0: - # Define a dictionary to map levels to logger methods - log_methods = { - 'DEBUG': logger.debug, - 'INFO': logger.info, - 'WARNING': logger.warning, - 'ERROR': logger.error, - 'CRITICAL': logger.critical, - } - - # Get the appropriate logging method, default to 'INFO' - log_method = log_methods.get(level.upper(), logger.info) - - if log_method == logger.info and level.upper() not in log_methods: - # Log a warning if the level is invalid - logger.warning(f'log level = {level}: not a valid level --> set to INFO') - - # Call the logging method - log_method(message) - - -def _make_description(mapping_path, update=False): - description = bufr.encoders.Description(mapping_path) - - if update: - # Define the variables to be added in a list of dictionaries - variables = [ - { - 'name': 'ObsType/windEastward', - 'source': 'variables/obstype_uwind', - 'units': '1', - 'longName': 'Observation Type based on Satellite-derived Wind Computation Method and Spectral Band', - }, - { - 'name': 'ObsType/windNorthward', - 'source': 'variables/obstype_vwind', - 'units': '1', - 'longName': 'Observation Type based on Satellite-derived Wind Computation Method and Spectral Band', - }, - { - 'name': 'ObsValue/windEastward', - 'source': 'variables/windEastward', - 'units': 'm s-1', - 'longName': 'Eastward Wind Component', - }, - { - 'name': 'ObsValue/windNorthward', - 'source': 'variables/windNorthward', - 'units': 'm s-1', - 'longName': 'Northward Wind Component', - }, - # MetaData/windGeneratingApplication will be inferred from variables/generatingApplication - # following a search for the proper variables/generatingApplication column - { - 'name': 'MetaData/windGeneratingApplication', - 'source': 'variables/windGeneratingApplication', - 'units': '1', - 'longName': 'Wind Generating Application', - }, - # MetaData/qualityInformationWithoutForecast will be inferred from variables/qualityInformation - # following a search for the proper variables/generatingApplication column - { - 'name': 'MetaData/qualityInformationWithoutForecast', - 'source': 'variables/qualityInformationWithoutForecast', - 'units': 'percent', - 'longName': 'Quality Information Without Forecast', - } - ] - - # Loop through each variable and add it to the description - for var in variables: - description.add_variable( - name=var['name'], - source=var['source'], - units=var['units'], - longName=var['longName'] - ) - - return description - - -def compute_wind_components(wdir, wspd): - """ - Compute the U and V wind components from wind direction and wind speed. - - Parameters: - wdir (array-like): Wind direction in degrees (meteorological convention: 0° = North, 90° = East). - wspd (array-like): Wind speed. - - Returns: - tuple: U and V wind components as numpy arrays with dtype float32. - """ - wdir_rad = np.radians(wdir) # Convert degrees to radians - u = -wspd * np.sin(wdir_rad) - v = -wspd * np.cos(wdir_rad) - - return u.astype(np.float32), v.astype(np.float32) - - -def _get_quality_information_and_generating_application(comm, gnap2D, pccf2D): - # For NOAA VIIRS data, qi w/o forecast (qifn) is packaged in same - # vector of qi with ga = 5 (EUMETSAT QI without forecast). Must - # conduct a search and extract the correct vector for gnap and qi - findQI = 5 - # 1. Find dimension-sizes of ga and qi (should be the same!) - gDim1, gDim2 = np.shape(gnap2D) - qDim1, qDim2 = np.shape(pccf2D) - logging(comm, 'INFO', 'Generating Application and Quality Information SEARCH') - logging(comm, 'DEBUG', f'Dimension size of GNAP ({gDim1},{gDim2})') - logging(comm, 'DEBUG', f'Dimension size of PCCF ({qDim1},{qDim2})') - # 2. Initialize gnap and qifn as None, and search for dimension of - # ga with values of 5. If the same column exists for qi, assign - # gnap to ga[:,i] and qifn to qi[:,i], else raise warning that no - # appropriate GNAP/PCCF combination was found - gnap = None - qifn = None - for i in range(gDim2): - if np.unique(gnap2D[:, i].squeeze()) == findQI: - if i <= qDim2: - logging(comm, 'INFO', f'GNAP/PCCF found for column {i}') - gnap = gnap2D[:, i].squeeze() - qifn = pccf2D[:, i].squeeze() - else: - logging(comm, 'INFO', f'ERROR: GNAP column {i} outside of PCCF dimension {qDim2}') - if (gnap is None) & (qifn is None): - raise ValueError(f'GNAP == {findQI} NOT FOUND OR OUT OF PCCF DIMENSION-RANGE, WILL FAIL!') - # If EE is needed, key search on np.unique(gnap2D[:,i].squeeze()) == 7 instead - # NOTE: Make sure to return np.float32 or np.int32 types as appropriate!!! - return gnap.astype(np.int32), qifn.astype(np.int32) - - -def _get_obs_type(swcm): - """ - Determine the observation type based on `swcm` and `chanfreq`. - - Parameters: - swcm (array-like): Switch mode values. - chanfreq (array-like): Channel frequency values (Hz). - - Returns: - numpy.ndarray: Observation type array. - - Raises: - ValueError: If any `obstype` is unassigned. - """ - - obstype = swcm.copy() - - # Use numpy vectorized operations - obstype = np.where(swcm == 1, 260, obstype) # IRLW - - if not np.any(np.isin(obstype, [260])): - raise ValueError("Error: Unassigned ObsType found ... ") - - return obstype - - -def _make_obs(comm, input_path, mapping_path): - - # Get container from mapping file first - logging(comm, 'INFO', 'Get container from bufr') - container = bufr.Parser(input_path, mapping_path).parse(comm) - logging(comm, 'DEBUG', f'container list (original): {container.list()}') - logging(comm, 'DEBUG', f'all_sub_categories = {container.all_sub_categories()}') - logging(comm, 'DEBUG', f'category map = {container.get_category_map()}') - - # Add new/derived data into container - for cat in container.all_sub_categories(): - - logging(comm, 'DEBUG', f'category = {cat}') - - satid = container.get('variables/satelliteId', cat) - if satid.size == 0: - logging(comm, 'WARNING', f'category {cat[0]} does not exist in input file') - paths = container.get_paths('variables/windComputationMethod', cat) - obstype = container.get('variables/windComputationMethod', cat) - container.add('variables/obstype_uwind', obstype, paths, cat) - container.add('variables/obstype_vwind', obstype, paths, cat) - - paths = container.get_paths('variables/windSpeed', cat) - wob = container.get('variables/windSpeed', cat) - container.add('variables/windEastward', wob, paths, cat) - container.add('variables/windNorthward', wob, paths, cat) - - paths = container.get_paths('variables/windComputationMethod', cat) - dummy = container.get('variables/windSpeed', cat) - container.add('variables/windGeneratingApplication', dummy, paths, cat) - container.add('variables/qualityInformationWithoutForecast', dummy, paths, cat) - - else: - # Add new variables: ObsType/windEastward & ObsType/windNorthward - swcm = container.get('variables/windComputationMethod', cat) - chanfreq = container.get('variables/sensorCentralFrequency', cat) - - logging(comm, 'DEBUG', f'swcm min/max = {swcm.min()} {swcm.max()}') - logging(comm, 'DEBUG', f'chanfreq min/max = {chanfreq.min()} {chanfreq.max()}') - - obstype = _get_obs_type(swcm) - - logging(comm, 'DEBUG', f'obstype = {obstype}') - logging(comm, 'DEBUG', f'obstype min/max = {obstype.min()} {obstype.max()}') - - paths = container.get_paths('variables/windComputationMethod', cat) - container.add('variables/obstype_uwind', obstype, paths, cat) - container.add('variables/obstype_vwind', obstype, paths, cat) - - # Add new variables: ObsValue/windEastward & ObsValue/windNorthward - wdir = container.get('variables/windDirection', cat) - wspd = container.get('variables/windSpeed', cat) - - logging(comm, 'DEBUG', f'wdir min/max = {wdir.min()} {wdir.max()}') - logging(comm, 'DEBUG', f'wspd min/max = {wspd.min()} {wspd.max()}') - - uob, vob = compute_wind_components(wdir, wspd) - - logging(comm, 'DEBUG', f'uob min/max = {uob.min()} {uob.max()}') - logging(comm, 'DEBUG', f'vob min/max = {vob.min()} {vob.max()}') - - paths = container.get_paths('variables/windSpeed', cat) - container.add('variables/windEastward', uob, paths, cat) - container.add('variables/windNorthward', vob, paths, cat) - - # Add new variables: MetaData/windGeneratingApplication and qualityInformationWithoutForecast - gnap2D = container.get('variables/generatingApplication', cat) - pccf2D = container.get('variables/qualityInformation', cat) - - gnap, qifn = _get_quality_information_and_generating_application(comm, gnap2D, pccf2D) - - logging(comm, 'DEBUG', f'gnap min/max = {gnap.min()} {gnap.max()}') - logging(comm, 'DEBUG', f'qifn min/max = {qifn.min()} {qifn.max()}') - - paths = container.get_paths('variables/windComputationMethod', cat) - container.add('variables/windGeneratingApplication', gnap, paths, cat) - container.add('variables/qualityInformationWithoutForecast', qifn, paths, cat) - - # Check - logging(comm, 'DEBUG', f'container list (updated): {container.list()}') - logging(comm, 'DEBUG', f'all_sub_categories {container.all_sub_categories()}') - - return container - - -def create_obs_group(input_path, mapping_path, category, env): - - comm = bufr.mpi.Comm(env["comm_name"]) - - description = _make_description(mapping_path, update=True) - - # Check the cache for the data and return it if it exists - logging(comm, 'DEBUG', f'Check if bufr.DataCache exists? {bufr.DataCache.has(input_path, mapping_path)}') - if bufr.DataCache.has(input_path, mapping_path): - container = bufr.DataCache.get(input_path, mapping_path) - logging(comm, 'INFO', f'Encode {category} from cache') - data = iodaEncoder(description).encode(container)[(category,)] - logging(comm, 'INFO', f'Mark {category} as finished in the cache') - bufr.DataCache.mark_finished(input_path, mapping_path, [category]) - logging(comm, 'INFO', f'Return the encoded data for {category}') - return data - - container = _make_obs(comm, input_path, mapping_path) - - # Gather data from all tasks into all tasks. Each task will have the complete record - logging(comm, 'INFO', f'Gather data from all tasks into all tasks') - container.all_gather(comm) - - logging(comm, 'INFO', f'Add container to cache') - # Add the container to the cache - bufr.DataCache.add(input_path, mapping_path, container.all_sub_categories(), container) - - # Encode the data - logging(comm, 'INFO', f'Encode {category}') - data = iodaEncoder(description).encode(container)[(category,)] - - logging(comm, 'INFO', f'Mark {category} as finished in the cache') - # Mark the data as finished in the cache - bufr.DataCache.mark_finished(input_path, mapping_path, [category]) +import bufr +from bufr.obs_builder import add_main_functions +from bufr_satwnd_amv_obs_builder import SatWndAmvObsBuilder, map_path - logging(comm, 'INFO', f'Return the encoded data for {category}') - return data +MAPPING_PATH = map_path('bufr_satwnd_amv_viirs.yaml') +FIND_QI = 5 -def create_obs_file(input_path, mapping_path, output_path): - comm = bufr.mpi.Comm("world") - container = _make_obs(comm, input_path, mapping_path) - container.gather(comm) +class SatWndAmvViirsObsBuilder(SatWndAmvObsBuilder): + def __init__(self): + super().__init__(MAPPING_PATH, log_name=os.path.basename(__file__)) - description = _make_description(mapping_path, update=True) + def make_obs(self, comm, input_path): + container = super().make_obs(comm, input_path) - # Encode the data - if comm.rank() == 0: - netcdfEncoder(description).encode(container, output_path) + # Add new/derived data into container + for cat in container.all_sub_categories(): + self._add_quality_info_and_gen_app(FIND_QI, container, cat) - logging(comm, 'INFO', f'Return the encoded data') + return container + def _make_description(self): + description = super()._make_description() + self._add_quality_info_and_gen_app_descriptions(description) -if __name__ == '__main__': + return description - start_time = time.time() + def _get_obs_type(self, swcm, chan_freq): + obstype = swcm.copy() - bufr.mpi.App(sys.argv) - comm = bufr.mpi.Comm("world") + # Use numpy vectorized operations + obstype = np.where(swcm == 1, 260, obstype) # IRLW - # Required input arguments as positional arguments - parser = argparse.ArgumentParser(description="Convert BUFR to NetCDF using a mapping file.") - parser.add_argument('input', type=str, help='Input BUFR file') - parser.add_argument('mapping', type=str, help='BUFR2IODA Mapping File') - parser.add_argument('output', type=str, help='Output NetCDF file') + if not np.any(np.isin(obstype, [260])): + raise ValueError("Error: Unassigned ObsType found ... ") - args = parser.parse_args() - infile = args.input - mapping = args.mapping - output = args.output + return obstype - create_obs_file(infile, mapping, output) - end_time = time.time() - running_time = end_time - start_time - logging(comm, 'INFO', f'Total running time: {running_time}') +# Add main functions create_obs_file and create_obs_group +add_main_functions(SatWndAmvViirsObsBuilder) diff --git a/dump/mapping/bufr_satwnd_amv_viirs_mapping.yaml b/dump/mapping/bufr_satwnd_amv_viirs.yaml similarity index 100% rename from dump/mapping/bufr_satwnd_amv_viirs_mapping.yaml rename to dump/mapping/bufr_satwnd_amv_viirs.yaml diff --git a/dump/mapping/bufr_satwnd_ascat.py b/dump/mapping/bufr_satwnd_ascat.py new file mode 100644 index 0000000..4fccb03 --- /dev/null +++ b/dump/mapping/bufr_satwnd_ascat.py @@ -0,0 +1,141 @@ +#!/usr/bin/env python3 + +import os +import numpy as np + +import bufr +from bufr.obs_builder import ObsBuilder, add_main_functions, map_path +from bufr.transforms import compute_wind_components + + +MAPPING_PATH = map_path('bufr_satwnd_ascat.yaml') + + +class BufrAscatObsBuilder(ObsBuilder): + """ + A builder class to generate satellite wind observations from ASCAT BUFR input. + + Inherits from :class:`bufr.obs_builder.ObsBuilder` and uses a mapping file to + extract wind speed and direction, compute wind vector components, and attach + observation metadata. + """ + + def __init__(self): + super().__init__(MAPPING_PATH, log_name=os.path.basename(__file__)) + + def make_obs(self, comm, input_path): + + # Get container from mapping file first + self.log.info('Get container from bufr') + container = super().make_obs(comm, input_path) + + self.log.debug(f'container list (original): {container.list()}') + self.log.debug(f'all_sub_categories = {container.all_sub_categories()}') + self.log.debug(f'category map = {container.get_category_map()}') + + # Add new/derived data into container + for cat in container.all_sub_categories(): + + self.log.debug(f'category = {cat}') + + satId = container.get('satelliteId', cat) + if not np.any(satId): + self.log.warning(f'category {cat[0]} does not exist in input file') + + paths = container.get_paths('windSpeedAt10M', cat) + dummy = container.get('windSpeedAt10M', cat) + container.add('windEastward', dummy, paths, cat) + container.add('windNorthward', dummy, paths, cat) + + paths = container.get_paths('satelliteId', cat) + dummy = container.get('satelliteId', cat) + container.add('obstype_windEastward', dummy, paths, cat) + container.add('obstype_windNorthward', dummy, paths, cat) + continue + + wdir = container.get('windDirectionAt10M', cat) + wspd = container.get('windSpeedAt10M', cat) + self.log.debug(f'wdir min/max = {wdir.min()} {wdir.max()}') + self.log.debug(f'wspd min/max = {wspd.min()} {wspd.max()}') + + uob, vob = compute_wind_components(wspd, wdir) + self.log.debug(f'uob min/max = {uob.min()} {uob.max()}') + self.log.debug(f'vob min/max = {vob.min()} {vob.max()}') + + paths = container.get_paths('windSpeedAt10M', cat) + container.add('windEastward', uob, paths, cat) + container.add('windNorthward', vob, paths, cat) + + obstype = self._get_obs_type(container, cat) + paths = container.get_paths('satelliteId', cat) + container.add('obstype_windEastward', obstype, paths, cat) + container.add('obstype_windNorthward', obstype, paths, cat) + + # Check + self.log.debug(f'container list (updated): {container.list()}') + self.log.debug(f'all_sub_categories {container.all_sub_categories()}') + + return container + + def _get_obs_type(self, container, category): + """ + Retrieve observation type values for wind components. + + :param container: Observation container from BUFR input. + :type container: bufr.DataContainer + :param category: Subcategory identifier (e.g., sensor/platform group). + :type category: tuple + + :returns: Array filled with observation type (290) for the given category. + :rtype: numpy.ndarray + """ + + satId = container.get('satelliteId', category) + + if not satId.size: + paths = container.get_paths('satelliteId', cat) + dummy = container.get('satelliteId', cat) + container.add('obstype_windEastward', dummy, paths, cat) + container.add('obstype_windNorthward', dummy, paths, cat) + return + + obstype = np.full_like(satId, 290) + + return obstype + + def _make_description(self): + description = super()._make_description() + self._add_new_variable_descriptions(description) + + return description + + def _add_new_variable_descriptions(self, description): + description.add_variables([ + { + 'name': 'ObsValue/windEastward', + 'source': 'windEastward', + 'units': 'm s-1', + 'longName': '10-meter U-Wind Component', + }, + { + 'name': 'ObsValue/windNorthward', + 'source': 'windNorthward', + 'units': 'm s-1', + 'longName': '10-meter V-Wind Component', + }, + { + 'name': 'ObsType/windEastward', + 'source': 'obstype_windEastward', + 'units': '1', + 'longName': 'Observation Type for Wind Components', + }, + { + 'name': 'ObsType/windNorthward', + 'source': 'obstype_windNorthward', + 'units': '1', + 'longName': 'Observation Type for Wind Components', + }]) + + +# Add main functions create_obs_file or create_obs_group +add_main_functions(BufrAscatObsBuilder) diff --git a/dump/mapping/bufr_satwnd_ascat.yaml b/dump/mapping/bufr_satwnd_ascat.yaml new file mode 100755 index 0000000..9012635 --- /dev/null +++ b/dump/mapping/bufr_satwnd_ascat.yaml @@ -0,0 +1,128 @@ +bufr: + subsets: + - NC012122 + + variables: + # MetaData + timestamp: + datetime: + year: "*/YEAR" + month: "*/MNTH" + day: "*/DAYS" + hour: "*/HOUR" + minute: "*/MINU" + second: "*/SECO" + + dataReceiptTime: + datetime: + year: "*/RCYR" + month: "*/RCMO" + day: "*/RCDY" + hour: "*/RCHR" + minute: "*/RCMI" + + latitude: + query: "*/CLAT" + + longitude: + query: "*/CLON" + + satelliteId: + query: "*/SAID" + + qualityFlags: + query: "*/WVCQ" + + # ObsValue - Wind Direction + windDirectionAt10M: + query: '*/WD10' + + # ObsValue - Wind Speed + windSpeedAt10M: + query: '*/WS10' + + splits: + satId: + category: + variable: satelliteId + map: + _3: metop-b + _4: metop-a + _5: metop-c + +encoder: +# type: netcdf + +# dimensions: + + globals: + - name: "platformCommonName" + type: string + value: "MetOp" + + - name: "platformLongDescription" + type: string + value: "Meteorological Operational Satellite" + + - name: "sensor" + type: string + value: "Advanced Scatterometer - ASCATW" + + - name: "source" + type: string + value: "NCEP BUFR Dump for satellite derived surface wind components from ASCAT (scatw)" + + - name: "providerFullName" + type: string + value: "National Environmental Satellite, Data, and Information Service" + + - name: "processingLevel" + type: string + value: "Level-2" + + - name: "converter" + type: string + value: "BUFR" + + variables: + # MetaData + - name: "MetaData/dateTime" + source: variables/timestamp + longName: "Observation Time" + units: "seconds since 1970-01-01T00:00:00Z" + + - name: "MetaData/dataReceiptTime" + source: variables/dataReceiptTime + longName: "Observation Receipt Time" + units: "seconds since 1970-01-01T00:00:00Z" + + - name: "MetaData/latitude" + source: variables/latitude + longName: "Latitude" + units: "degrees_north" + range: [ -90, 90 ] + + - name: "MetaData/longitude" + source: variables/longitude + longName: "Longitude" + units: "degrees_east" + range: [ -180, 180 ] + + - name: "MetaData/satelliteIdentifier" + source: variables/satelliteId + longName: "Satellite Identifier" + + - name: "MetaData/qualityFlags" + source: variables/qualityFlags + longName: "Quality Flags" + units: "1" + + - name: "ObsValue/windDirection" + source: variables/windDirectionAt10M + longName: "10-meter Wind Direction" + units: "degree" + + - name: "ObsValue/windSpeed" + source: variables/windSpeedAt10M + longName: "10-meter Wind Speed" + units: "m s-1" diff --git a/dump/mapping/bufr_sfcshp_prepbufr.py b/dump/mapping/bufr_sfcshp_prepbufr.py deleted file mode 100755 index 8055289..0000000 --- a/dump/mapping/bufr_sfcshp_prepbufr.py +++ /dev/null @@ -1,311 +0,0 @@ -#!/usr/bin/env python3 -import os -import sys -import bufr -import argparse -import copy -import numpy as np -import numpy.ma as ma -import math -import calendar -import time -from datetime import datetime -from pyioda.ioda.Engines.Bufr import Encoder as iodaEncoder -from bufr.encoders.netcdf import Encoder as netcdfEncoder -from wxflow import Logger - - -# Initialize Logger -# Get log level from the environment variable, default to 'INFO it not set -log_level = os.getenv('LOG_LEVEL', 'INFO') -logger = Logger('bufr_sfcshp_prepbufr.py', level=log_level, colored_log=False) - - -def logging(comm, level, message): - """ - Logs a message to the console or log file, based on the specified logging level. - - This function ensures that logging is only performed by the root process (`rank 0`) - in a distributed computing environment. The function maps the logging level to - appropriate logger methods and defaults to the 'INFO' level if an invalid level is provided. - - Parameters: - comm: object - The communicator object, typically from a distributed computing framework - (e.g., MPI). It must have a `rank()` method to determine the process rank. - level: str - The logging level as a string. Supported levels are: - - 'DEBUG' - - 'INFO' - - 'WARNING' - - 'ERROR' - - 'CRITICAL' - If an invalid level is provided, a warning will be logged, and the level - will default to 'INFO'. - message: str - The message to be logged. - - Behavior: - - Logs messages only on the root process (`comm.rank() == 0`). - - Maps the provided logging level to a method of the logger object. - - Defaults to 'INFO' and logs a warning if an invalid logging level is given. - - Supports standard logging levels for granular control over log verbosity. - - Example: - >>> logging(comm, 'DEBUG', 'This is a debug message.') - >>> logging(comm, 'ERROR', 'An error occurred!') - - Notes: - - Ensure that a global `logger` object is configured before using this function. - - The `comm` object should conform to MPI-like conventions (e.g., `rank()` method). - """ - - if comm.rank() == 0: - # Define a dictionary to map levels to logger methods - log_methods = { - 'DEBUG': logger.debug, - 'INFO': logger.info, - 'WARNING': logger.warning, - 'ERROR': logger.error, - 'CRITICAL': logger.critical, - } - - # Get the appropriate logging method, default to 'INFO' - log_method = log_methods.get(level.upper(), logger.info) - - if log_method == logger.info and level.upper() not in log_methods: - # Log a warning if the level is invalid - logger.warning(f'log level = {level}: not a valid level --> set to INFO') - - # Call the logging method - log_method(message) - - -def _make_description(mapping_path, cycle_time, update=False): - description = bufr.encoders.Description(mapping_path) - - ReferenceTime = np.int64(calendar.timegm(time.strptime(str(int(cycle_time)), '%Y%m%d%H'))) - - if update: - # Define the variables to be added in a list of dictionaries - variables = [ - { - 'name': 'MetaData/sequenceNumber', - 'source': 'variables/sequenceNumber', - 'units': '1', - 'longName': 'Sequence Number (Obs Subtype)', - }, - ] - - # Loop through each variable and add it to the description - for var in variables: - description.add_variable( - name=var['name'], - source=var['source'], - units=var['units'], - longName=var['longName'] - ) - - # description.add_global(name='Reference_time', value=str(ReferenceTime)) - - return description - - -def _compute_sequence_number(typ, t29): - """ - Compute sequenceNumber - - Parameters: - typ: observation Type (obsType) - t29: data dump report type - - Returns: - Masked array of sequenceNumber values - """ - - sequenceNumber = np.zeros(typ.shape, dtype=np.int32) - for i in range(len(typ)): - if (typ[i] == 180 or typ[i] == 280): - if (t29[i] > 555 and t29[i] < 565): - sequenceNumber[i] = 0 - else: - sequenceNumber[i] = 1 - - return sequenceNumber - - -def _compute_datetime(cycleTimeSinceEpoch, dhr): - """ - Compute dateTime using the cycleTimeSinceEpoch and Cycle Time - minus Cycle Time - - Parameters: - cycleTimeSinceEpoch: Time of cycle in Epoch Time - dhr: Observation Time Minus Cycle Time - - Returns: - Masked array of dateTime values - """ - - int64_fill_value = np.int64(0) - - dateTime = np.zeros(dhr.shape, dtype=np.int64) - for i in range(len(dateTime)): - if ma.is_masked(dhr[i]): - continue - else: - dateTime[i] = np.int64(dhr[i]*3600) + cycleTimeSinceEpoch - - dateTime = ma.array(dateTime) - dateTime = ma.masked_values(dateTime, int64_fill_value) - - return dateTime - - -def _make_obs(comm, input_path, mapping_path, cycle_time): - """ - Create the ioda sfcshp prepbufr observations: - - reads values - - adds sequenceNum - - Parameters - ---------- - comm: object - The communicator object (e.g., MPI) - input_path: str - The input bufr file - mapping_path: str - The input bufr2ioda mapping file - cycle_time: str - The cycle in YYYYMMDDHH format - """ - - # Get container from mapping file first - logging(comm, 'INFO', 'Get container from bufr') - container = bufr.Parser(input_path, mapping_path).parse(comm) - - logging(comm, 'DEBUG', f'container list (original): {container.list()}') - logging(comm, 'DEBUG', f'Change longitude range from [0,360] to [-180,180]') - lon = container.get('variables/longitude') - lon_paths = container.get_paths('variables/longitude') - lon[lon > 180] -= 360 - lon = ma.round(lon, decimals=2) - logging(comm, 'DEBUG', f'longitude new max/min: ${lon.max()}, ${lon.min()}') - - logging(comm, 'DEBUG', f'Do DateTime calculation') - otmct = container.get('variables/obsTimeMinusCycleTime') - otmct_paths = container.get_paths('variables/obsTimeMinusCycleTime') - otmct2 = np.array(otmct) - cycleTimeSinceEpoch = np.int64(calendar.timegm(time.strptime(str(int(cycle_time)), '%Y%m%d%H'))) - dateTime = _compute_datetime(cycleTimeSinceEpoch, otmct2) - min_dateTime_ge_zero = min(x for x in dateTime if x >= 0) - logging(comm, 'DEBUG', f'dateTime min/max = {min_dateTime_ge_zero} {dateTime.max()}') - - logging(comm, 'DEBUG', f'Do sequenceNumber (Obs SubType) calculation') - typ = container.get('variables/observationType') - typ_paths = container.get_paths('variables/observationType') - t29 = container.get('variables/obssubtype') - t29_paths = container.get_paths('variables/obssubtype') - seqNum = _compute_sequence_number(typ, t29) - logging(comm, 'DEBUG', f' sequenceNum min/max = {seqNum.min()} {seqNum.max()}') - - logging(comm, 'DEBUG', f'Do tsen and tv calculation') - tpc = container.get('variables/temperatureEventCode') - tob = container.get('variables/airTemperatureObsValue') - tob_paths = container.get_paths('variables/airTemperatureObsValue') - tsen = np.full(tob.shape[0], tob.fill_value) - tsen = np.where(((tpc >= 1) & (tpc < 8)), tob, tsen) - tvo = np.full(tob.shape[0], tob.fill_value) - tvo = np.where((tpc == 8), tob, tvo) - - logging(comm, 'DEBUG', f'Do tsen and tv QM calculations') - tobqm = container.get('variables/airTemperatureQualityMarker') - tsenqm = np.full(tobqm.shape[0], tobqm.fill_value) - tsenqm = np.where(((tpc >= 1) & (tpc < 8)), tobqm, tsenqm) - tvoqm = np.full(tobqm.shape[0], tobqm.fill_value) - tvoqm = np.where((tpc == 8), tobqm, tvoqm) - - logging(comm, 'DEBUG', f'Do tsen and tv ObsError calculations') - toboe = container.get('variables/airTemperatureObsError') - tsenoe = np.full(toboe.shape[0], toboe.fill_value) - tsenoe = np.where(((tpc >= 1) & (tpc < 8)), toboe, tsenoe) - tvooe = np.full(toboe.shape[0], toboe.fill_value) - tvooe = np.where((tpc == 8), toboe, tvooe) - - logging(comm, 'DEBUG', f'Update variables in container') - container.replace('variables/longitude', lon) - container.replace('variables/timestamp', dateTime) - container.replace('variables/airTemperatureObsValue', tsen) - container.replace('variables/airTemperatureQualityMarker', tsenqm) - container.replace('variables/airTemperatureObsError', tsenoe) - container.replace('variables/virtualTemperatureObsValue', tvo) - container.replace('variables/virtualTemperatureQualityMarker', tvoqm) - container.replace('variables/virtualTemperatureObsError', tvooe) - - logging(comm, 'DEBUG', f'Add variables to container') - container.add('variables/sequenceNumber', seqNum, typ_paths) - - # Check - logging(comm, 'DEBUG', f'container list (updated): {container.list()}') - - return container - - -def create_obs_group(input_path, mapping_path, cycle_time, env): - - comm = bufr.mpi.Comm(env["comm_name"]) - - logging(comm, 'INFO', f'Make description and make obs') - description = _make_description(mapping_path, cycle_time, update=True) - container = _make_obs(comm, input_path, mapping_path, cycle_time) - - # Gather data from all tasks into all tasks. Each task will have the complete record - logging(comm, 'INFO', f'Gather data from all tasks into all tasks') - container.all_gather(comm) - - logging(comm, 'INFO', f'Encode the data') - data = next(iter(iodaEncoder(description).encode(container).values())) - - logging(comm, 'INFO', f'Return the encoded data.') - return data - - -def create_obs_file(input_path, mapping_path, output_path, cycle_time): - - comm = bufr.mpi.Comm("world") - container = _make_obs(comm, input_path, mapping_path, cycle_time) - container.gather(comm) - - description = _make_description(mapping_path, cycle_time, update=True) - - # Encode the data - if comm.rank() == 0: - netcdfEncoder(description).encode(container, output_path) - - logging(comm, 'INFO', f'Return the encoded data') - - -if __name__ == '__main__': - start_time = time.time() - - bufr.mpi.App(sys.argv) - comm = bufr.mpi.Comm("world") - - # Required input arguments as positional arguments - parser = argparse.ArgumentParser(description="Convert BUFR to NetCDF using a mapping file.") - parser.add_argument('input', type=str, help='Input BUFR file') - parser.add_argument('mapping', type=str, help='BUFR2IODA Mapping File') - parser.add_argument('output', type=str, help='Output NetCDF file') - parser.add_argument('cycle_time', type=str, help='cycle time in YYYYMMDDHH format') - - args = parser.parse_args() - infile = args.input - mapping = args.mapping - output = args.output - cycle_time = args.cycle_time - - create_obs_file(infile, mapping, output, cycle_time) - - end_time = time.time() - running_time = end_time - start_time - logging(comm, 'INFO', f'Total running time: {running_time}') diff --git a/dump/mapping/bufr_sfcsno.py b/dump/mapping/bufr_sfcsno.py index cc1a5ce..bc208f2 100755 --- a/dump/mapping/bufr_sfcsno.py +++ b/dump/mapping/bufr_sfcsno.py @@ -1,190 +1,37 @@ #!/usr/bin/env python3 -import sys import os -import argparse -import time import numpy as np -import bufr -from pyioda.ioda.Engines.Bufr import Encoder as iodaEncoder -from bufr.encoders.netcdf import Encoder as netcdfEncoder -from wxflow import Logger - -# Initialize Logger -# Get log level from the environment variable, default to 'INFO it not set -log_level = os.getenv('LOG_LEVEL', 'INFO') -logger = Logger('BUFR2IODA_sfcsno.py', level=log_level, colored_log=False) - - -def logging(comm, level, message): - """ - Logs a message to the console or log file, based on the specified logging level. - - This function ensures that logging is only performed by the root process (`rank 0`) - in a distributed computing environment. The function maps the logging level to - appropriate logger methods and defaults to the 'INFO' level if an invalid level is provided. - - Parameters: - comm: object - The communicator object, typically from a distributed computing framework - (e.g., MPI). It must have a `rank()` method to determine the process rank. - level: str - The logging level as a string. Supported levels are: - - 'DEBUG' - - 'INFO' - - 'WARNING' - - 'ERROR' - - 'CRITICAL' - If an invalid level is provided, a warning will be logged, and the level - will default to 'INFO'. - message: str - The message to be logged. - - Behavior: - - Logs messages only on the root process (`comm.rank() == 0`). - - Maps the provided logging level to a method of the logger object. - - Defaults to 'INFO' and logs a warning if an invalid logging level is given. - - Supports standard logging levels for granular control over log verbosity. - - Example: - >>> logging(comm, 'DEBUG', 'This is a debug message.') - >>> logging(comm, 'ERROR', 'An error occurred!') - - Notes: - - Ensure that a global `logger` object is configured before using this function. - - The `comm` object should conform to MPI-like conventions (e.g., `rank()` method). - """ - - if comm.rank() == 0: - # Define a dictionary to map levels to logger methods - log_methods = { - 'DEBUG': logger.debug, - 'INFO': logger.info, - 'WARNING': logger.warning, - 'ERROR': logger.error, - 'CRITICAL': logger.critical, - } - - # Get the appropriate logging method, default to 'INFO' - log_method = log_methods.get(level.upper(), logger.info) - - if log_method == logger.info and level.upper() not in log_methods: - # Log a warning if the level is invalid - logger.warning(f'log level = {level}: not a valid level --> set to INFO') - - # Call the logging method - log_method(message) - - -def _mask_container(container, mask): - - new_container = bufr.DataContainer() - for var_name in container.list(): - var = container.get(var_name) - paths = container.get_paths(var_name) - new_container.add(var_name, var[mask], paths) - - return new_container - - -def _make_description(mapping_path, update=False): - - description = bufr.encoders.Description(mapping_path) +import numpy.ma as ma - return description - - -def _make_obs(comm, input_path, mapping_path): - """ - Create the ioda snow depth observations: - - reads state of ground (sogr) and snow depth (snod) - - applys sogr conditions to the missing snod values - - removes the filled/missing snow values and creates the masked container - - Parameters - ---------- - comm: object - The communicator object (e.g., MPI) - input_path: str - The input bufr file - mapping_path: str - The input bufr2ioda mapping file - """ - - # Get container from mapping file first - logging(comm, 'INFO', 'Get container from bufr') - container = bufr.Parser(input_path, mapping_path).parse(comm) - - logging(comm, 'DEBUG', f'container list (original): {container.list()}') - - # Add new/derived data into container - sogr = container.get('variables/groundState') - snod = container.get('variables/totalSnowDepth') - snod[(sogr <= 11.0) & snod.mask] = 0.0 - snod[(sogr == 15.0) & snod.mask] = 0.0 - snod.mask = (snod < 0.0) | snod.mask - container.replace('variables/totalSnowDepth', snod) - snod_upd = container.get('variables/totalSnowDepth') - - masked_container = _mask_container(container, (~snod.mask)) - - return masked_container - - -def create_obs_group(input_path, mapping_path, env): - - comm = bufr.mpi.Comm(env["comm_name"]) - - description = _make_description(mapping_path, update=False) - container = _make_obs(comm, input_path, mapping_path) - - # Gather data from all tasks into all tasks. Each task will have the complete record - logging(comm, 'INFO', f'Gather data from all tasks into all tasks') - container.all_gather(comm) - - # Encode the data - logging(comm, 'INFO', f'Encode data') - data = next(iter(iodaEncoder(mapping_path).encode(container).values())) - - logging(comm, 'INFO', f'Return the encoded data') - - return data - - -def create_obs_file(input_path, mapping_path, output_path): - - comm = bufr.mpi.Comm("world") - container = _make_obs(comm, input_path, mapping_path) - container.gather(comm) +import bufr +from bufr.obs_builder import ObsBuilder, add_main_functions - description = _make_description(mapping_path, update=False) - # Encode the data - if comm.rank() == 0: - netcdfEncoder(description).encode(container, output_path) +MAPPING_PATH = map_path('bufr_sfcsno.yaml') - logging(comm, 'INFO', f'Return the encoded data') +class BufrSfcsnoObsBuilder(ObsBuilder): + def __init__(self): + super().__init__(MAPPING_PATH, log_name=os.path.basename(__file__)) -if __name__ == '__main__': + def make_obs(self, comm, input_path): + self.log.info('Get container from bufr') + container = super().make_obs(comm, input_path) - start_time = time.time() + self.log.debug(f'container list (original): {container.list()}') - bufr.mpi.App(sys.argv) - comm = bufr.mpi.Comm("world") + # Add new/derived data into container + sogr = container.get('groundState') + snod = container.get('totalSnowDepth') + snod[(sogr <= 11.0) & snod.mask] = 0.0 + snod[(sogr == 15.0) & snod.mask] = 0.0 + snod.mask = (snod < 0.0) | snod.mask + container.replace('totalSnowDepth', snod) + snod_upd = container.get('totalSnowDepth') - # Required input arguments as positional arguments - parser = argparse.ArgumentParser(description="Convert BUFR to NetCDF using a mapping file.") - parser.add_argument('input', type=str, help='Input BUFR file') - parser.add_argument('mapping', type=str, help='BUFR2IODA Mapping File') - parser.add_argument('output', type=str, help='Output NetCDF file') + container.apply_mask(~snod.mask) - args = parser.parse_args() - mapping = args.mapping - infile = args.input - output = args.output + return container - create_obs_file(infile, mapping, output) - end_time = time.time() - running_time = end_time - start_time - logging(comm, 'INFO', f'Total running time: {running_time}') +add_main_functions(BufrSfcsnoObsBuilder) diff --git a/dump/mapping/bufr_sfcsno_mapping.yaml b/dump/mapping/bufr_sfcsno.yaml similarity index 100% rename from dump/mapping/bufr_sfcsno_mapping.yaml rename to dump/mapping/bufr_sfcsno.yaml diff --git a/dump/mapping/bufr_ssmis.py b/dump/mapping/bufr_ssmis.py new file mode 100755 index 0000000..af7d76c --- /dev/null +++ b/dump/mapping/bufr_ssmis.py @@ -0,0 +1,176 @@ +#!/usr/bin/env python3 +import os +import numpy as np +import numpy.ma as ma + +import bufr +from bufr.obs_builder import ObsBuilder, add_main_functions, map_path +from bufr.obs_builder import nprocs_per_task, add_dummy_variable +from bufr.transforms import compute_solar_angles +from datetime import datetime + +MAPPING_PATH = map_path('bufr_ssmis.yaml') + + +class BufrSsmisObsBuilder(ObsBuilder): + """ + Class for building observations from ssmis BUFR data. + + This class extends `ObsBuilder` to include specific logic for processing + SSMIS data such as solar angles and satellite ascending/descending orbits. + + :param mapping_path: Path to the mapping file. + :type mapping_path: str + """ + + def __init__(self): + super().__init__(MAPPING_PATH, log_name=os.path.basename(__file__)) + + def make_obs(self, comm, input_path): + + # Get container from mapping file first + self.log.info('Get container from bufr') + container = super().make_obs(comm, input_path) + + self.log.debug(f'container list (original): {container.list()}') + self.log.debug(f'all_sub_categories = {container.all_sub_categories()}') + self.log.debug(f'category map = {container.get_category_map()}') + + # Add new/derived data into container + for cat in container.all_sub_categories(): + + self.log.debug(f'category = {cat}') + + satId = container.get('satelliteId', cat) + if not np.any(satId): + self.log.warning(f'category {cat[0]} does not exist in input file') + + self._add_sensor_zenith_and_solar_angles(container, cat) + self._add_satellite_ascend_descent_orbit(container, cat) + + # Check + self.log.debug(f'container list (updated): {container.list()}') + self.log.debug(f'all_sub_categories {container.all_sub_categories()}') + + return container + + def _make_description(self): + description = super()._make_description() + self._add_new_variable_descriptions(description) + + return description + + def _add_new_variable_descriptions(self, description): + description.add_variables([ + { + 'name': 'MetaData/satelliteAscendingFlag', + 'source': 'satelliteAscendingFlag', + 'longName': 'Satellite Ascending/Descending Orbit Flag (Ascend:1; Descend:-1)', + }, + { + 'name': 'MetaData/sensorZenithAngle', + 'source': 'sensorZenithAngle', + 'units': 'degree', + 'longName': 'Sensor Zenith Angle', + }, + { + 'name': 'MetaData/sensorAzimuthAngle', + 'source': 'sensorAzimuthAngle', + 'units': 'degree', + 'longName': 'Sensor Azimuth Angle', + }, + { + 'name': 'MetaData/solarZenithAngle', + 'source': 'solarZenithAngle', + 'units': 'degree', + 'longName': 'Solar Zenith Angle', + }, + { + 'name': 'MetaData/solarAzimuthAngle', + 'source': 'solarAzimuthAngle', + 'units': 'degree', + 'longName': 'Solar Azimuth Angle', + }]) + + def _add_satellite_ascend_descent_orbit(self, container, category): + """ + Determine satellite orbit type (ascending or descending) based on latitude changes. + + :param container: Observation data container. + :type container: Container + :param category: Data category to process. + :type category: str + """ + + satId = container.get('satelliteId', category) + + if not satId.size: + add_dummy_variable(container, 'satelliteAscendingFlag', category, 'fieldOfViewNumber') + return + + # Get data from container + # ephemeris data - latitude values in order of time + first_lat = container.get('latitude1', category) + self.log.debug(f'first_lat min/max = {first_lat.min()} {first_lat.max()}') + second_lat = container.get('latitude2', category) + self.log.debug(f'second_lat min/max = {second_lat.min()} {second_lat.max()}') + fovn = container.get('fieldOfViewNumber', category) + self.log.debug(f'fovn min/max = {fovn.min()} {fovn.max()}') + + # Determine ascending/descending mode + # Compare latitude between the first and second records + orbit = np.where(second_lat > first_lat, 1, -1).astype(np.int32) + + self.log.debug(f'orbit min/max = {orbit.min()} {orbit.max()}') + + paths = container.get_paths('fieldOfViewNumber', category) + self.log.debug(f'paths = {paths}') + container.add('satelliteAscendingFlag', orbit, paths, category) + + def _add_sensor_zenith_and_solar_angles(self, container, category): + """ + Compute and add solar zenith and azimuth angles to the observation container. + + :param container: Observation data container. + :type container: Container + :param category: Data category to process. + :type category: str + """ + + satId = container.get('satelliteId', category) + if not satId.size: + add_dummy_variable(container, 'solarZenithAngle', category, 'latitude') + add_dummy_variable(container, 'solarAzimuthAngle', category, 'latitude') + add_dummy_variable(container, 'sensorZenithAngle', category, 'latitude') + add_dummy_variable(container, 'sensorAzimuthAngle', category, 'latitude') + return + + # Prepare input arrays + unix_times = container.get('timestamp', category) + latitudes = container.get('latitude', category) + longitudes = container.get('longitude', category) + self.log.debug(f'latitudes min/max = {latitudes.min()} {latitudes.max()}') + self.log.debug(f'longitudes min/max = {longitudes.min()} {longitudes.max()}') + self.log.debug(f'unix_times min/max = {unix_times.min()} {unix_times.max()}') + + # Calculate solar angles + zenith_angles, azimuth_angles = compute_solar_angles(latitudes, longitudes, unix_times) + + self.log.debug(f'zenith_angles min/max = {zenith_angles.min()} {zenith_angles.max()}') + self.log.debug(f'azimuth_angles min/max = {azimuth_angles.min()} {azimuth_angles.max()}') + + # Add solar angles + paths = container.get_paths('latitude', category) + self.log.debug(f'paths = {paths}') + container.add('solarZenithAngle', zenith_angles, paths, category) + container.add('solarAzimuthAngle', azimuth_angles, paths, category) + + # Add sensor angles + sensor_zenith = np.full_like(latitudes, 53.0) + sensor_azimuth = np.full_like(latitudes, latitudes.fill_value) + container.add('sensorZenithAngle', sensor_zenith, paths, category) + container.add('sensorAzimuthAngle', sensor_azimuth, paths, category) + + +# Add main functions create_obs_file or create_obs_group +add_main_functions(BufrSsmisObsBuilder) diff --git a/dump/mapping/bufr_ssmis.yaml b/dump/mapping/bufr_ssmis.yaml new file mode 100644 index 0000000..2b4712c --- /dev/null +++ b/dump/mapping/bufr_ssmis.yaml @@ -0,0 +1,196 @@ +bufr: + variables: + # MetaData + timestamp: + datetime: + year: "*/YEAR" + month: "*/MNTH" + day: "*/DAYS" + hour: "*/HOUR" + minute: "*/MINU" + second: "*/SECO" + + latitude: + query: "*/CLAT" + + latitude1: + query: "*/SATEPHEM{1}/CLATH" + + latitude2: + query: "*/SATEPHEM{2}/CLATH" + + longitude: + query: "*/CLON" + + satelliteId: + query: "*/SAID" + + scanLineNumber: + query: "*/SLNM" + + fieldOfViewNumber: + query: "*/FOVN" + + surfaceFlag: + query: "*/SFLG" + + rainFlag: + query: "*/RFLAG" + + sensorChannelNumber: + query: "*/SSMISCHN/CHNM" + + sensorScanAngle: + sensorScanAngle: + fieldOfViewNumber: "*/FOVN" + scanStart: 0.0 + scanStep: 0.0 + sensor: ssmis + + remappedBT: + remappedBrightnessTemperature: + sensor: ssmis + method: 1 + satelliteId: "*/SAID" + latitude: "*/CLAT" + longitude: "*/CLON" + fieldOfViewNumber: "*/FOVN" + rainFlag: "*/RFLAG" + sensorChannelNumber: "*/SSMISCHN/CHNM" + brightnessTemperature: "*/SSMISCHN/TMBR" + obsTime: + year: "*/YEAR" + month: "*/MNTH" + day: "*/DAYS" + hour: "*/HOUR" + minute: "*/MINU" + second: "*/SECO" + + splits: + satId: + category: + variable: satelliteId + map: + _249: f16 + _285: f17 + _286: f18 + +encoder: + dimensions: + - name: Channel + source: variables/sensorChannelNumber + path: "*/SSMISCHN" + + globals: + - name: "description" + type: string + value: "SSMIS Level-1B calibrated radiance data with precipitation-contaminated FOVs filtered out and spatial averaging applied to reduce spatial noise" + + - name: "platformCommonName" + type: string + value: "DSMP" + + - name: "platformLongDescription" + type: string + value: "Defense Meteorological Satellite Program" + + - name: "sensorCommonName" + type: string + value: "SSMIS" + + - name: "sensorLongDescription" + type: string + value: "Special Sensor Microwave Imager/Sounder (SSMIS) - 21 frequencies, 24 channels MW radiometer covering the 54 and 183 GHz bands; conical scanning" + + - name: "source" + type: string + value: "U.S. National Weather Service, National Centres for Environmental Prediction (NCEP))" + + - name: "sourceFiles" + type: string + value: "NCEP BUFR Dump - ssmisu NC021201" + + - name: "processingLevel" + type: string + value: "Level-1B" + + - name: "converter" + type: string + value: "bufr-query" + + variables: + # MetaData + - name: "MetaData/dateTime" + source: variables/timestamp + longName: "Datetime" + units: "seconds since 1970-01-01T00:00:00Z" + + - name: "MetaData/latitude" + source: variables/latitude + longName: "Latitude" + units: "degree_north" + range: [-90, 90] + +# - name: "MetaData/latitudeSatEphem1" +# source: variables/latitude1 +# longName: "Latitude from SATEPHEM Record 1 (high resolution)" +# units: "degree_north" +# range: [-90, 90] + +# - name: "MetaData/latitudeSatEphem2" +# source: variables/latitude2 +# longName: "Latitude from SATEPHEM Record 2 (high resolution)" +# units: "degree_north" +# range: [-90, 90] + + - name: "MetaData/longitude" + source: variables/longitude + longName: "Longitude" + units: "degree_east" + + - name: "MetaData/satelliteIdentifier" + source: variables/satelliteId + longName: "Satellite Identifier" + + - name: "MetaData/sensorScanPosition" + source: variables/fieldOfViewNumber + longName: "Field of View Number" + + - name: "MetaData/sensorChannelNumber" + source: variables/sensorChannelNumber + longName: "Sensor Channel Number" + + - name: "MetaData/sensorViewAngle" + source: variables/sensorScanAngle + longName: "Sensor View Angle" + units: "degree" + + - name: "MetaData/earthSurfaceType" + source: variables/surfaceFlag + longName: "Earth Surface Type" + + - name: "MetaData/rainFlag" + source: variables/rainFlag + longName: "Rain Flag (1:rain; 0:no rain (clear); -1:no determined" + + - name: "MetaData/scanLineNumber" + source: variables/scanLineNumber + longName: "Scan Line Number" + +# # ObsValue +# # Brightness Temperature +# - name: "ObsValue/brightnessTemperature" +# source: variables/brightnessTemperature +# longName: "Brightness Temperature" +# units: "K" +# range: [120, 500] +# chunks: [10000, 24] + + # ObsValue + # Remapped Brightness Temperature + - name: "ObsValue/brightnessTemperature" + source: variables/remappedBT + longName: "3-by-3 Averaged Brightness Temperature" + units: "K" + range: [120, 500] + chunks: [10000, 24] diff --git a/dump/mapping/prepbufr_adpsfc.py b/dump/mapping/prepbufr_adpsfc.py new file mode 100755 index 0000000..313aff2 --- /dev/null +++ b/dump/mapping/prepbufr_adpsfc.py @@ -0,0 +1,81 @@ +#!/usr/bin/env python3 +import calendar +import os +import numpy as np +import numpy.ma as ma + +import bufr +from bufr.obs_builder import add_main_functions +from prepbufr_obs_builder import PrepbufrObsBuilder, map_path + + +MAPPING_PATH = map_path('prepbufr_adpsfc.yaml') + + +class AdpsfcPrepbufrObsBuilder(PrepbufrObsBuilder): + def __init__(self): + super().__init__(MAPPING_PATH, log_name=os.path.basename(__file__)) + + def _make_description(): + description = super()._make_description() + + description.add_variables([ + { + 'name': "MetaData/dateTime", + 'source': 'timestamp', + 'units': "seconds since 1970-01-01T00:00:00Z", + 'longName': "Observation Time" + }, + { + 'name': 'MetaData/sequenceNumber', + 'source': 'sequenceNumber', + 'units': '1', + 'longName': 'Sequence Number (Obs Subtype)', + } + ]) + + return description + + def make_obs(self, comm, input_path): + """ + Create the ioda adpsfc prepbufr observations: + - reads values + - adds sequenceNum + + Parameters + ---------- + comm: object + The communicator object (e.g., MPI) + input_path: str + The input bufr file + """ + + container = super().make_obs(comm, input_path) + + # Get container from mapping file first + self.log.info('Get container from bufr') + container = bufr.Parser(input_path, mapping_path).parse(comm) + + self.log.debug(f'container list (original): {container.list()}') + + self.log.debug(f'Do DateTime calculation') + otmct = container.get('obsTimeMinusCycleTime') + otmct_paths = container.get_paths('obsTimeMinusCycleTime') + otmct2 = np.array(otmct) + + self._add_timestamp(container, self._get_reference_time(input_path)) + + self.log.debug(f'Make an array of 0s for MetaData/sequenceNumber') + sequenceNum = np.zeros(lon.shape, dtype=np.int32) + self.log.debug(f' sequenceNummin/max = {sequenceNum.min()} {sequenceNum.max()}') + + self.log.debug(f'Add variables to container') + container.add('sequenceNumber', sequenceNum, lon_paths) + + # Check + self.log.debug(f'container list (updated): {container.list()}') + + return container + + +add_main_functions(AdpsfcPrepbufrObsBuilder) diff --git a/dump/mapping/bufr_adpsfc_prepbufr_mapping.yaml b/dump/mapping/prepbufr_adpsfc.yaml similarity index 98% rename from dump/mapping/bufr_adpsfc_prepbufr_mapping.yaml rename to dump/mapping/prepbufr_adpsfc.yaml index 6d02794..8499669 100755 --- a/dump/mapping/bufr_adpsfc_prepbufr_mapping.yaml +++ b/dump/mapping/prepbufr_adpsfc.yaml @@ -19,23 +19,33 @@ bufr: transforms: - scale: 3600 referenceTime: "2021-08-01T00:00:00Z" + prepbufrDataLevelCategory: query: "*/CAT" + obsTimeMinusCycleTime: query: "*/DHR" + longitude: query: "*/XOB" + transforms: + - wrap: [ -180.0, 180.0 ] + latitude: query: "*/YOB" + stationIdentification: query: "*/SID" + pressure: query: "*/P___INFO/P__EVENT{1}/POB" transforms: - scale: 100 + height: query: "*/Z___INFO/Z__EVENT{1}/ZOB" type: float + stationElevation: query: "*/ELV" type: float diff --git a/dump/mapping/prepbufr_obs_builder.py b/dump/mapping/prepbufr_obs_builder.py new file mode 100644 index 0000000..9eceead --- /dev/null +++ b/dump/mapping/prepbufr_obs_builder.py @@ -0,0 +1,81 @@ +#!/usr/bin/env python3 + +import os +import re +import numpy as np +from pathlib import Path + +from datetime import datetime + +import bufr +from bufr.obs_builder import ObsBuilder + + +def map_path(map_file_name): + script_dir = os.path.dirname(os.path.abspath(__file__)) + return os.path.join(script_dir, map_file_name) + + +class PrepbufrObsBuilder(ObsBuilder): + def __init__(self, mapping_path, log_name=os.path.basename(__file__)): + super().__init__(mapping_path, log_name=log_name) + + def _get_reference_time(self, input_path) -> np.datetime64: + path_components = Path(input_path).parts + + dump = re.compile(r'\w+\.(?P\d{4})(?P\d{2})(?P\d{2})') + test = re.compile(r'(?P\d{4})(?P\d{2})(?P\d{2})(?P\d{2})') + + for idx, component in enumerate(reversed(path_components[:-1])): + dump_match = dump.match(component) + test_match = test.match(component) + + if dump_match: + ref_time = datetime(year=int(dump_match.group('year')), + month=int(dump_match.group('month')), + day=int(dump_match.group('day')), + hour=int(path_components[-1*(idx+1) + 1])) + break + elif test_match: + ref_time = datetime(year=int(test_match.group('year')), + month=int(test_match.group('month')), + day=int(test_match.group('day')), + hour=int(test_match.group('hour'))) + break + else: + print(f'Reference date not found in path.') + ref_time = datetime(year=2020, month=1, day=1) + + return np.datetime64(ref_time) + + def _compute_datetime(self, cycleTimeSinceEpoch, dhr): + """ + Compute dateTime using the cycleTimeSinceEpoch and Cycle Time + minus Cycle Time + + Parameters: + cycleTimeSinceEpoch: Time of cycle in Epoch Time + dhr: Observation Time Minus Cycle Time + + Returns: + Masked array of dateTime values + """ + + int64_fill_value = np.int64(0) + + dateTime = np.zeros(dhr.shape, dtype=np.int64) + for i in range(len(dateTime)): + if ma.is_masked(dhr[i]): + continue + else: + dateTime[i] = np.int64(dhr[i]*3600) + cycleTimeSinceEpoch + + dateTime = ma.array(dateTime) + dateTime = ma.masked_values(dateTime, int64_fill_value) + + return dateTime + + def _add_timestamp(self, container: bufr.DataContainer, reference_time: np.datetime64) -> np.array: + cycle_times = np.array([3600 * t for t in container.get('obsTimeMinusCycleTime')]).astype('timedelta64[s]') + time = (reference_time + cycle_times).astype('datetime64[s]').astype('int64') + container.add('timestamp', time, ['*']) diff --git a/dump/mapping/prepbufr_sfcshp.py b/dump/mapping/prepbufr_sfcshp.py new file mode 100755 index 0000000..e5656f8 --- /dev/null +++ b/dump/mapping/prepbufr_sfcshp.py @@ -0,0 +1,132 @@ +#!/usr/bin/env python3 +import calendar +import os +import numpy as np +import numpy.ma as ma + +import bufr +from bufr.obs_builder import add_main_functions +from prepbufr_obs_builder import PrepbufrObsBuilder, map_path + + +MAPPING_PATH = map_path('bufr_sfcshp_prepbufr.yaml') + + +class SfcshpPrepbufrObsBuilder(PrepbufrObsBuilder): + def __init__(self): + super().__init__(MAPPING_PATH, log_name=os.path.basename(__file__)) + + def _make_description(self): + description = super._make_description() + + description.add_variables([ + { + 'name': "MetaData/dateTime", + 'source': 'timestamp', + 'units': "seconds since 1970-01-01T00:00:00Z", + 'longName': "Observation Time" + }, + { + 'name': 'MetaData/sequenceNumber', + 'source': 'sequenceNumber', + 'units': '1', + 'longName': 'Sequence Number (Obs Subtype)', + } + ]) + + return description + + def make_obs(self, comm, input_path): + """ + Create the ioda sfcshp prepbufr observations: + - reads values + - adds sequenceNum + + Parameters + ---------- + comm: object + The communicator object (e.g., MPI) + input_path: str + The input bufr file + """ + + # Get container from mapping file first + self.log.info('Get container from bufr') + container = self.super().make_obs(comm, input_path) + + self.log.debug(f'container list (original): {container.list()}') + + self.log.debug(f'Do DateTime calculation') + self._add_timestamp(container, self._get_reference_time(input_path)) + + self.log.debug(f'Do sequenceNumber (Obs SubType) calculation') + typ = container.get('observationType') + typ_paths = container.get_paths('observationType') + t29 = container.get('obssubtype') + t29_paths = container.get_paths('obssubtype') + seqNum = self._compute_sequence_number(typ, t29) + self.log.debug(f' sequenceNum min/max = {seqNum.min()} {seqNum.max()}') + + self.log.debug(f'Do tsen and tv calculation') + tpc = container.get('temperatureEventCode') + tob = container.get('airTemperatureObsValue') + tob_paths = container.get_paths('airTemperatureObsValue') + tsen = np.full(tob.shape[0], tob.fill_value) + tsen = np.where(((tpc >= 1) & (tpc < 8)), tob, tsen) + tvo = np.full(tob.shape[0], tob.fill_value) + tvo = np.where((tpc == 8), tob, tvo) + + self.log.debug(f'Do tsen and tv QM calculations') + tobqm = container.get('airTemperatureQualityMarker') + tsenqm = np.full(tobqm.shape[0], tobqm.fill_value) + tsenqm = np.where(((tpc >= 1) & (tpc < 8)), tobqm, tsenqm) + tvoqm = np.full(tobqm.shape[0], tobqm.fill_value) + tvoqm = np.where((tpc == 8), tobqm, tvoqm) + + self.log.debug(f'Do tsen and tv ObsError calculations') + toboe = container.get('airTemperatureObsError') + tsenoe = np.full(toboe.shape[0], toboe.fill_value) + tsenoe = np.where(((tpc >= 1) & (tpc < 8)), toboe, tsenoe) + tvooe = np.full(toboe.shape[0], toboe.fill_value) + tvooe = np.where((tpc == 8), toboe, tvooe) + + self.log.debug(f'Update variables in container') + container.replace('airTemperatureObsValue', tsen) + container.replace('airTemperatureQualityMarker', tsenqm) + container.replace('airTemperatureObsError', tsenoe) + container.replace('virtualTemperatureObsValue', tvo) + container.replace('virtualTemperatureQualityMarker', tvoqm) + container.replace('virtualTemperatureObsError', tvooe) + + self.log.debug('DEBUG', f'Add variables to container') + container.add('sequenceNumber', seqNum, typ_paths) + + # Check + self.log.debug(f'container list (updated): {container.list()}') + + return container + + def _compute_sequence_number(self, typ, t29): + """ + Compute sequenceNumber + + Parameters: + typ: observation Type (obsType) + t29: data dump report type + + Returns: + Masked array of sequenceNumber values + """ + + sequenceNumber = np.zeros(typ.shape, dtype=np.int32) + for i in range(len(typ)): + if (typ[i] == 180 or typ[i] == 280): + if (t29[i] > 555 and t29[i] < 565): + sequenceNumber[i] = 0 + else: + sequenceNumber[i] = 1 + + return sequenceNumber + + +add_main_functions(SfcshpPrepbufrObsBuilder) diff --git a/dump/mapping/bufr_sfcshp_prepbufr_mapping.yaml b/dump/mapping/prepbufr_sfcshp.yaml similarity index 98% rename from dump/mapping/bufr_sfcshp_prepbufr_mapping.yaml rename to dump/mapping/prepbufr_sfcshp.yaml index e667866..77c7a36 100755 --- a/dump/mapping/bufr_sfcshp_prepbufr_mapping.yaml +++ b/dump/mapping/prepbufr_sfcshp.yaml @@ -13,18 +13,15 @@ bufr: query: "*/TYP" # MetaData - timestamp: - timeoffset: - timeOffset: "*/DHR" - transforms: - - scale: 3600 - referenceTime: "2021-08-01T00:00:00Z" prepbufrDataLevelCategory: query: "*/CAT" obsTimeMinusCycleTime: query: "*/DHR" longitude: query: "*/XOB" + transforms: + - wrap: [ -180.0, 180.0 ] + latitude: query: "*/YOB" stationIdentification: diff --git a/ush/test/config/bufr_bufr4backend_ahicsr.yaml b/ush/test/config/bufr_bufr4backend_ahicsr.yaml new file mode 100644 index 0000000..0e56a76 --- /dev/null +++ b/ush/test/config/bufr_bufr4backend_ahicsr.yaml @@ -0,0 +1,39 @@ +time window: + begin: "2021-07-31T21:00:00Z" + end: "2021-08-01T03:00:00Z" + bound to include: begin + +observations: +- obs space: + name: "ahicsr-h8" + simulated variables: [cloudAmount] + obsdatain: + engine: + type: bufr + obsfile: "./testinput/2021080100/gdas.t00z.ahicsr.tm00.bufr_d" + mapping file: "./bufr_ahicsr.yaml" + category: ["h8"] + cache categories: # optional + - ["h8"] + - ["h9"] + obsdataout: + engine: + type: H5File + obsfile: "./testoutput/2021080100/bufr4backend/gdas.t00z.ahicsr_h8.tm00.nc" + +- obs space: + name: "ahicsr-h9" + simulated variables: [cloudAmount] + obsdatain: + engine: + type: bufr + obsfile: "./testinput/2021080100/gdas.t00z.ahicsr.tm00.bufr_d" + mapping file: "./bufr_ahicsr.yaml" + category: ["h9"] + cache categories: # optional + - ["h8"] + - ["h9"] + obsdataout: + engine: + type: H5File + obsfile: "./testoutput/2021080100/bufr4backend/gdas.t00z.ahicsr_h9.tm00.nc" diff --git a/ush/test/config/bufr_bufr4backend_satwnd_amv_ahi.yaml b/ush/test/config/bufr_bufr4backend_satwnd_amv_ahi.yaml index 56af062..49bab9a 100644 --- a/ush/test/config/bufr_bufr4backend_satwnd_amv_ahi.yaml +++ b/ush/test/config/bufr_bufr4backend_satwnd_amv_ahi.yaml @@ -10,7 +10,7 @@ observations: engine: type: bufr obsfile: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping file: "./bufr_satwnd_amv_ahi_mapping.yaml" + mapping file: "./bufr_satwnd_amv_ahi.yaml" category: ["h8"] cache categories: # optional - ["h8"] @@ -27,7 +27,7 @@ observations: engine: type: bufr obsfile: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping file: "./bufr_satwnd_amv_ahi_mapping.yaml" + mapping file: "./bufr_satwnd_amv_ahi.yaml" category: ["h9"] cache categories: # optional - ["h8"] diff --git a/ush/test/config/bufr_bufr4backend_satwnd_amv_avhrr.yaml b/ush/test/config/bufr_bufr4backend_satwnd_amv_avhrr.yaml index f721553..dd764e5 100644 --- a/ush/test/config/bufr_bufr4backend_satwnd_amv_avhrr.yaml +++ b/ush/test/config/bufr_bufr4backend_satwnd_amv_avhrr.yaml @@ -10,7 +10,7 @@ observations: engine: type: bufr obsfile: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping file: "./bufr_satwnd_amv_avhrr_mapping.yaml" + mapping file: "./bufr_satwnd_amv_avhrr.yaml" category: ["metop-1a"] cache categories: # optional - ["metop-1a"] @@ -31,7 +31,7 @@ observations: engine: type: bufr obsfile: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping file: "./bufr_satwnd_amv_avhrr_mapping.yaml" + mapping file: "./bufr_satwnd_amv_avhrr.yaml" category: ["metop-1b"] cache categories: # optional - ["metop-1a"] @@ -52,7 +52,7 @@ observations: engine: type: bufr obsfile: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping file: "./bufr_satwnd_amv_avhrr_mapping.yaml" + mapping file: "./bufr_satwnd_amv_avhrr.yaml" category: ["metop-1c"] cache categories: # optional - ["metop-1a"] @@ -73,7 +73,7 @@ observations: engine: type: bufr obsfile: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping file: "./bufr_satwnd_amv_avhrr_mapping.yaml" + mapping file: "./bufr_satwnd_amv_avhrr.yaml" category: ["noaa-15"] cache categories: # optional - ["metop-1a"] @@ -94,7 +94,7 @@ observations: engine: type: bufr obsfile: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping file: "./bufr_satwnd_amv_avhrr_mapping.yaml" + mapping file: "./bufr_satwnd_amv_avhrr.yaml" category: ["noaa-18"] cache categories: # optional - ["metop-1a"] @@ -115,7 +115,7 @@ observations: engine: type: bufr obsfile: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping file: "./bufr_satwnd_amv_avhrr_mapping.yaml" + mapping file: "./bufr_satwnd_amv_avhrr.yaml" category: ["noaa-19"] cache categories: # optional - ["metop-1a"] diff --git a/ush/test/config/bufr_bufr4backend_satwnd_amv_leogeo.yaml b/ush/test/config/bufr_bufr4backend_satwnd_amv_leogeo.yaml index ca40eb7..c594a32 100644 --- a/ush/test/config/bufr_bufr4backend_satwnd_amv_leogeo.yaml +++ b/ush/test/config/bufr_bufr4backend_satwnd_amv_leogeo.yaml @@ -10,7 +10,7 @@ observations: engine: type: bufr obsfile: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping file: "./bufr_satwnd_amv_leogeo_mapping.yaml" + mapping file: "./bufr_satwnd_amv_leogeo.yaml" category: ["multi"] cache categories: # optional - ["multi"] diff --git a/ush/test/config/bufr_bufr4backend_satwnd_amv_modis.yaml b/ush/test/config/bufr_bufr4backend_satwnd_amv_modis.yaml index 1f3d24f..23cdfd0 100644 --- a/ush/test/config/bufr_bufr4backend_satwnd_amv_modis.yaml +++ b/ush/test/config/bufr_bufr4backend_satwnd_amv_modis.yaml @@ -10,7 +10,7 @@ observations: engine: type: bufr obsfile: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping file: "./bufr_satwnd_amv_modis_mapping.yaml" + mapping file: "./bufr_satwnd_amv_modis.yaml" category: ["terra"] cache categories: # optional - ["terra"] @@ -27,7 +27,7 @@ observations: engine: type: bufr obsfile: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping file: "./bufr_satwnd_amv_modis_mapping.yaml" + mapping file: "./bufr_satwnd_amv_modis.yaml" category: ["aqua"] cache categories: # optional - ["terra"] diff --git a/ush/test/config/bufr_bufr4backend_satwnd_amv_seviri.yaml b/ush/test/config/bufr_bufr4backend_satwnd_amv_seviri.yaml index 40db3b2..3664760 100644 --- a/ush/test/config/bufr_bufr4backend_satwnd_amv_seviri.yaml +++ b/ush/test/config/bufr_bufr4backend_satwnd_amv_seviri.yaml @@ -10,7 +10,7 @@ observations: engine: type: bufr obsfile: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping file: "./bufr_satwnd_amv_seviri_mapping.yaml" + mapping file: "./bufr_satwnd_amv_seviri.yaml" category: ["m8"] cache categories: # optional - ["m8"] @@ -29,7 +29,7 @@ observations: engine: type: bufr obsfile: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping file: "./bufr_satwnd_amv_seviri_mapping.yaml" + mapping file: "./bufr_satwnd_amv_seviri.yaml" category: ["m9"] cache categories: # optional - ["m8"] @@ -48,7 +48,7 @@ observations: engine: type: bufr obsfile: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping file: "./bufr_satwnd_amv_seviri_mapping.yaml" + mapping file: "./bufr_satwnd_amv_seviri.yaml" category: ["m10"] cache categories: # optional - ["m8"] @@ -67,7 +67,7 @@ observations: engine: type: bufr obsfile: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping file: "./bufr_satwnd_amv_seviri_mapping.yaml" + mapping file: "./bufr_satwnd_amv_seviri.yaml" category: ["m11"] cache categories: # optional - ["m8"] diff --git a/ush/test/config/bufr_bufr4backend_satwnd_amv_viirs.yaml b/ush/test/config/bufr_bufr4backend_satwnd_amv_viirs.yaml index 5f53e5e..ca0f441 100644 --- a/ush/test/config/bufr_bufr4backend_satwnd_amv_viirs.yaml +++ b/ush/test/config/bufr_bufr4backend_satwnd_amv_viirs.yaml @@ -10,7 +10,7 @@ observations: engine: type: bufr obsfile: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping file: "./bufr_satwnd_amv_viirs_mapping.yaml" + mapping file: "./bufr_satwnd_amv_viirs.yaml" category: ["npp"] cache categories: # optional - ["npp"] @@ -28,7 +28,7 @@ observations: engine: type: bufr obsfile: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping file: "./bufr_satwnd_amv_viirs_mapping.yaml" + mapping file: "./bufr_satwnd_amv_viirs.yaml" category: ["n20"] cache categories: # optional - ["npp"] @@ -46,7 +46,7 @@ observations: engine: type: bufr obsfile: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping file: "./bufr_satwnd_amv_viirs_mapping.yaml" + mapping file: "./bufr_satwnd_amv_viirs.yaml" category: ["n21"] cache categories: # optional - ["npp"] diff --git a/ush/test/config/bufr_bufr4backend_sfcshp_prepbufr.yaml b/ush/test/config/bufr_bufr4backend_sfcshp_prepbufr.yaml index ee0b0c9..497a46c 100644 --- a/ush/test/config/bufr_bufr4backend_sfcshp_prepbufr.yaml +++ b/ush/test/config/bufr_bufr4backend_sfcshp_prepbufr.yaml @@ -11,7 +11,7 @@ observations: engine: type: bufr obsfile: "./testinput/2021080100/gdas.t00z.prepbufr" - mapping file: "./bufr_sfcshp_prepbufr_mapping.yaml" + mapping file: "./bufr_sfcshp_prepbufr.yaml" obsdataout: engine: type: H5File diff --git a/ush/test/config/bufr_bufr4backend_sfcshp_prepbufr_mpi4.yaml b/ush/test/config/bufr_bufr4backend_sfcshp_prepbufr_mpi4.yaml index 9f042e2..d4b5e06 100644 --- a/ush/test/config/bufr_bufr4backend_sfcshp_prepbufr_mpi4.yaml +++ b/ush/test/config/bufr_bufr4backend_sfcshp_prepbufr_mpi4.yaml @@ -11,7 +11,7 @@ observations: engine: type: bufr obsfile: "./testinput/2021080100/gdas.t00z.prepbufr" - mapping file: "./bufr_sfcshp_prepbufr_mapping.yaml" + mapping file: "./bufr_sfcshp_prepbufr.yaml" obsdataout: engine: type: H5File diff --git a/ush/test/config/bufr_script4backend_ahicsr.yaml b/ush/test/config/bufr_script4backend_ahicsr.yaml new file mode 100644 index 0000000..63ccd76 --- /dev/null +++ b/ush/test/config/bufr_script4backend_ahicsr.yaml @@ -0,0 +1,41 @@ +time window: + begin: "2021-07-31T21:00:00Z" + end: "2021-08-01T03:00:00Z" + bound to include: begin + +observations: +- obs space: + name: "ahicsr-h8" + observed variables: [cloudAmount] + derived variables: [cloudAmount] + simulated variables: [cloudAmount] + obsdatain: + engine: + type: script + script file: "bufr_ahicsr.py" + args: + input_path: "./testinput/2021080100/gdas.t00z.ahicsr.tm00.bufr_d" + mapping_path: "./bufr_ahicsr.yaml" + category: "h8" + obsdataout: + engine: + type: H5File + obsfile: "./testoutput/2021080100/script4backend/gdas.t00z.ahicsr_h8.tm00.nc" + +- obs space: + name: "ahicsr-h9" + observed variables: [cloudAmount] + derived variables: [cloudAmount] + simulated variables: [cloudAmount] + obsdatain: + engine: + type: script + script file: "bufr_ahicsr.py" + args: + input_path: "./testinput/2021080100/gdas.t00z.ahicsr.tm00.bufr_d" + mapping_path: "./bufr_ahicsr.yaml" + category: "h9" + obsdataout: + engine: + type: H5File + obsfile: "./testoutput/2021080100/script4backend/gdas.t00z.ahicsr_h9.tm00.nc" diff --git a/ush/test/config/bufr_script4backend_satwnd_amv_ahi.yaml b/ush/test/config/bufr_script4backend_satwnd_amv_ahi.yaml index 9f1a1d8..86ef32b 100644 --- a/ush/test/config/bufr_script4backend_satwnd_amv_ahi.yaml +++ b/ush/test/config/bufr_script4backend_satwnd_amv_ahi.yaml @@ -14,7 +14,7 @@ observations: script file: "bufr_satwnd_amv_ahi.py" args: input_path: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping_path: "./bufr_satwnd_amv_ahi_mapping.yaml" + mapping_path: "./bufr_satwnd_amv_ahi.yaml" category: "h8" obsdataout: engine: @@ -32,7 +32,7 @@ observations: script file: "bufr_satwnd_amv_ahi.py" args: input_path: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping_path: "./bufr_satwnd_amv_ahi_mapping.yaml" + mapping_path: "./bufr_satwnd_amv_ahi.yaml" category: "h9" obsdataout: engine: diff --git a/ush/test/config/bufr_script4backend_satwnd_amv_avhrr.yaml b/ush/test/config/bufr_script4backend_satwnd_amv_avhrr.yaml index c8c2f0f..61a66de 100644 --- a/ush/test/config/bufr_script4backend_satwnd_amv_avhrr.yaml +++ b/ush/test/config/bufr_script4backend_satwnd_amv_avhrr.yaml @@ -14,7 +14,7 @@ observations: script file: "bufr_satwnd_amv_avhrr.py" args: input_path: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping_path: "./bufr_satwnd_amv_avhrr_mapping.yaml" + mapping_path: "./bufr_satwnd_amv_avhrr.yaml" category: "metop-1a" obsdataout: engine: @@ -32,7 +32,7 @@ observations: script file: "bufr_satwnd_amv_avhrr.py" args: input_path: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping_path: "./bufr_satwnd_amv_avhrr_mapping.yaml" + mapping_path: "./bufr_satwnd_amv_avhrr.yaml" category: "metop-1b" obsdataout: engine: @@ -50,7 +50,7 @@ observations: script file: "bufr_satwnd_amv_avhrr.py" args: input_path: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping_path: "./bufr_satwnd_amv_avhrr_mapping.yaml" + mapping_path: "./bufr_satwnd_amv_avhrr.yaml" category: "metop-1c" obsdataout: engine: @@ -68,7 +68,7 @@ observations: script file: "bufr_satwnd_amv_avhrr.py" args: input_path: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping_path: "./bufr_satwnd_amv_avhrr_mapping.yaml" + mapping_path: "./bufr_satwnd_amv_avhrr.yaml" category: "noaa-15" obsdataout: engine: @@ -86,7 +86,7 @@ observations: script file: "bufr_satwnd_amv_avhrr.py" args: input_path: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping_path: "./bufr_satwnd_amv_avhrr_mapping.yaml" + mapping_path: "./bufr_satwnd_amv_avhrr.yaml" category: "noaa-18" obsdataout: engine: @@ -104,7 +104,7 @@ observations: script file: "bufr_satwnd_amv_avhrr.py" args: input_path: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping_path: "./bufr_satwnd_amv_avhrr_mapping.yaml" + mapping_path: "./bufr_satwnd_amv_avhrr.yaml" category: "noaa-19" obsdataout: engine: diff --git a/ush/test/config/bufr_script4backend_satwnd_amv_leogeo.yaml b/ush/test/config/bufr_script4backend_satwnd_amv_leogeo.yaml index cef8af0..1baa637 100644 --- a/ush/test/config/bufr_script4backend_satwnd_amv_leogeo.yaml +++ b/ush/test/config/bufr_script4backend_satwnd_amv_leogeo.yaml @@ -14,7 +14,7 @@ observations: script file: "bufr_satwnd_amv_leogeo.py" args: input_path: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping_path: "./bufr_satwnd_amv_leogeo_mapping.yaml" + mapping_path: "./bufr_satwnd_amv_leogeo.yaml" category: "multi" obsdataout: engine: diff --git a/ush/test/config/bufr_script4backend_satwnd_amv_modis.yaml b/ush/test/config/bufr_script4backend_satwnd_amv_modis.yaml index 8c25f6a..18bcc84 100644 --- a/ush/test/config/bufr_script4backend_satwnd_amv_modis.yaml +++ b/ush/test/config/bufr_script4backend_satwnd_amv_modis.yaml @@ -14,7 +14,7 @@ observations: script file: "bufr_satwnd_amv_modis.py" args: input_path: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping_path: "./bufr_satwnd_amv_modis_mapping.yaml" + mapping_path: "./bufr_satwnd_amv_modis.yaml" category: "terra" obsdataout: engine: @@ -32,7 +32,7 @@ observations: script file: "bufr_satwnd_amv_modis.py" args: input_path: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping_path: "./bufr_satwnd_amv_modis_mapping.yaml" + mapping_path: "./bufr_satwnd_amv_modis.yaml" category: "aqua" obsdataout: engine: diff --git a/ush/test/config/bufr_script4backend_satwnd_amv_seviri.yaml b/ush/test/config/bufr_script4backend_satwnd_amv_seviri.yaml index 9f40ea4..9c744a1 100644 --- a/ush/test/config/bufr_script4backend_satwnd_amv_seviri.yaml +++ b/ush/test/config/bufr_script4backend_satwnd_amv_seviri.yaml @@ -14,7 +14,7 @@ observations: script file: "bufr_satwnd_amv_seviri.py" args: input_path: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping_path: "./bufr_satwnd_amv_seviri_mapping.yaml" + mapping_path: "./bufr_satwnd_amv_seviri.yaml" category: "m8" obsdataout: engine: @@ -32,7 +32,7 @@ observations: script file: "bufr_satwnd_amv_seviri.py" args: input_path: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping_path: "./bufr_satwnd_amv_seviri_mapping.yaml" + mapping_path: "./bufr_satwnd_amv_seviri.yaml" category: "m9" obsdataout: engine: @@ -50,7 +50,7 @@ observations: script file: "bufr_satwnd_amv_seviri.py" args: input_path: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping_path: "./bufr_satwnd_amv_seviri_mapping.yaml" + mapping_path: "./bufr_satwnd_amv_seviri.yaml" category: "m10" obsdataout: engine: @@ -68,7 +68,7 @@ observations: script file: "bufr_satwnd_amv_seviri.py" args: input_path: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping_path: "./bufr_satwnd_amv_seviri_mapping.yaml" + mapping_path: "./bufr_satwnd_amv_seviri.yaml" category: "m11" obsdataout: engine: diff --git a/ush/test/config/bufr_script4backend_satwnd_amv_viirs.yaml b/ush/test/config/bufr_script4backend_satwnd_amv_viirs.yaml index ddb1958..72fcc42 100644 --- a/ush/test/config/bufr_script4backend_satwnd_amv_viirs.yaml +++ b/ush/test/config/bufr_script4backend_satwnd_amv_viirs.yaml @@ -14,7 +14,7 @@ observations: script file: "bufr_satwnd_amv_viirs.py" args: input_path: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping_path: "./bufr_satwnd_amv_viirs_mapping.yaml" + mapping_path: "./bufr_satwnd_amv_viirs.yaml" category: "npp" obsdataout: engine: @@ -32,7 +32,7 @@ observations: script file: "bufr_satwnd_amv_viirs.py" args: input_path: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping_path: "./bufr_satwnd_amv_viirs_mapping.yaml" + mapping_path: "./bufr_satwnd_amv_viirs.yaml" category: "n20" obsdataout: engine: @@ -50,7 +50,7 @@ observations: script file: "bufr_satwnd_amv_viirs.py" args: input_path: "./testinput/2021080100/gdas.t00z.satwnd.tm00.bufr_d" - mapping_path: "./bufr_satwnd_amv_viirs_mapping.yaml" + mapping_path: "./bufr_satwnd_amv_viirs.yaml" category: "n21" obsdataout: engine: diff --git a/ush/test/config/bufr_script4backend_sfcshp_prepbufr.yaml b/ush/test/config/bufr_script4backend_sfcshp_prepbufr.yaml index aabab8d..2009c98 100644 --- a/ush/test/config/bufr_script4backend_sfcshp_prepbufr.yaml +++ b/ush/test/config/bufr_script4backend_sfcshp_prepbufr.yaml @@ -15,7 +15,7 @@ observations: script file: "bufr_sfcshp_prepbufr.py" args: input_path: "testinput/2021080100/gdas.t00z.prepbufr" - mapping_path: "./bufr_sfcshp_prepbufr_mapping.yaml" + mapping_path: "./bufr_sfcshp_prepbufr.yaml" cycle_time: "2021080100" obsdataout: engine: diff --git a/ush/test/config/bufr_script4backend_sfcshp_prepbufr_mpi4.yaml b/ush/test/config/bufr_script4backend_sfcshp_prepbufr_mpi4.yaml index d60c53a..3a938bf 100644 --- a/ush/test/config/bufr_script4backend_sfcshp_prepbufr_mpi4.yaml +++ b/ush/test/config/bufr_script4backend_sfcshp_prepbufr_mpi4.yaml @@ -15,7 +15,7 @@ observations: script file: "bufr_sfcshp_prepbufr.py" args: input_path: "testinput/2021080100/gdas.t00z.prepbufr" - mapping_path: "./bufr_sfcshp_prepbufr_mapping.yaml" + mapping_path: "./bufr_sfcshp_prepbufr.yaml" cycle_time: "2021080100" obsdataout: engine: