Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add the "state of ground" to bring in the extra snow depth reports #1368

Draft
wants to merge 20 commits into
base: develop
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
20 commits
Select commit Hold shift + click to select a range
1c8096f
Add state of ground to increase the snow depth reports.
jiaruidong2017 Nov 9, 2024
365675d
Update the defined index with zero snow depth.
jiaruidong2017 Nov 13, 2024
aee4dc2
Merge branch 'develop' into feature/snow-sogr
jiaruidong2017 Nov 13, 2024
0c02522
Move python code to parm/gdas/snow/ directory to aviod extra copy.
jiaruidong2017 Nov 13, 2024
e57ac67
Modif to only replace the missing/filled snod values.
jiaruidong2017 Nov 13, 2024
4376049
Merge branch 'develop' into feature/snow-sogr
jiaruidong2017 Nov 13, 2024
244693c
Address reviewer's comments.
jiaruidong2017 Nov 14, 2024
76a5e94
Fix py norms error.
jiaruidong2017 Nov 14, 2024
2a72e4d
Add the mapping_path indicating the script backend requirements.
jiaruidong2017 Nov 15, 2024
a2d87fd
Merge branch 'develop' into feature/snow-sogr
jiaruidong2017 Nov 15, 2024
fbd1bb9
Add comments
jiaruidong2017 Nov 15, 2024
db46314
Merge branch 'develop' into feature/snow-sogr
jiaruidong2017 Nov 22, 2024
91f160e
Update the missing value for elevation from bufr file.
jiaruidong2017 Nov 22, 2024
dc8d9b8
Merge branch 'develop' into feature/snow-sogr
jiaruidong2017 Dec 20, 2024
d7bf5d7
Merge branch 'develop' into feature/snow-sogr
jiaruidong2017 Dec 24, 2024
3c1dbf7
Update the scripts and sorc/bufr-query.
jiaruidong2017 Dec 24, 2024
55b54cd
Merge branch 'develop' into feature/snow-sogr
jiaruidong2017 Jan 13, 2025
5cd21ab
Remove the negative snow depth values.
jiaruidong2017 Jan 18, 2025
47bfb19
Merge branch 'develop' into feature/snow-sogr
jiaruidong2017 Jan 18, 2025
bb85eb6
Fix a logic error
jiaruidong2017 Jan 18, 2025
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions parm/snow/jcb-base.yaml.j2
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,7 @@ crtm_coefficient_path: "{{ DATA }}/crtm/"

# Naming conventions for observational files
snow_obsdataroot_path: "{{COMIN_OBS}}"
snow_script_path: "{{snow_script_path}}"
snow_obsdatain_path: "{{snow_obsdatain_path}}"
snow_obsdatain_prefix: "{{OPREFIX}}"
snow_obsdatain_suffix: ".tm00.bufr_d"
Expand Down
18 changes: 11 additions & 7 deletions parm/snow/obs/config/bufr_sfcsno_mapping.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -17,20 +17,17 @@ bufr:
query: "[*/CLON, */CLONH]"
stationIdentification:
query: "*/RPID"

stationElevation:
query: "[*/SELV, */HSMSL]"
type: float

# ObsValue
totalSnowDepth:
query: "[*/SNWSQ1/TOSD, */MTRMSC/TOSD, */STGDSNDM/TOSD]"
transforms:
- scale: 1000.0
Comment on lines 27 to 28
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Based on IODA data convention, the unit of totalSnowDepth is m. We should remove the conversion of it from meter to millimeter.

Copy link
Collaborator Author

@jiaruidong2017 jiaruidong2017 Jan 13, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as above to keep use the conversion.

filters:
- bounding:
variable: totalSnowDepth
lowerBound: 0
upperBound: 10000000
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For my understanding:

  • is this introducing a 10 m limit on snow depth obs?
  • if they are above 10 m, are they discarded, or set to 10.?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I previously put this large number (10000m) here for removing the missing values. Here, we don't need this any more because we need to set the missing snod values to 0 when sogr satisfies the defined conditions.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We might still have missing values though wouldn't we?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The missing values will be removed after applying sogr conditions.

groundState:
query: "[*/GRDSQ1/SOGR, */STGDSNDM/SOGR]"

encoder:
variables:
Expand Down Expand Up @@ -65,11 +62,18 @@ encoder:
coordinates: "longitude latitude"
source: variables/stationIdentification
longName: "Identification of Observing Location"
units: "m"
units: "index"
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jiaruidong2017 The unit of stationIdentification is not defined (empty) in the IODA convention table, and it is left for users to decide. Have you discussed the unit for stationIdentification with JCSDA (Ben Ruston)?
Let's leave it as "index" for now and I will open an issue in IODA repository and check if we can use "index" as unit for this variable.


# ObsValue
- name: "ObsValue/totalSnowDepth"
coordinates: "longitude latitude"
source: variables/totalSnowDepth
longName: "Total Snow Depth"
units: "mm"
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copy link
Collaborator Author

@jiaruidong2017 jiaruidong2017 Jan 13, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We have to keep using the mm here in the observation files to match the unit used in the UFS model for the snow depth.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jiaruidong2017 Is it possible to read in the totalSnowDepth in m unit from the IODA/NetCDF file and then convert it to mm in the UFS model?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is possible, but will require additional efforts. We can discuss this issue with the physics team.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jiaruidong2017 I guess only a one-line change is necessary for UFS. You just need to convert totalSnowDepth from m to mm where it is read from IODA obs file.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This issue (mismatch between units proscribed by JEDI and those used in the UFS) is bigger than this PR - as all other snow depth IODA files are in mm too. If we change these observations to mm, we'll need to change the others. I suggest that we leave this in mm for now, and create a separate issue to address the unit mismatch (by introducing a unit transform somewhere).

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed. Let's keep it as "mm" for now.


- name: "ObsValue/groundState"
coordinates: "longitude latitude"
source: variables/groundState
longName: "STATE OF THE GROUND"
units: "index"

2 changes: 1 addition & 1 deletion sorc/bufr-query
Submodule bufr-query updated 39 files
+4 −2 CMakeLists.txt
+5 −5 cmake/Targets.cmake
+4 −2 core/CMakeLists.txt
+8 −0 core/include/bufr/BufrParser.h
+25 −2 core/include/bufr/DataContainer.h
+523 −1 core/include/bufr/DataObject.h
+5 −7 core/include/bufr/DataProvider.h
+8 −3 core/include/bufr/File.h
+3 −0 core/include/bufr/NcepDataProvider.h
+3 −0 core/include/bufr/WmoDataProvider.h
+31 −3 core/include/bufr/encoders/Description.h
+66 −1 core/src/bufr/BufrReader/BufrParser.cpp
+51 −11 core/src/bufr/BufrReader/Query/DataProvider/DataProvider.cpp
+8 −0 core/src/bufr/BufrReader/Query/DataProvider/NcepDataProvider.cpp
+7 −0 core/src/bufr/BufrReader/Query/DataProvider/WmoDataProvider.cpp
+12 −5 core/src/bufr/BufrReader/Query/File.cpp
+56 −0 core/src/bufr/DataContainer.cpp
+96 −8 core/src/encoders/Description.cpp
+11 −2 core/src/encoders/netcdf/Encoder.cpp
+51 −0 docs/python_api.rst
+1 −0 docs/uml/BUFR_BigPicture.puml
+25 −28 python/CMakeLists.txt
+17 −10 python/DataObjectFunctions.cpp
+4 −4 python/DataObjectFunctions.h
+4 −0 python/py_bufr.cpp
+17 −1 python/py_data_container.cpp
+80 −1 python/py_encoder_description.cpp
+2 −1 python/py_file.cpp
+49 −0 python/py_mpi.cpp
+23 −0 python/py_mpi.h
+15 −2 python/py_parser.cpp
+12 −8 test/CMakeLists.txt
+2 −2 test/testinput/bufrtest_fieldname_validation.py
+88 −0 test/testinput/bufrtest_python_mpi_test.py
+65 −17 test/testinput/bufrtest_python_test.py
+3 −0 tools/bufr2netcdf/CMakeLists.txt
+140 −44 tools/bufr2netcdf/bufr2netcdf.cpp
+1 −1 tools/build_scripts/bufr_query_cpp_lint.py
+11 −3 tools/show_queries/QueryPrinter.cpp
190 changes: 190 additions & 0 deletions ush/ioda/bufr2ioda/bufr_sfcsno.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,190 @@
#!/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)

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)

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}')