Regional Masked Statistics
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 [1]:
Copied!
# -- Import required packages -- #
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from nemo_cookbook import NEMODataTree
xr.set_options(display_style="text")
# -- Import required packages -- #
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from nemo_cookbook import NEMODataTree
xr.set_options(display_style="text")
Out[1]:
<xarray.core.options.set_options at 0x3156afb60>
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 url to domain_cfg ancillary file:
domain_url = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/domain_cfg"
# Open eORCA1 NEMO model domain_cfg:
ds_domain = xr.open_zarr(domain_url, consolidated=True, chunks={})
ds_domain
# Define url to domain_cfg ancillary file:
domain_url = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/domain_cfg"
# Open eORCA1 NEMO model domain_cfg:
ds_domain = xr.open_zarr(domain_url, consolidated=True, chunks={})
ds_domain
Out[3]:
<xarray.Dataset> Size: 709MB
Dimensions: (y: 331, x: 360, nav_lev: 75)
Coordinates:
* nav_lev (nav_lev) int64 600B 0 1 2 3 4 5 6 7 ... 68 69 70 71 72 73 74
* x (x) int64 3kB 0 1 2 3 4 5 6 7 ... 353 354 355 356 357 358 359
* y (y) int64 3kB 0 1 2 3 4 5 6 7 ... 324 325 326 327 328 329 330
Data variables: (12/49)
atlmsk (y, x) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>
bathy_metry (y, x) float32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
bottom_level (y, x) int32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
e1f (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
e1t (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
e1u (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
... ...
top_level (y, x) int32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
umask (nav_lev, y, x) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>
umaskutil (y, x) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>
vmask (nav_lev, y, x) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>
vmaskutil (y, x) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>
wmask (nav_lev, y, x) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>
Attributes:
CfgName: UNKNOWN
CfgIndex: -999
Iperio: 1
Jperio: 0
NFold: 1
NFtype: F
VertCoord: zps
IsfCav: 0
file_name: mesh_mask.nc
TimeStamp: 01/03/2025 22:19:49 -0000
Next, we will import the sea surface temperature and sea surface salinity stored at T-points in a single dataset.
In [4]:
Copied!
# Define url to T-grid NEMO model output files:
gridT_url = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/T1m"
# Open NEMO model grid dataset, including vertical grid cell thicknesses (m) and meridional velocities (m/s):
ds_gridT = xr.open_zarr(gridT_url, consolidated=True, chunks={})
ds_gridT
# Define url to T-grid NEMO model output files:
gridT_url = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/T1m"
# Open NEMO model grid dataset, including vertical grid cell thicknesses (m) and meridional velocities (m/s):
ds_gridT = xr.open_zarr(gridT_url, consolidated=True, chunks={})
ds_gridT
Out[4]:
<xarray.Dataset> Size: 467GB
Dimensions: (time_counter: 577, y: 331, x: 360, deptht: 75,
axis_nbounds: 2)
Coordinates:
* deptht (deptht) float32 300B 0.5058 1.556 ... 5.902e+03
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...
Dimensions without coordinates: y, x, axis_nbounds
Data variables: (12/74)
berg_latent_heat_flux (time_counter, y, x) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
deptht_bounds (deptht, axis_nbounds) float32 600B dask.array<chunksize=(25, 2), meta=np.ndarray>
e3t (time_counter, deptht, y, x) float32 21GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
empmr (time_counter, y, x) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
evs (time_counter, y, x) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
ficeberg (time_counter, y, x) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
... ...
ttrd_qns_li (time_counter, y, x) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
vohfcisf (time_counter, deptht, y, x) float32 21GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
vohflisf (time_counter, deptht, y, x) float32 21GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
vowflisf (time_counter, deptht, y, x) float32 21GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
zos (time_counter, y, x) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
zossq (time_counter, y, x) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
Attributes:
name: OUTPUT/eORCA1_1m_grid_T
description: ocean T grid variables
title: ocean T grid variables
Conventions: CF-1.6
timeStamp: 2024-Dec-09 10:18:50 GMT
uuid: d4f9020a-a621-4a95-81da-0b477f67ed38
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:
datasets = {"parent": {"domain": ds_domain, "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", read_mask=True)
nemo
# Define dictionary of grid datasets defining eORCA1 parent model domain with no child/grand-child nests:
datasets = {"parent": {"domain": ds_domain, "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", read_mask=True)
nemo
Out[5]:
<xarray.DataTree 'NEMO model'>
Group: /
│ Dimensions: (time_counter: 577, axis_nbounds: 2)
│ 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:...
│ Dimensions without coordinates: axis_nbounds
│ Data variables:
│ time_centered_bounds (time_counter, axis_nbounds) datetime64[ns] 9kB dask.array<chunksize=(1, 2), meta=np.ndarray>
│ time_counter_bounds (time_counter, axis_nbounds) datetime64[ns] 9kB dask.array<chunksize=(1, 2), meta=np.ndarray>
│ Attributes:
│ nftype: F
│ iperio: True
├── Group: /gridT
│ Dimensions: (time_counter: 577, j: 331, i: 360, k: 75,
│ axis_nbounds: 2)
│ Coordinates:
│ * deptht (k) float32 300B 0.5058 1.556 ... 5.698e+03 5.902e+03
│ time_centered (time_counter) datetime64[ns] 5kB dask.array<chunksize=(1,), meta=np.ndarray>
│ gphit (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ glamt (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ * k (k) int64 600B 1 2 3 4 5 6 7 ... 69 70 71 72 73 74 75
│ * j (j) int64 3kB 1 2 3 4 5 6 ... 326 327 328 329 330 331
│ * i (i) int64 3kB 1 2 3 4 5 6 ... 355 356 357 358 359 360
│ Dimensions without coordinates: axis_nbounds
│ Data variables: (12/80)
│ berg_latent_heat_flux (time_counter, j, i) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
│ deptht_bounds (k, axis_nbounds) float32 600B dask.array<chunksize=(25, 2), meta=np.ndarray>
│ e3t (time_counter, k, j, i) float32 21GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
│ empmr (time_counter, j, i) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
│ evs (time_counter, j, i) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
│ ficeberg (time_counter, j, i) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
│ ... ...
│ e1t (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ e2t (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ top_level (j, i) int32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ bottom_level (j, i) int32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ tmask (k, j, i) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>
│ tmaskutil (j, i) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ Attributes:
│ nftype: F
│ iperio: True
├── Group: /gridU
│ Dimensions: (j: 331, i: 360, k: 75)
│ Coordinates:
│ gphiu (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ glamu (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ * k (k) int64 600B 1 2 3 4 5 6 7 8 9 ... 68 69 70 71 72 73 74 75
│ * j (j) int64 3kB 1 2 3 4 5 6 7 8 ... 325 326 327 328 329 330 331
│ * i (i) float64 3kB 1.5 2.5 3.5 4.5 ... 357.5 358.5 359.5 360.5
│ Data variables:
│ e1u (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ e2u (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ umask (k, j, i) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>
│ umaskutil (j, i) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ Attributes:
│ nftype: F
│ iperio: True
├── Group: /gridV
│ Dimensions: (j: 331, i: 360, k: 75)
│ Coordinates:
│ gphiv (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ glamv (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ * k (k) int64 600B 1 2 3 4 5 6 7 8 9 ... 68 69 70 71 72 73 74 75
│ * j (j) float64 3kB 1.5 2.5 3.5 4.5 ... 328.5 329.5 330.5 331.5
│ * i (i) int64 3kB 1 2 3 4 5 6 7 8 ... 354 355 356 357 358 359 360
│ Data variables:
│ e1v (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ e2v (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ vmask (k, j, i) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>
│ vmaskutil (j, i) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ Attributes:
│ nftype: F
│ iperio: True
├── Group: /gridW
│ Dimensions: (j: 331, i: 360, k: 75)
│ Coordinates:
│ gphit (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ glamt (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ * k (k) float64 600B 0.5 1.5 2.5 3.5 4.5 ... 71.5 72.5 73.5 74.5
│ * j (j) int64 3kB 1 2 3 4 5 6 7 8 ... 325 326 327 328 329 330 331
│ * i (i) int64 3kB 1 2 3 4 5 6 7 8 ... 354 355 356 357 358 359 360
│ Data variables:
│ e1t (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ e2t (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ wmask (k, j, i) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>
│ wmaskutil (j, i) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ Attributes:
│ nftype: F
│ iperio: True
└── Group: /gridF
Dimensions: (j: 331, i: 360, k: 75)
Coordinates:
gphif (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
glamf (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
* k (k) int64 600B 1 2 3 4 5 6 7 8 9 ... 68 69 70 71 72 73 74 75
* j (j) float64 3kB 1.5 2.5 3.5 4.5 ... 328.5 329.5 330.5 331.5
* i (i) float64 3kB 1.5 2.5 3.5 4.5 ... 357.5 358.5 359.5 360.5
Data variables:
e1f (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
e2f (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
fmask (k, j, i) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>
fmaskutil (j, i) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>
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 can 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_clipped = nemo.clip_grid(grid='gridT', bbox=bbox)
nemo_clipped
# 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_clipped = nemo.clip_grid(grid='gridT', bbox=bbox)
nemo_clipped
Out[6]:
<xarray.DataTree 'NEMO model'>
Group: /
│ Dimensions: (time_counter: 577, axis_nbounds: 2)
│ 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:...
│ Dimensions without coordinates: axis_nbounds
│ Data variables:
│ time_centered_bounds (time_counter, axis_nbounds) datetime64[ns] 9kB dask.array<chunksize=(1, 2), meta=np.ndarray>
│ time_counter_bounds (time_counter, axis_nbounds) datetime64[ns] 9kB dask.array<chunksize=(1, 2), meta=np.ndarray>
│ Attributes:
│ nftype: F
│ iperio: True
├── Group: /gridT
│ Dimensions: (time_counter: 577, j: 83, i: 90, k: 75,
│ axis_nbounds: 2)
│ Coordinates:
│ * deptht (k) float32 300B 0.5058 1.556 ... 5.698e+03 5.902e+03
│ time_centered (time_counter) datetime64[ns] 5kB dask.array<chunksize=(1,), meta=np.ndarray>
│ gphit (j, i) float64 60kB 20.54 20.54 20.55 ... 71.07 70.74
│ glamt (j, i) float64 60kB -79.5 -78.5 -77.5 ... 31.15 32.2
│ * k (k) int64 600B 1 2 3 4 5 6 7 ... 69 70 71 72 73 74 75
│ * j (j) int64 664B 223 224 225 226 ... 302 303 304 305
│ * i (i) int64 720B 208 209 210 211 ... 294 295 296 297
│ Dimensions without coordinates: axis_nbounds
│ Data variables: (12/80)
│ berg_latent_heat_flux (time_counter, j, i) float32 17MB dask.array<chunksize=(1, 83, 90), meta=np.ndarray>
│ deptht_bounds (k, axis_nbounds, j, i) float32 4MB dask.array<chunksize=(25, 2, 83, 90), meta=np.ndarray>
│ e3t (time_counter, k, j, i) float32 1GB dask.array<chunksize=(1, 25, 83, 90), meta=np.ndarray>
│ empmr (time_counter, j, i) float32 17MB dask.array<chunksize=(1, 83, 90), meta=np.ndarray>
│ evs (time_counter, j, i) float32 17MB dask.array<chunksize=(1, 83, 90), meta=np.ndarray>
│ ficeberg (time_counter, j, i) float32 17MB dask.array<chunksize=(1, 83, 90), meta=np.ndarray>
│ ... ...
│ e1t (j, i) float64 60kB dask.array<chunksize=(83, 90), meta=np.ndarray>
│ e2t (j, i) float64 60kB dask.array<chunksize=(83, 90), meta=np.ndarray>
│ top_level (j, i) int32 30kB dask.array<chunksize=(83, 90), meta=np.ndarray>
│ bottom_level (j, i) int32 30kB dask.array<chunksize=(83, 90), meta=np.ndarray>
│ tmask (k, j, i) float32 2MB dask.array<chunksize=(75, 83, 90), meta=np.ndarray>
│ tmaskutil (j, i) float32 30kB dask.array<chunksize=(83, 90), meta=np.ndarray>
│ Attributes:
│ nftype: F
│ iperio: False
├── Group: /gridU
│ Dimensions: (j: 331, i: 360, k: 75)
│ Coordinates:
│ gphiu (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ glamu (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ * k (k) int64 600B 1 2 3 4 5 6 7 8 9 ... 68 69 70 71 72 73 74 75
│ * j (j) int64 3kB 1 2 3 4 5 6 7 8 ... 325 326 327 328 329 330 331
│ * i (i) float64 3kB 1.5 2.5 3.5 4.5 ... 357.5 358.5 359.5 360.5
│ Data variables:
│ e1u (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ e2u (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ umask (k, j, i) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>
│ umaskutil (j, i) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ Attributes:
│ nftype: F
│ iperio: True
├── Group: /gridV
│ Dimensions: (j: 331, i: 360, k: 75)
│ Coordinates:
│ gphiv (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ glamv (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ * k (k) int64 600B 1 2 3 4 5 6 7 8 9 ... 68 69 70 71 72 73 74 75
│ * j (j) float64 3kB 1.5 2.5 3.5 4.5 ... 328.5 329.5 330.5 331.5
│ * i (i) int64 3kB 1 2 3 4 5 6 7 8 ... 354 355 356 357 358 359 360
│ Data variables:
│ e1v (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ e2v (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ vmask (k, j, i) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>
│ vmaskutil (j, i) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ Attributes:
│ nftype: F
│ iperio: True
├── Group: /gridW
│ Dimensions: (j: 331, i: 360, k: 75)
│ Coordinates:
│ gphit (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ glamt (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ * k (k) float64 600B 0.5 1.5 2.5 3.5 4.5 ... 71.5 72.5 73.5 74.5
│ * j (j) int64 3kB 1 2 3 4 5 6 7 8 ... 325 326 327 328 329 330 331
│ * i (i) int64 3kB 1 2 3 4 5 6 7 8 ... 354 355 356 357 358 359 360
│ Data variables:
│ e1t (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ e2t (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ wmask (k, j, i) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>
│ wmaskutil (j, i) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ Attributes:
│ nftype: F
│ iperio: True
└── Group: /gridF
Dimensions: (j: 331, i: 360, k: 75)
Coordinates:
gphif (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
glamf (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
* k (k) int64 600B 1 2 3 4 5 6 7 8 9 ... 68 69 70 71 72 73 74 75
* j (j) float64 3kB 1.5 2.5 3.5 4.5 ... 328.5 329.5 330.5 331.5
* i (i) float64 3kB 1.5 2.5 3.5 4.5 ... 357.5 358.5 359.5 360.5
Data variables:
e1f (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
e2f (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
fmask (k, j, i) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>
fmaskutil (j, i) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>
Attributes:
nftype: F
iperio: True
In [7]:
Copied!
# -- Plotting time-mean sea surface temperature for the regional subdomain -- #
nemo_clipped['gridT/tos_con'].mean(dim='time_counter').plot()
# -- Plotting time-mean sea surface temperature for the regional subdomain -- #
nemo_clipped['gridT/tos_con'].mean(dim='time_counter').plot()
Out[7]:
<matplotlib.collections.QuadMesh at 0x341396120>
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 [8]:
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 [9]:
Copied!
# Masking T-grid using polygon coordinates:
mask_spg = nemo_clipped.mask_with_polygon(grid='gridT', lon_poly=lon_poly, lat_poly=lat_poly)
# -- Plotting SPG polygon sub-domain -- #
plt.figure()
plt.pcolormesh(nemo_clipped['gridT/glamt'], nemo_clipped['gridT/gphit'], nemo_clipped['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)')
# Axes labels:
plt.xlabel('Longitude', fontdict={'size':12, 'weight':'bold'})
plt.ylabel('Latitude', fontdict={'size':12, 'weight':'bold'})
plt.show()
# Masking T-grid using polygon coordinates:
mask_spg = nemo_clipped.mask_with_polygon(grid='gridT', lon_poly=lon_poly, lat_poly=lat_poly)
# -- Plotting SPG polygon sub-domain -- #
plt.figure()
plt.pcolormesh(nemo_clipped['gridT/glamt'], nemo_clipped['gridT/gphit'], nemo_clipped['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)')
# Axes labels:
plt.xlabel('Longitude', fontdict={'size':12, 'weight':'bold'})
plt.ylabel('Latitude', fontdict={'size':12, 'weight':'bold'})
plt.show()
/var/folders/z2/j_dr250s42x34hk63_rp4bm80000gq/T/ipykernel_65138/336358065.py:6: UserWarning: The input coordinates to pcolormesh are interpreted as cell centers, but are not monotonically increasing or decreasing. This may lead to incorrectly calculated cell edges, in which case, please supply explicit cell edges to pcolormesh. plt.pcolormesh(nemo_clipped['gridT/glamt'], nemo_clipped['gridT/gphit'], nemo_clipped['gridT/tos_con'].mean(dim='time_counter'), cmap='RdBu_r')
In [10]:
Copied!
# -- Plotting time-mean sea surface temperature for the SPG polygon -- #
plt.figure()
plt.pcolormesh(nemo_clipped['gridT/glamt'],
nemo_clipped['gridT/gphit'],
nemo_clipped['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)')
# Axes labels:
plt.title('eORCA1-JRA55v1:\nTime-Mean SST for SPG Region', fontdict={'size':12, 'weight':'bold'})
plt.xlabel('Longitude', fontdict={'size':12, 'weight':'bold'})
plt.xlim([-70, 5])
plt.ylabel('Latitude', fontdict={'size':12, 'weight':'bold'})
plt.ylim([40, 70])
plt.show()
# -- Plotting time-mean sea surface temperature for the SPG polygon -- #
plt.figure()
plt.pcolormesh(nemo_clipped['gridT/glamt'],
nemo_clipped['gridT/gphit'],
nemo_clipped['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)')
# Axes labels:
plt.title('eORCA1-JRA55v1:\nTime-Mean SST for SPG Region', fontdict={'size':12, 'weight':'bold'})
plt.xlabel('Longitude', fontdict={'size':12, 'weight':'bold'})
plt.xlim([-70, 5])
plt.ylabel('Latitude', fontdict={'size':12, 'weight':'bold'})
plt.ylim([40, 70])
plt.show()
/var/folders/z2/j_dr250s42x34hk63_rp4bm80000gq/T/ipykernel_65138/1564322244.py:3: UserWarning: The input coordinates to pcolormesh are interpreted as cell centers, but are not monotonically increasing or decreasing. This may lead to incorrectly calculated cell edges, in which case, please supply explicit cell edges to pcolormesh. plt.pcolormesh(nemo_clipped['gridT/glamt'],
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 [11]:
Copied!
# Calculating the area weighted-mean sea surface temperature in the SPG region:
sst_wmean = nemo_clipped["gridT/tos_con"].masked_statistic(lon_poly=lon_poly,
lat_poly=lat_poly,
statistic="weighted_mean",
dims=["i", "j"]
)
sst_wmean = sst_wmean.compute()
# Calculating the area weighted-mean sea surface temperature in the SPG region:
sst_wmean = nemo_clipped["gridT/tos_con"].masked_statistic(lon_poly=lon_poly,
lat_poly=lat_poly,
statistic="weighted_mean",
dims=["i", "j"]
)
sst_wmean = sst_wmean.compute()
In [13]:
Copied!
plt.figure(figsize=(8, 4))
# Plotting the area weighted-mean sea surface temperature time series for the SPG region:
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')
# Axes labels:
plt.title('eORCA1 JRA55v1:\nWeighted Mean SST for SPG Region', fontdict={'size':12, 'weight':'bold'})
plt.xlabel('Time', fontdict={'size':12, 'weight':'bold'})
plt.ylabel('Sea Surface Temperature (°C)', fontdict={'size':12, 'weight':'bold'})
plt.figure(figsize=(8, 4))
# Plotting the area weighted-mean sea surface temperature time series for the SPG region:
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')
# Axes labels:
plt.title('eORCA1 JRA55v1:\nWeighted Mean SST for SPG Region', fontdict={'size':12, 'weight':'bold'})
plt.xlabel('Time', fontdict={'size':12, 'weight':'bold'})
plt.ylabel('Sea Surface Temperature (°C)', fontdict={'size':12, 'weight':'bold'})
Out[13]:
Text(0, 0.5, 'Sea Surface Temperature (°C)')
