Skip to content

Commit 3901db1

Browse files
VeckoTheGeckoVeckoTheGecko
committed
Add nemo example data and comparison tooling
1 parent 0866440 commit 3901db1

2 files changed

Lines changed: 245 additions & 1 deletion

File tree

‎parcels/_datasets/structured/circulation_model.py‎

Lines changed: 153 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,5 @@
1+
import numpy as np
2+
import pandas as pd
13
import xarray as xr
24

35
from . import T, X, Y, Z
@@ -12,7 +14,157 @@ def _nemo_data() -> xr.Dataset:
1214
1315
https://www.mercator-ocean.eu/en/solutions-expertise/accessing-digital-data/product-details/?offer=4217979b-2662-329a-907c-602fdc69c3a3&system=d35404e4-40d3-59d6-3608-581c9495d86a
1416
"""
15-
...
17+
# Using data from lorenz.
18+
# Mesh file: /storage/shared/oceanparcels/input_data/MOi/domain_ORCA0083-N006/domain_ORCA0083-N006/PSY4V3R1_mesh_hgr.nc
19+
# Data files: /storage/shared/oceanparcels/input_data/MOi/GLO12/psy4v3r1-daily_{U,V}_*.nc
20+
# used modulefile for reference: "/storage/shared/oceanparcels/input_data/MOi/psy4v3r1/create_fieldset2D.py"
21+
22+
# scp "lorenz:/storage/shared/oceanparcels/input_data/MOi/GLO12/psy4v3r1-daily_{U,V,W,T}_2007-01-0{1,2}.nc" data-v4/nemo/field
23+
24+
time_counter_data = pd.date_range(start="2007-01-01T12:00:00", periods=T, freq="D")
25+
y_data = np.arange(1, Y + 1)
26+
x_data = np.arange(1, X + 1)
27+
deptht_data = np.linspace(0.494, 5.728e03, Z)
28+
29+
# Create the dataset
30+
return xr.Dataset(
31+
data_vars={
32+
"sotkeavmu1": (
33+
("time_counter", "y", "x"),
34+
np.random.rand(T, Y, X).astype(np.float64),
35+
{
36+
"units": "m2 s-1",
37+
"valid_min": np.float64(0.0),
38+
"valid_max": np.float64(100.0),
39+
"long_name": "Vertical Eddy Viscosity U 1m",
40+
"standard_name": "ocean_vertical_eddy_viscosity_u_1m",
41+
"short_name": "sotkeavmu1",
42+
"online_operation": "N/A",
43+
"interval_operation": np.int64(86400),
44+
"interval_write": np.int64(86400),
45+
"associate": "time_counter nav_lat nav_lon",
46+
},
47+
),
48+
"sotkeavmu15": (
49+
("time_counter", "y", "x"),
50+
np.random.rand(T, Y, X).astype(np.float64),
51+
{
52+
"units": "m2 s-1",
53+
"valid_min": np.float64(0.0),
54+
"valid_max": np.float64(100.0),
55+
"long_name": "Vertical Eddy Viscosity U 15m",
56+
"standard_name": "ocean_vertical_eddy_viscosity_u_15m",
57+
"short_name": "sotkeavmu15",
58+
"online_operation": "N/A",
59+
"interval_operation": np.int64(86400),
60+
"interval_write": np.int64(86400),
61+
"associate": "time_counter nav_lat nav_lon",
62+
},
63+
),
64+
"sotkeavmu30": (
65+
("time_counter", "y", "x"),
66+
np.random.rand(T, Y, X).astype(np.float64),
67+
{
68+
"units": "m2 s-1",
69+
"valid_min": np.float64(0.0),
70+
"valid_max": np.float64(100.0),
71+
"long_name": "Vertical Eddy Viscosity U 30m",
72+
"standard_name": "ocean_vertical_eddy_viscosity_u_30m",
73+
"short_name": "sotkeavmu30",
74+
"online_operation": "N/A",
75+
"interval_operation": np.int64(86400),
76+
"interval_write": np.int64(86400),
77+
"associate": "time_counter nav_lat nav_lon",
78+
},
79+
),
80+
"sotkeavmu50": (
81+
("time_counter", "y", "x"),
82+
np.random.rand(T, Y, X).astype(np.float64),
83+
{
84+
"units": "m2 s-1",
85+
"valid_min": np.float64(0.0),
86+
"valid_max": np.float64(100.0),
87+
"long_name": "Vertical Eddy Viscosity U 50m",
88+
"standard_name": "ocean_vertical_eddy_viscosity_u_50m",
89+
"short_name": "sotkeavmu50",
90+
"online_operation": "N/A",
91+
"interval_operation": np.int64(86400),
92+
"interval_write": np.int64(86400),
93+
"associate": "time_counter nav_lat nav_lon",
94+
},
95+
),
96+
"vozocrtx": (
97+
("time_counter", "deptht", "y", "x"),
98+
np.random.rand(T, Z, Y, X).astype(np.float64),
99+
{
100+
"units": "m s-1",
101+
"valid_min": np.float64(-10.0),
102+
"valid_max": np.float64(10.0),
103+
"long_name": "Zonal velocity",
104+
"standard_name": "sea_water_x_velocity",
105+
"short_name": "vozocrtx",
106+
"online_operation": "N/A",
107+
"interval_operation": np.int64(86400),
108+
"interval_write": np.int64(86400),
109+
"associate": "time_counter deptht nav_lat nav_lon",
110+
},
111+
),
112+
},
113+
coords={
114+
"nav_lon": (
115+
("y", "x"),
116+
np.random.rand(Y, X).astype(np.float32),
117+
{
118+
"units": "degrees_east",
119+
"valid_min": np.float32(-179.99984754002182),
120+
"valid_max": np.float32(179.999842386314),
121+
"long_name": "Longitude",
122+
"nav_model": "Default grid",
123+
"standard_name": "longitude",
124+
},
125+
),
126+
"nav_lat": (
127+
("y", "x"),
128+
np.random.rand(Y, X).astype(np.float32),
129+
{
130+
"units": "degrees_north",
131+
"valid_min": np.float32(-77.0104751586914),
132+
"valid_max": np.float32(89.9591064453125),
133+
"long_name": "Latitude",
134+
"nav_model": "Default grid",
135+
"standard_name": "latitude",
136+
},
137+
),
138+
"x": (("x",), x_data, {"standard_name": "projection_x_coordinate", "axis": "X", "units": "1"}),
139+
"y": (("y",), y_data, {"standard_name": "projection_y_coordinate", "axis": "Y", "units": "1"}),
140+
"time_counter": (
141+
("time_counter",),
142+
time_counter_data,
143+
{"standard_name": "time", "long_name": "Time axis", "axis": "T", "time_origin": "1950-JAN-01 00:00:00"},
144+
),
145+
"deptht": (
146+
("deptht",),
147+
deptht_data,
148+
{
149+
"units": "m",
150+
"positive": "down",
151+
"valid_min": np.float64(0.4940253794193268),
152+
"valid_max": np.float64(5727.91650390625),
153+
"long_name": "Vertical T levels",
154+
"standard_name": "depth",
155+
"axis": "Z",
156+
},
157+
),
158+
},
159+
attrs={
160+
"Conventions": "CF-1.0",
161+
"file_name": "ORCA12_LIM-T00_y2021m09d27_gridU.nc",
162+
"institution": "MERCATOR OCEAN",
163+
"source": "NEMO",
164+
"TimeStamp": "2021-OCT-03 18:27:01 GMT-0000",
165+
"references": "http://www.mercator-ocean.eu",
166+
},
167+
)
16168

17169

18170
def _hycom_data() -> xr.Dataset:

‎parcels/_datasets/utils.py‎

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -65,3 +65,95 @@ def dataset_repr_diff(ds1: xr.Dataset, ds2: xr.Dataset) -> str:
6565

6666
diff = difflib.ndiff(repr1.splitlines(keepends=True), repr2.splitlines(keepends=True))
6767
return "".join(diff)
68+
69+
70+
def compare_datasets(ds1, ds2, ds1_name="Dataset 1", ds2_name="Dataset 2"):
71+
print(f"Comparing {ds1_name} and {ds2_name}\n")
72+
73+
# Compare dataset attributes
74+
print("Dataset Attributes Comparison:")
75+
if ds1.attrs == ds2.attrs:
76+
print(" Dataset attributes are identical.")
77+
else:
78+
print(" Dataset attributes differ.")
79+
for attr_name in set(ds1.attrs.keys()) | set(ds2.attrs.keys()):
80+
if attr_name not in ds1.attrs:
81+
print(f" Attribute '{attr_name}' only in {ds2_name}")
82+
elif attr_name not in ds2.attrs:
83+
print(f" Attribute '{attr_name}' only in {ds1_name}")
84+
elif ds1.attrs[attr_name] != ds2.attrs[attr_name]:
85+
print(f" Attribute '{attr_name}' differs:")
86+
print(f" {ds1_name}: {ds1.attrs[attr_name]}")
87+
print(f" {ds2_name}: {ds2.attrs[attr_name]}")
88+
print("-" * 30)
89+
90+
# Compare dimensions
91+
print("Dimensions Comparison:")
92+
ds1_dims = set(ds1.dims)
93+
ds2_dims = set(ds2.dims)
94+
if ds1_dims == ds2_dims:
95+
print(" Dimension names are identical.")
96+
else:
97+
print(" Dimension names differ:")
98+
print(f" {ds1_name} dims: {sorted(list(ds1_dims))}")
99+
print(f" {ds2_name} dims: {sorted(list(ds2_dims))}")
100+
101+
# For common dimensions, compare order (implicit by comparing coordinate values for sortedness)
102+
# and size (though size is parameterized and expected to be different)
103+
for dim_name in ds1_dims.intersection(ds2_dims):
104+
print(f" Dimension '{dim_name}':")
105+
# Sizes will differ due to DIM_SIZE, so we don't strictly compare them.
106+
print(f" {ds1_name} size: {ds1.dims[dim_name]}, {ds2_name} size: {ds2.dims[dim_name]}")
107+
# Check if coordinates associated with dimensions are sorted (increasing)
108+
if dim_name in ds1.coords and dim_name in ds2.coords:
109+
is_ds1_sorted = np.all(np.diff(ds1[dim_name].values) >= 0) if len(ds1[dim_name].values) > 1 else True
110+
is_ds2_sorted = np.all(np.diff(ds2[dim_name].values) >= 0) if len(ds2[dim_name].values) > 1 else True
111+
if is_ds1_sorted == is_ds2_sorted:
112+
print(f" Order for '{dim_name}' is consistent (both sorted: {is_ds1_sorted})")
113+
else:
114+
print(
115+
f" Order for '{dim_name}' differs: {ds1_name} sorted: {is_ds1_sorted}, {ds2_name} sorted: {is_ds2_sorted}"
116+
)
117+
print("-" * 30)
118+
119+
# Compare variables (name, attributes, dimensions used)
120+
print("Variables Comparison:")
121+
ds1_vars = set(ds1.variables.keys())
122+
ds2_vars = set(ds2.variables.keys())
123+
124+
if ds1_vars == ds2_vars:
125+
print(" Variable names are identical.")
126+
else:
127+
print(" Variable names differ:")
128+
print(f" {ds1_name} vars: {sorted(list(ds1_vars - ds2_vars))}")
129+
print(f" {ds2_name} vars: {sorted(list(ds2_vars - ds1_vars))}")
130+
print(f" Common vars: {sorted(list(ds1_vars.intersection(ds2_vars)))}")
131+
132+
for var_name in ds1_vars.intersection(ds2_vars):
133+
print(f" Variable '{var_name}':")
134+
var1 = ds1[var_name]
135+
var2 = ds2[var_name]
136+
137+
# Compare attributes
138+
if var1.attrs == var2.attrs:
139+
print(" Attributes are identical.")
140+
else:
141+
print(" Attributes differ.")
142+
for attr_name in set(var1.attrs.keys()) | set(var2.attrs.keys()):
143+
if attr_name not in var1.attrs:
144+
print(f" Attribute '{attr_name}' only in {ds2_name}'s '{var_name}'")
145+
elif attr_name not in var2.attrs:
146+
print(f" Attribute '{attr_name}' only in {ds1_name}'s '{var_name}'")
147+
elif var1.attrs[attr_name] != var2.attrs[attr_name]:
148+
print(f" Attribute '{attr_name}' differs for '{var_name}':")
149+
print(f" {ds1_name}: {var1.attrs[attr_name]}")
150+
print(f" {ds2_name}: {var2.attrs[attr_name]}")
151+
152+
# Compare dimensions used by the variable
153+
if var1.dims == var2.dims:
154+
print(f" Dimensions used are identical: {var1.dims}")
155+
else:
156+
print(" Dimensions used differ:")
157+
print(f" {ds1_name}: {var1.dims}")
158+
print(f" {ds2_name}: {var2.dims}")
159+
print("=" * 30 + " End of Comparison " + "=" * 30)

0 commit comments

Comments
 (0)