Skip to content

Commit e17815f

Browse files
ScottWalesScottWales
committed
Working grid_stat and plotting
1 parent 3f46b55 commit e17815f

16 files changed

Lines changed: 411 additions & 53 deletions

File tree

‎docs/source/reference/workflow/grid_stat.rst‎

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,35 @@
1+
Comparing against Gridded Data
2+
==============================
3+
4+
METplus can perform comparison against gridded datasets using the ``grid_stat``
5+
tool. Enable ``grid_stat`` comparisons by adding to your ``rose-suite.conf``::
6+
7+
RUN_METPLUS_GRID_STAT = True
8+
9+
Observations have to be configured for your site so that METplus knows the
10+
correct source data paths.
11+
12+
Select which observations to run using e.g.::
13+
14+
METPLUS_GRID_STAT_OBS = ["GPM"]
15+
16+
Output
17+
------
18+
19+
The workflow creates spatial plots of the model field, observation field and
20+
difference between model and observation. ``grid_stat`` handles regridding both
21+
data sources to a common grid.
22+
23+
Observation Types
24+
-----------------
25+
26+
GPM
27+
^^^
28+
29+
https://gpm.nasa.gov/data/imerg
30+
31+
Compares model field ``stratiform_rainfall_flux`` against observation field
32+
``precipitation_flux``.
33+
34+
Available sites:
35+
* NCI

‎docs/source/reference/workflow/index.rst‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,3 +7,4 @@ Reference information for the CSET workflow
77
:maxdepth: 1
88

99
observations
10+
grid_stat

‎src/CSET/cset_workflow/app/grid_stat_plot/bin/grid_stat_plot.py‎

Lines changed: 197 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,197 @@
1+
#!/usr/bin/env python3
2+
3+
"""Create CSET plots from MET grid_stat output."""
4+
5+
import argparse
6+
import json
7+
import logging
8+
import os
9+
from pathlib import Path
10+
from typing import Iterable
11+
12+
import iris.coords
13+
import iris.cube
14+
15+
from CSET.operators.plot import spatial_pcolormesh_plot
16+
17+
log = logging.getLogger(__name__)
18+
19+
20+
def standardise_names(cube: iris.cube.Cube) -> iris.cube.Cube:
21+
"""
22+
Convert the name MET gives fields to something more standard.
23+
24+
From a name like
25+
{type}_{var0}_{levels0}(_{var1}_{levels1})?_{region}
26+
27+
Sets:
28+
var_name to '{type}_{var0}'
29+
long_name to '{var0}_{long_type}'
30+
standard_name to '{var0}' (for OBS and FCST fields only)
31+
32+
Modifies the cube in-place and returns it
33+
"""
34+
assert cube.var_name is not None
35+
name = cube.var_name
36+
levels = cube.attributes["level"].split(" and ")
37+
mask_region = cube.attributes["masking_region"]
38+
type = name.split("_")[0]
39+
40+
level_names = [lev.replace("*", "all").replace(",", "_") for lev in levels]
41+
42+
name_middle = name[len(type) + 1 : -(len(mask_region) + 1)]
43+
44+
var0 = name_middle[: name_middle.index(level_names[0]) - 1]
45+
46+
cube.var_name = f"{type}_{var0}_{mask_region}"
47+
cube.attributes["type"] = type
48+
49+
if type == "OBS":
50+
cube.long_name = f"{var0}_observation"
51+
cube.standard_name = var0
52+
cube.attributes["verification_type"] = "Observation Value"
53+
elif type == "FCST":
54+
cube.long_name = f"{var0}_forecast"
55+
cube.standard_name = var0
56+
cube.attributes["verification_type"] = "Forecast Value"
57+
elif type == "DIFF":
58+
cube.long_name = f"{var0}_difference"
59+
cube.attributes["model_var"] = var0
60+
cube.attributes["verification_type"] = cube.attributes["Difference"]
61+
62+
return cube
63+
64+
65+
def postproc_cube(cube: iris.cube.Cube, field: str, filename: str):
66+
"""
67+
Prepare the MET data for processing.
68+
69+
Use this as a callback when loading Iris cubes from grid_stat output.
70+
71+
Fixes the names and sets up time coordinates properly
72+
"""
73+
cube = standardise_names(cube)
74+
75+
# Grab timestamps
76+
init_time_ut = int(cube.attributes.pop("init_time_ut"))
77+
valid_time_ut = int(cube.attributes.pop("valid_time_ut"))
78+
79+
# Clean up attributes that might not match between different files
80+
cube.attributes.pop("init_time")
81+
cube.attributes.pop("valid_time")
82+
cube.attributes.pop("level")
83+
cube.attributes.pop("name")
84+
cube.attributes.pop("FileOrigins")
85+
86+
cube.add_aux_coord(
87+
iris.coords.DimCoord(
88+
init_time_ut,
89+
standard_name="forecast_reference_time",
90+
units="seconds since 1970-01-01 00:00Z",
91+
)
92+
)
93+
cube.add_aux_coord(
94+
iris.coords.DimCoord(
95+
valid_time_ut, standard_name="time", units="seconds since 1970-01-01 00:00Z"
96+
)
97+
)
98+
99+
if "invalid_units" in cube.attributes:
100+
cube.units = cube.attributes["invalid_units"].split(" and ")[0]
101+
102+
103+
def load_grid_stat(paths: Iterable[Path]) -> iris.cube.CubeList:
104+
"""Load processed grid_stat netcdf output as iris cubes."""
105+
cubes = iris.cube.CubeList()
106+
for path in paths:
107+
cubes.extend(iris.load(path, callback=postproc_cube))
108+
return cubes.merge()
109+
110+
111+
def plot_grid_stat(
112+
cube: iris.cube.Cube,
113+
webdir: Path,
114+
*,
115+
style_file: Path | str | None = None,
116+
plot_resolution: int | None = None,
117+
):
118+
"""
119+
Create plots from grid_stat output.
120+
121+
Parameters
122+
----------
123+
cube:
124+
Input data cube
125+
webdir:
126+
Base web path for this cycle
127+
style_file:
128+
Colorbar definition JSON file
129+
plot_resolution:
130+
Plot resolution in pixels per inch
131+
"""
132+
model = cube.attributes["model"]
133+
obstype = cube.attributes["obtype"]
134+
assert cube.var_name is not None
135+
outdir = webdir / "grid_stat" / f"{model}_vs_{obstype}" / cube.var_name
136+
outdir.mkdir(exist_ok=True, parents=True)
137+
os.chdir(outdir)
138+
139+
log.info("writing %s to %s", cube.name(), outdir)
140+
141+
# Set up CSET metadata
142+
meta = {
143+
"category": "Gridded Verification",
144+
"title": f"{cube.long_name}",
145+
"case_date": cube.coord("forecast_reference_time")
146+
.as_string_arrays(fmt="%Y%m%dT%H%MZ")
147+
.points[0],
148+
"MODEL_NAME": model,
149+
"OBSERVATION_TYPE": obstype,
150+
"GRID_STAT_TYPE": cube.attributes["verification_type"],
151+
"SUBAREA_EXTENT": cube.attributes["masking_region"],
152+
}
153+
154+
if cube.attributes["type"] == "OBS":
155+
meta["title"] = f"{obstype} {cube.long_name}"
156+
elif cube.attributes["type"] == "FCST":
157+
meta["title"] = f"{model} {cube.long_name}"
158+
elif cube.attributes["type"] == "DIFF":
159+
meta["title"] = f"{model} vs {obstype} {cube.long_name}"
160+
161+
with open("meta.json", "wt") as f:
162+
json.dump(meta, f)
163+
164+
cube.coord("latitude").guess_bounds()
165+
cube.coord("longitude").guess_bounds()
166+
spatial_pcolormesh_plot(cube, f"grid_stat_DarwinCTL_vs_GPM_{cube.var_name}.png")
167+
168+
169+
def main():
170+
"""CLI entry point."""
171+
parser = argparse.ArgumentParser(description=__doc__)
172+
parser.add_argument(
173+
"--web-dir", help="directory to store plots", type=Path, required=True
174+
)
175+
parser.add_argument("--style-file", help="colorbar definitions", type=Path)
176+
parser.add_argument("--plot-resolution", help="plot resolution", type=int)
177+
parser.add_argument(
178+
"input", help="grid_stat directories to process", type=Path, nargs="+"
179+
)
180+
args = parser.parse_args()
181+
182+
logging.basicConfig(level=logging.INFO)
183+
184+
for grid_stat_dir in args.input:
185+
cubes = load_grid_stat(grid_stat_dir.glob("**/*.nc"))
186+
187+
for cube in cubes:
188+
plot_grid_stat(
189+
cube,
190+
args.web_dir,
191+
style_file=args.style_file,
192+
plot_resolution=args.plot_resolution,
193+
)
194+
195+
196+
if __name__ == "__main__":
197+
main()

‎src/CSET/cset_workflow/app/grid_stat_plot/rose-app.conf‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,2 @@
1+
[command]
2+
default = grid_stat_plot.py --web-dir "${CYLC_WORKFLOW_SHARE_DIR}/web/plots/${CYLC_TASK_CYCLE_POINT}" "${CYLC_TASK_SHARE_CYCLE_DIR}/grid_stat"

‎src/CSET/cset_workflow/app/metplus_grid_stat/file/nci-gadi/gpm.conf‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,12 @@
1+
[config]
2+
3+
OBTYPE = GPM
4+
5+
# Daily files contain times 0015 to 2345
6+
# Make sure valid time 0000 uses the previous day's file by adding a shift of -1s to valid times
7+
OBS_GRID_STAT_INPUT_DIR = /g/data/dp9/da/verification/satellite/products/gpm/netcdf/imerg/NRTlate
8+
OBS_GRID_STAT_INPUT_TEMPLATE = {valid?fmt=%Y}/gpm_imerg_NRTlate_V07B_{valid?fmt=%Y%m%d?shift=-1}.nc
9+
10+
OBS_VAR1_NAME = precipitation_flux
11+
FCST_VAR1_NAME = stratiform_rainfall_flux
12+
BOTH_VAR1_LEVELS = "({valid?fmt=%Y%m%d_%H%M%S},*,*)"

‎src/CSET/cset_workflow/app/metplus_grid_stat/file/nci-gadi/gridstat.conf‎

Lines changed: 56 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,56 @@
1+
# Common METplus grid_stat configuration to make it work with CSET
2+
3+
[dir]
4+
MET_INSTALL_DIR = {ENV[CONDA_PREFIX]}
5+
6+
[config]
7+
PROCESS_LIST = GridStat
8+
9+
OUTPUT_BASE = {ENV[CYLC_TASK_SHARE_CYCLE_DIR]}
10+
LOG_DIR = {ENV[CYLC_TASK_LOG_DIR]}
11+
12+
MODEL = {ENV[MODEL_NAME]}
13+
14+
# [!] Set this for your site/obs type in a site file
15+
# OBTYPE =
16+
17+
###
18+
# Time Info
19+
# LOOP_BY options are INIT, VALID, RETRO, and REALTIME
20+
# If set to INIT or RETRO:
21+
# INIT_TIME_FMT, INIT_BEG, INIT_END, and INIT_INCREMENT must also be set
22+
# If set to VALID or REALTIME:
23+
# VALID_TIME_FMT, VALID_BEG, VALID_END, and VALID_INCREMENT must also be set
24+
# LEAD_SEQ is the list of forecast leads to process
25+
# https://metplus.readthedocs.io/en/latest/Users_Guide/systemconfiguration.html#timing-control
26+
###
27+
28+
LOOP_BY = INIT
29+
30+
INIT_TIME_FMT = %Y%m%dT%H%MZ
31+
INIT_BEG = {ENV[CYLC_TASK_CYCLE_POINT]}
32+
INIT_END = {ENV[CYLC_TASK_CYCLE_POINT]}
33+
INIT_INCREMENT = 1H
34+
35+
LEAD_SEQ = {ENV[LEAD_SEQ]}
36+
37+
###
38+
# File I/O
39+
# https://metplus.readthedocs.io/en/latest/Users_Guide/systemconfiguration.html#directory-and-filename-template-info
40+
###
41+
42+
FCST_GRID_STAT_INPUT_DIR = {ENV[CYLC_WORKFLOW_SHARE_DIR]}/cycle/{init?fmt=%Y%m%dT%H%MZ}/metplus_fcst
43+
FCST_GRID_STAT_INPUT_TEMPLATE = {ENV[MODEL_ID]}.nc
44+
45+
# [!] Set this for your site/obs type in a site file
46+
# OBS_GRID_STAT_INPUT_DIR =
47+
# OBS_GRID_STAT_INPUT_TEMPLATE =
48+
49+
GRID_STAT_OUTPUT_DIR = {OUTPUT_BASE}
50+
GRID_STAT_OUTPUT_TEMPLATE = grid_stat/{MODEL}/{OBTYPE}
51+
GRID_STAT_OUTPUT_PREFIX = {MODEL}_vs_{OBTYPE}_{CURRENT_OBS_NAME}
52+
53+
# Default regridding
54+
GRID_STAT_REGRID_TO_GRID = FCST
55+
GRID_STAT_REGRID_METHOD = BILIN
56+
GRID_STAT_REGRID_WIDTH = 2

‎src/CSET/cset_workflow/app/metplus_grid_stat/opt/rose-app-nci-gadi.conf‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
[command]
2+
default = for TYPE in $GRID_STAT_TYPES; do app_env_wrapper run_metplus.py "$SITE/common.conf" "$SITE/$TYPE.conf"; done
3+
4+
[env]
5+
!CONDA_VENV_LOCATION=/dev/null
6+
MET_INSTALL_DIR=/bom-ngm/conda
7+
8+
# TODO: put in top-level
9+
GRID_STAT_TYPES=gpm

src/CSET/cset_workflow/app/metplus_point_stat/file/nci/pointstat.conf renamed to src/CSET/cset_workflow/app/metplus_point_stat/file/nci-gadi/pointstat.conf

Lines changed: 13 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -1,24 +1,15 @@
1+
# Common METplus point_stat configuration to make it work with CSET
2+
13
[dir]
2-
FCST_BASE = {ENV[METPLUS_FCST_DIR]}
3-
OUTPUT_BASE = {ENV[CYLC_TASK_WORK_DIR]}
44
MET_INSTALL_DIR = {ENV[CONDA_PREFIX]}
55

66
[config]
7+
PROCESS_LIST = PointStat
78

9+
OUTPUT_BASE = {ENV[CYLC_TASK_SHARE_CYCLE_DIR]}
810
LOG_DIR = {ENV[CYLC_TASK_LOG_DIR]}
911

10-
# Documentation for this use case can be found at
11-
# https://metplus.readthedocs.io/en/latest/generated/met_tool_wrapper/ASCII2NC/ASCII2NC_python_embedding.html
12-
13-
# For additional information, please see the METplus Users Guide.
14-
# https://metplus.readthedocs.io/en/latest/Users_Guide
15-
16-
###
17-
# Processes to run
18-
# https://metplus.readthedocs.io/en/latest/Users_Guide/systemconfiguration.html#process-list
19-
###
20-
21-
PROCESS_LIST = PointStat
12+
MODEL = {ENV[MODEL_NAME]}
2213

2314
###
2415
# Time Info
@@ -32,25 +23,24 @@ PROCESS_LIST = PointStat
3223
###
3324

3425
LOOP_BY = INIT
35-
INIT_TIME_FMT = %Y%m%dT%H
36-
INIT_BEG = {ENV[TASK_START_TIME]}
37-
INIT_END = {ENV[TASK_START_TIME]}
26+
27+
INIT_TIME_FMT = %Y%m%dT%H%MZ
28+
INIT_BEG = {ENV[CYLC_TASK_CYCLE_POINT]}
29+
INIT_END = {ENV[CYLC_TASK_CYCLE_POINT]}
3830
INIT_INCREMENT = 1H
39-
LEAD_SEQ = begin_end_incr(0,{ENV[FORECAST_LENGTH]},1)
4031

41-
# Number of seconds to shift times in the fcst file (try half the time increment)
42-
FCST_SHIFT = 1800
32+
LEAD_SEQ = {ENV[LEAD_SEQ]}
4333

4434
###
4535
# File I/O
4636
# https://metplus.readthedocs.io/en/latest/Users_Guide/systemconfiguration.html#directory-and-filename-template-info
4737
###
4838

49-
FCST_POINT_STAT_INPUT_DIR = {FCST_BASE}
50-
FCST_POINT_STAT_INPUT_TEMPLATE = m1.nc
39+
FCST_POINT_STAT_INPUT_DIR = {ENV[CYLC_WORKFLOW_SHARE_DIR]}/cycle/{init?fmt=%Y%m%dT%H%MZ}/metplus_fcst
40+
FCST_POINT_STAT_INPUT_TEMPLATE = {ENV[MODEL_ID]}.nc
5141

5242
OBS_POINT_STAT_INPUT_DIR = {ENV[CYLC_WORKFLOW_SHARE_DIR]}/obs_nc
5343
OBS_POINT_STAT_INPUT_TEMPLATE = {valid?fmt=%Y%m%dT%H}.nc
5444

55-
POINT_STAT_OUTPUT_DIR = {ENV[CYLC_TASK_SHARE_CYCLE_DIR]}
45+
POINT_STAT_OUTPUT_DIR = {OUTPUT_BASE}
5646
POINT_STAT_OUTPUT_TEMPLATE = Point_Stat_{ENV[MODEL_NAME]}

src/CSET/cset_workflow/app/metplus_point_stat/file/nci/surface.conf renamed to src/CSET/cset_workflow/app/metplus_point_stat/file/nci-gadi/surface.conf

File renamed without changes.

‎src/CSET/cset_workflow/app/metplus_point_stat/opt/rose-app-nci-gadi.conf‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
[command]
2-
default=app_env_wrapper run_metplus.py nci/pointstat.conf nci/surface.conf
2+
default=app_env_wrapper run_metplus.py "$SITE/pointstat.conf" "$SITE/surface.conf"
33

4-
# Don't inherit rose environment
5-
[!env]
4+
[env]
5+
!CONDA_VENV_LOCATION=/dev/null
6+
MET_INSTALL_DIR=/bom-ngm/conda

0 commit comments

Comments
 (0)