Comparing MOM6 output with in-situ ocean observations (small domain, serial)

In this tutorial we will get familiar with the basics of using model2obs to interpolate MOM6 output onto the space of in-situ observations stored in CrocoLake. By the end of it, you will learn to:

explore data in CrocoLake
interpolate MOM6 output into CrocoLake observations space
display the results on an interactive map
display the results on an interactive 2D plot

Reminder: we’re looking at the machinery, not really at the science (that is up to you!)

This version uses output from a local CESM-MOM6 run around the Hawaiian islands (case hawaii_pmo, 3 daily snapshots: 2023-06-16, 2023-06-17 and 2023-06-18, on a grid spanning 20.03-24.98 N and 155.03-159.98 W). The CrocoLake observation files were pre-generated and are already in place.

Model data is read from /glade/derecho/scratch/emilanese/m2o_data/hawaii/mom6, and everything model2obs produces is stored under /glade/derecho/scratch/emilanese/m2o_data/hawaii/m2o.

CrocoLake¶

CrocoLake is a flexible parquet dataset that is fast and easy to access and filter. A copy of the CrocoLake version with physical variables ships with model2obs and is pointed at by $CROCOLAKE_PATH, so let’s set it up:

import os

crocolake_path = os.path.expandvars('$CROCOLAKE_PATH')
print(f'Reading CrocoLake from: {crocolake_path}')

We now have a look at the CrocoLake data that we are going to work with: we load the observations in the region and time covered by the model output, and plot their positions on a map.

The region boundaries used here (LAT0, LAT1, LON0, LON1) are the same ones that were used to build the obs_seq files: deliberately a bit wider than the model domain, so that some observations fall outside the grid.

%%time
import cartopy.crs as ccrs
import dask.dataframe as dd
import datetime
import matplotlib.pyplot as plt

# region of interest: slightly wider than the Hawaii model domain (20.03-24.98 N,
# 155.03-159.98 W), so that some observations fall outside the model grid
LAT0, LAT1 = 18, 27
LON0, LON1 = -162, -153     # CrocoLake stores longitude in the -180:180 convention

cl = dd.read_parquet(
    crocolake_path,
    columns=["DB_NAME", "LATITUDE", "LONGITUDE", "PRES", "TEMP", "PSAL"],
    filters=[("LATITUDE", ">", LAT0), ("LATITUDE", "<", LAT1),
             ("LONGITUDE", ">", LON0), ("LONGITUDE", "<", LON1),
             ("JULD", ">", datetime.datetime(2023, 6, 16)),
             ("JULD", "<", datetime.datetime(2023, 6, 19))],
).compute()
print(f"{len(cl)} observations from {cl['DB_NAME'].unique().tolist()}, "
      f"{len(cl.groupby(['LATITUDE', 'LONGITUDE']))} profiles")

ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([LON0, LON1, LAT0, LAT1])
ax.coastlines()
ax.gridlines(draw_labels=True)
ax.plot(cl["LONGITUDE"], cl["LATITUDE"], 'o', ms=5, label="CrocoLake profiles")

Running the workflow¶

The obs_seq files that model2obs will consume have already been generated from CrocoLake (one file per model snapshot, see data assimilation tutorial for how to build them) and are in the folder declared as obs_seq_in_folder in model2obs.config_tutorial_hawaii.yaml. Their names sort in the same order as the model output files, which is what model2obs relies on when pairing them up. Note: per model2obs design, obs_seq_in_folder must contain only obs_seq files.

Running the workflow to interpolate the model is quite simple: create a WorkflowModelObs instance using a configuration file, and then call its run() method. For this tutorial we will use the configuration file model2obs.config_tutorial_hawaii.yaml. A template file to use as reference is also provided at config_templates/config_template_MOM6.yaml.

While running, model2obs generates temporary input files that tell DART’s perfect_model_obs executable where to find MOM6 and CrocoLake data to perform the interpolation.

from model2obs.workflows import WorkflowModelObs

# Create and run workflow to interpolate MOM6 model onto CrocoLake obs space
workflow_crocolake = WorkflowModelObs.from_config_file('model2obs.config_tutorial_hawaii.yaml')
workflow_crocolake.run(clear_output=True) #use flag clear_output=True if you want to re-run it and automatically clean all previous output

If you look through the log, you’ll see a passage that says:

           Convex hull boundaries: 
            Longitude: 200.03° to 204.97° (0-360 convention)
            Latitude: 20.02° to 24.98°
           Number of observations before filtering to convex hull: 4980
           Number of observations after filtering to convex hull: 2968
           Percentage of original observations retained: 59.6%

This is model2obs automatically removing observations that fall outside the MOM6 domain.

Output and bias estimate¶

model2obs generates a parquet dataset that contains the values of the CrocoLake observations, the MOM6 model data interpolated onto the observations space, and some basic statistics.

The dataframe by default is stored to model_obs_df (as a dask dataframe), and model2obs offers tools to access all data, only the successful interpolations, or only the failed interpolations. For now, let’s load only the interpolations that succeeded (the output message from the previous cell reports how many of the total they are):

good_model_obs_df = workflow_crocolake.get_good_model_obs_df(compute=True) # compute=True triggers the compute of the dask dataframe, returning a pandas dataframe with data loaded in memory
good_model_obs_df.head()                                                   # displays first 5 rows in the dataframe

With this table, we can already look at the overall bias of our model with respect to the observations. We do so by computing the mean of the difference variable, which is defined as obs - interpolated_model:

bias = good_model_obs_df.groupby('type')['difference'].agg(
    n='count',
    bias='mean',
    std='std',
    rmse=lambda d: (d ** 2).mean() ** 0.5,
)
bias

For this three-day Hawaii run the temperature bias is about -0.13 degC with an RMSE of 0.42 degC, i.e. the model is slightly warmer than the Argo observations, and the systematic part of the misfit is small compared to its scatter. The salinity numbers are tiny because CrocoLake stores salinity in kg/kg in this dataset, so read them as a relative rather than an absolute offset.

Displaying the interactive map¶

Loading the interactive map to explore the successful interpolations is as simple as importing the widget and passing the dataframe to it. We also pass a MapConfig to frame the map on the Hawaii region:

from model2obs.viz import InteractiveWidgetMap, MapConfig

# Frame the map on the Hawaii region
map_config = MapConfig(
    map_extent=(LON0, LON1, LAT0, LAT1)  #(lon_min, lon_max, lat_min, lat_max)
)

# Create an interactive map widget to visualize model-observation comparisons
# The widget provides controls for selecting variables, observation types, and time ranges
widget = InteractiveWidgetMap(good_model_obs_df, config=map_config)
widget.setup()

Comparing a single profile¶

The map above mixes observations coming from different platforms, depths and times. To judge how the model does against a single vertical profile, we can isolate the observations that belong to the same profile, i.e. those sharing the same latitude, longitude and time, and then use the interactive profile widget.

Let’s first see how many distinct profiles we have, and pick the one with the most observations:

# a CrocoLake profile is identified by a unique (latitude, longitude, time) triplet
profiles = (
    good_model_obs_df
    .groupby(['latitude', 'longitude', 'time'])
    .size()
    .sort_values(ascending=False)
)
print(f'Number of distinct profiles: {len(profiles)}')
print(profiles)

# select the profile with the largest number of observations
lat0, lon0, time0 = profiles.index[0]
profile_df = good_model_obs_df[
    (good_model_obs_df['latitude'] == lat0)
    & (good_model_obs_df['longitude'] == lon0)
    & (good_model_obs_df['time'] == time0)
]
print(f'\nSelected profile at ({lat0} N, {lon0} E) on {time0}: '
      f'{len(profile_df)} observations of types {sorted(profile_df["type"].unique())}')
profile_df.head()

Interactive profile¶

InteractiveWidgetProfile plots any pair of columns of the dataframe against each other. Plotting obs (or interpolated_model) versus vertical gives us the vertical profile; the y-axis is inverted by default, so that depth grows downwards.

Use the dropdowns to swap the x-axis between obs, interpolated_model and difference, and the Types selector to switch between temperature and salinity (ctrl+click to select more than one type at a time):

from model2obs.viz import InteractiveWidgetProfile, ProfileConfig

profile_config = ProfileConfig(
    figure_size=(7, 7),
    marker_size=5,
    marker_alpha=0.6,
    invert_yaxis=True,   # depth grows downwards
    grid=True
)

widget_profile = InteractiveWidgetProfile(profile_df, x='obs', y='vertical', config=profile_config)
widget_profile.setup()

1:1 comparison¶

Putting the observations on the x-axis and the interpolated model on the y-axis instead gives the classic 1:1 (scatter) comparison: the closer the points fall to the diagonal, the better the model reproduces the observations. Here we also turn off the y-axis inversion, since neither axis is a depth anymore:

profile_config_1to1 = ProfileConfig(
    figure_size=(7, 7),
    marker_size=5,
    marker_alpha=0.6,
    invert_yaxis=False,  # not a depth axis anymore
    grid=True
)

widget_1to1 = InteractiveWidgetProfile(
    profile_df, x='obs', y='interpolated_model', config=profile_config_1to1
)
widget_1to1.setup()

Interpolation errors¶

The dataframe built by model2obs during WorkflowModelObs.run() contains the column interpolated_model_QC, which stores information about the quality of the interpolation. The value is set by DART’s perfect_model_obs program, and you can find more information about it here and here. In general, for this workflow we want QC≤2, and indeed the method get_good_model_obs_df() that we used earlier uses this criterion.

In this Hawaii run, however all interpolations were successful. The cells below are therefore left commented out: they are suggestions to follow whenever your own run reports a non-zero number of failed interpolations. Uncomment and run them in that case.

# failed_model_obs_df = workflow_crocolake.get_failed_model_obs_df(compute=True) # compute=True triggers the compute of the dask dataframe, returning a pandas dataframe with data loaded in memory
# failed_model_obs_df.head()                                                     # displays first 5 rows in the dataframe

# widget_failed = InteractiveWidgetMap(failed_model_obs_df, config=map_config)
# widget_failed.setup()

Had there been failures, the two cells above would have shown:

when an interpolation fails, a value of -888888 is assigned to the interpolated model, so never feed a raw dataframe to your statistics: always filter on the QC first (or use get_good_model_obs_df());
from head(), values of interpolated_model_QC greater than 1000;

The unique values of the QC flag are the quickest way to classify the failures:

# failed_model_obs_df['interpolated_model_QC'].unique()

Values greater than 1000 are given as 1000 + failure_code, with failure_code from here. The QC flag 1018 indicates that one or more grid points required for the interpolation are not in the basin, so the interpolation cannot be performed: this is the flag you get for observations outside the domain or close to the sea/land border, and it is the one you would cross-check against the map above.

Finally, if you want to load all the data at once (successful and failed interpolations together), you can use the following command:

# Load the parquet dataset generated by the workflow above
model_obs_df = workflow_crocolake.get_all_model_obs_df(compute=True) # compute=True triggers the compute of the dask dataframe, returning a pandas dataframe with data loaded in memory
model_obs_df.head() # displays first 5 rows in the dataframe