Masked Statistics

Calculating Statistics for Regional Sub-Domains¶

Description:¶

This recipe shows how to calculate statistics for regional sub-domains using monthly-mean outputs from the National Oceanography Centre Near-Present-Day global eORCA1 configuration of NEMO forced using JRA55-do from 1976-2024.

For more details on this model configuration and the available outputs, users can explore the Near-Present-Day documentation here.

In [9]:

Copied!





# -- Import required packages -- #
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from nemo_cookbook import NEMODataTree
# -- Import required packages -- #
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from nemo_cookbook import NEMODataTree

Using Dask¶

Optional: Connect Client to Dask Local Cluster to run analysis in parallel.

Note that, although using Dask is not strictly necessary for this simple example using eORCA1, if we wanted to generalise this recipe to eORCA025 or eORCA12 outputs, using Dask would be essential to avoid unnecessary slow calculations using only a single process.

In [ ]:

Copied!





# -- Initialise Dask Local Cluster -- #
import os
import dask
from dask.distributed import Client, LocalCluster

# Update temporary directory for Dask workers:
dask.config.set({'temporary_directory': f"{os.getcwd()}/dask_tmp",
                 'local_directory': f"{os.getcwd()}/dask_tmp"
                 })

# Create Local Cluster:
cluster = LocalCluster(n_workers=4, threads_per_worker=3, memory_limit='5GB')
client = Client(cluster)
client
# -- Initialise Dask Local Cluster -- #
import os
import dask
from dask.distributed import Client, LocalCluster

# Update temporary directory for Dask workers:
dask.config.set({'temporary_directory': f"{os.getcwd()}/dask_tmp",
                 'local_directory': f"{os.getcwd()}/dask_tmp"
                 })

# Create Local Cluster:
cluster = LocalCluster(n_workers=4, threads_per_worker=3, memory_limit='5GB')
client = Client(cluster)
client

Accessing NEMO Model Data¶

Let's begin by loading the grid variables for our eORCA1 JRA-55 model from the JASMIN Object Store.

Alternatively, you can replace the domain_filepath below with a file path to your domain_cfg.nc file and read this with xarray's open_dataset() function.

In [3]:

Copied!

# Define directory path to ancillary files:
domain_filepath = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/domain"

# Open eORCA1 NEMO model domain_cfg:
ds_domain = xr.open_zarr(domain_filepath, consolidated=True, chunks={})

ds_domain
# Define directory path to ancillary files:
domain_filepath = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/domain"

# Open eORCA1 NEMO model domain_cfg:
ds_domain = xr.open_zarr(domain_filepath, consolidated=True, chunks={})

ds_domain

Out[3]:

<xarray.Dataset> Size: 669MB
Dimensions:        (y: 331, x: 360, t: 1, z: 75)
Dimensions without coordinates: y, x, t, z
Data variables: (12/58)
    atlmsk         (y, x) float32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
    bathy_metry    (t, y, x) float32 477kB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
    bottom_level   (t, y, x) int32 477kB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
    e1f            (t, y, x) float64 953kB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
    e1t            (t, y, x) float64 953kB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
    e1u            (t, y, x) float64 953kB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
    ...             ...
    nav_lev        (z) float32 300B dask.array<chunksize=(75,), meta=np.ndarray>
    nav_lon        (y, x) float32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
    pacmsk         (y, x) float32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
    socmsk         (y, x) float32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
    time_counter   (t) float64 8B dask.array<chunksize=(1,), meta=np.ndarray>
    top_level      (t, y, x) int32 477kB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
Attributes:
    NCO:      netCDF Operators version 5.2.2 (Homepage = http://nco.sf.net, C...
    history:  Fri Oct  4 10:58:13 2024: ncks -O --4 --no_abc --cnk_plc=xpl --...

Next, we will import the sea surface temperature and sea surface salinity stored at T-points in a single dataset.

Typically, NEMO model outputs defined on T-grid points are stored together in netCDF files. In this case, you can replace xr.merge() with a single call to xarray's open_dataset() function passing the file path to your _gridT.nc file(s).

In [4]:

Copied!

# Define directory path to model output files:
output_dir = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/T1m"

# Construct NEMO model grid dataset, including vertical grid cell thicknesses (m) and meridional velocities (m/s):
ds_gridT = xr.merge([xr.open_zarr(f"{output_dir}/{var}", consolidated=True, chunks={})[var] for var in ['tos_con', 'sos_abs']], compat="override")

ds_gridT
# Define directory path to model output files:
output_dir = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/T1m"

# Construct NEMO model grid dataset, including vertical grid cell thicknesses (m) and meridional velocities (m/s):
ds_gridT = xr.merge([xr.open_zarr(f"{output_dir}/{var}", consolidated=True, chunks={})[var] for var in ['tos_con', 'sos_abs']], compat="override")

ds_gridT

Out[4]:

<xarray.Dataset> Size: 552MB
Dimensions:        (y: 331, x: 360, time_counter: 577)
Coordinates:
    nav_lat        (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
    nav_lon        (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
    time_centered  (time_counter) datetime64[ns] 5kB dask.array<chunksize=(1,), meta=np.ndarray>
  * time_counter   (time_counter) datetime64[ns] 5kB 1976-01-16T12:00:00 ... ...
Dimensions without coordinates: y, x
Data variables:
    tos_con        (time_counter, y, x) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
    sos_abs        (time_counter, y, x) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
Attributes:
    cell_methods:        time: mean (interval: 3600 s)
    interval_operation:  3600 s
    interval_write:      1 month
    long_name:           sea_surface_conservative_temperature
    online_operation:    average
    standard_name:       sea_surface_temperature
    units:               degC

Creating a NEMODataTree¶

Next, let's create a NEMODataTree to store our domain and V-grid variables for the eORCA1 model.

In [5]:

Copied!





# Define dictionary of grid datasets defining eORCA1 parent model domain with no child/grand-child nests:
# Note: domain_cfg z-dimension is expected to be named 'nav_lev'.
datasets = {"parent": {"domain": ds_domain.rename({"z": "nav_lev"}), "gridT": ds_gridT}}

# Initialise a new NEMODataTree whose parent domain is zonally periodic & north-folding on F-points:
nemo = NEMODataTree.from_datasets(datasets=datasets, iperio=True, nftype="F")

nemo
# Define dictionary of grid datasets defining eORCA1 parent model domain with no child/grand-child nests:
# Note: domain_cfg z-dimension is expected to be named 'nav_lev'.
datasets = {"parent": {"domain": ds_domain.rename({"z": "nav_lev"}), "gridT": ds_gridT}}

# Initialise a new NEMODataTree whose parent domain is zonally periodic & north-folding on F-points:
nemo = NEMODataTree.from_datasets(datasets=datasets, iperio=True, nftype="F")

nemo

Out[5]:

<xarray.DatasetView> Size: 9kB
Dimensions:        (time_counter: 577)
Coordinates:
    time_centered  (time_counter) datetime64[ns] 5kB dask.array<chunksize=(1,), meta=np.ndarray>
  * time_counter   (time_counter) datetime64[ns] 5kB 1976-01-16T12:00:00 ... ...
Data variables:
    *empty*
Attributes:
    nftype:   F
    iperio:   True

Defining a Regional Sub-Domain using a Bounding Box¶

Now we have constructed our NEMODataTree, let's start by defining a regional sub-domain using a geographical bounding box.

By using the clip_grid() method, we permanently modify the size of the specfied grid stored in our NEMODataTree.

Alternatively, we can use clip_domain() to clip all of the grids associated with a given NEMO model domain to a given bounding box.

In [6]:

Copied!





# Define bounding box (lon_min, lon_max, lat_min, lat_max)
bbox = (-80, 10, 20, 70)

# Clip eORCA1 model T-grid to bounding box:
nemo.clip_grid(grid='/gridT', bbox=bbox)

# Plotting time-mean sea surface temperature for the regional sub-domain:
nemo['/gridT']['tos_con'].mean(dim='time_counter').plot()
# Define bounding box (lon_min, lon_max, lat_min, lat_max)
bbox = (-80, 10, 20, 70)

# Clip eORCA1 model T-grid to bounding box:
nemo.clip_grid(grid='/gridT', bbox=bbox)

# Plotting time-mean sea surface temperature for the regional sub-domain:
nemo['/gridT']['tos_con'].mean(dim='time_counter').plot()

Out[6]:

<matplotlib.collections.QuadMesh at 0x11df692b0>

No description has been provided for this image

Defining a Regional Sub-Domain using a Polygon¶

Next, let's define a more complex regional sub-domain by constructing a mask using a polygon. Since we have already clipped the T-grid of our NEMODataTree parent domain, we will define a polygon comprised of longitude-latitude coordinates within this region.

We will use the Overturning in the Subpolar North Atlantic Program (OSNAP) observational array coordinates made available via the JASMIN Object Store to construct a polygon enclosing the North Atlantic subpolar gyre.

In [7]:

Copied!





# Open OSNAP gridded observations dataset: 
ds_osnap = xr.open_zarr("https://noc-msm-o.s3-ext.jc.rl.ac.uk/ocean-obs/OSNAP/OSNAP_Gridded_TSV_201408_202006_2023")

# Define a closed polygon which includes both the OSNAP West & East arrays:
lon_poly = np.concatenate([ds_osnap['LONGITUDE'].values, np.array([ds_osnap['LONGITUDE'][-1], ds_osnap['LONGITUDE'][0]])])
lat_poly = np.concatenate([ds_osnap['LATITUDE'].values, np.array([ds_osnap['LATITUDE'][0], ds_osnap['LATITUDE'][0]])])
# Open OSNAP gridded observations dataset: 
ds_osnap = xr.open_zarr("https://noc-msm-o.s3-ext.jc.rl.ac.uk/ocean-obs/OSNAP/OSNAP_Gridded_TSV_201408_202006_2023")

# Define a closed polygon which includes both the OSNAP West & East arrays:
lon_poly = np.concatenate([ds_osnap['LONGITUDE'].values, np.array([ds_osnap['LONGITUDE'][-1], ds_osnap['LONGITUDE'][0]])])
lat_poly = np.concatenate([ds_osnap['LATITUDE'].values, np.array([ds_osnap['LATITUDE'][0], ds_osnap['LATITUDE'][0]])])

Now we have defined our polygon, we can use the mask_with_polygon() method to return the boolean mask classifying whether NEMO model grid points are inside (True) or outside (False) the polygon

In [18]:

Copied!





# Masking T-grid using polygon coordinates:
mask_spg = nemo.mask_with_polygon(grid='/gridT', lon_poly=lon_poly, lat_poly=lat_poly)

# Plotting SPG polygon sub-domain:
plt.figure()
plt.pcolormesh(nemo['/gridT']['glamt'], nemo['/gridT']['gphit'], nemo['/gridT']['tos_con'].mean(dim='time_counter'), cmap='RdBu_r')
plt.plot(lon_poly, lat_poly, color='0.1', lw=2)
plt.colorbar(label='Sea Surface Temperature (°C)')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.show()

# Plotting time-mean sea surface temperature for the SPG polygon sub-domain:
plt.figure()
plt.pcolormesh(nemo['/gridT']['glamt'], nemo['/gridT']['gphit'], nemo['/gridT']['tos_con'].mean(dim='time_counter').where(mask_spg), cmap='RdBu_r')
plt.plot(lon_poly, lat_poly, color='0.1', lw=2)
plt.colorbar(label='Sea Surface Temperature (°C)')
plt.xlabel('Longitude')
plt.xlim([-70, 5])
plt.ylabel('Latitude')
plt.ylim([40, 70])
plt.show()
# Masking T-grid using polygon coordinates:
mask_spg = nemo.mask_with_polygon(grid='/gridT', lon_poly=lon_poly, lat_poly=lat_poly)

# Plotting SPG polygon sub-domain:
plt.figure()
plt.pcolormesh(nemo['/gridT']['glamt'], nemo['/gridT']['gphit'], nemo['/gridT']['tos_con'].mean(dim='time_counter'), cmap='RdBu_r')
plt.plot(lon_poly, lat_poly, color='0.1', lw=2)
plt.colorbar(label='Sea Surface Temperature (°C)')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.show()

# Plotting time-mean sea surface temperature for the SPG polygon sub-domain:
plt.figure()
plt.pcolormesh(nemo['/gridT']['glamt'], nemo['/gridT']['gphit'], nemo['/gridT']['tos_con'].mean(dim='time_counter').where(mask_spg), cmap='RdBu_r')
plt.plot(lon_poly, lat_poly, color='0.1', lw=2)
plt.colorbar(label='Sea Surface Temperature (°C)')
plt.xlabel('Longitude')
plt.xlim([-70, 5])
plt.ylabel('Latitude')
plt.ylim([40, 70])
plt.show()

Calculating statistics for a Regional Sub-Domain¶

Finally, let's use our North Atlantic subpolar gyre polygon to calculate statistics for this regional sub-domain of the eORCA1 model.

Given a closed polygon, we can use the masked_statistic() method to calculate statistics of a specified variable in the masked sub-domain

In [17]:

Copied!





# Calculating the area weighted-mean sea surface temperature in the SPG region:
sst_wmean = nemo.masked_statistic(grid="/gridT",
                                  var="tos_con",
                                  lon_poly=lon_poly,
                                  lat_poly=lat_poly,
                                  statistic="weighted_mean",
                                  dims=["i", "j"]
                                  )

sst_wmean.plot(lw=1, color='0.1', alpha=0.3)
sst_wmean.rolling(time_counter=12, center=True).mean().plot(lw=3, color='0.1')
plt.title('Area Weighted Mean SST for SPG Region', fontsize=12, fontweight='bold')
plt.xlabel('Time', fontsize=12)
plt.ylabel('Sea Surface Temperature (°C)', fontsize=12)
# Calculating the area weighted-mean sea surface temperature in the SPG region:
sst_wmean = nemo.masked_statistic(grid="/gridT",
                                  var="tos_con",
                                  lon_poly=lon_poly,
                                  lat_poly=lat_poly,
                                  statistic="weighted_mean",
                                  dims=["i", "j"]
                                  )

sst_wmean.plot(lw=1, color='0.1', alpha=0.3)
sst_wmean.rolling(time_counter=12, center=True).mean().plot(lw=3, color='0.1')
plt.title('Area Weighted Mean SST for SPG Region', fontsize=12, fontweight='bold')
plt.xlabel('Time', fontsize=12)
plt.ylabel('Sea Surface Temperature (°C)', fontsize=12)

Out[17]:

Text(0, 0.5, 'Sea Surface Temperature (°C)')

In [ ]: