Note
Go to the end to download the full example code.
Read variables from a model along the trajectories of drifters.
import xarray as xr
import numpy as np
import cf_xarray as _
import pyproj
import trajan as ta
import matplotlib.pyplot as plt
Open drifter dataset from a CSV file
ds = ta.read_csv('bug05_pos.csv.xz',
lon='Longitude',
lat='Latitude',
time='Time',
name='Device')
print(ds)
2024-09-02 11:26:10 fv-az1272-597 trajan.trajectory_accessor[3440] DEBUG Detecting time-dimension for "obs"..
2024-09-02 11:26:10 fv-az1272-597 trajan.trajectory_accessor[3440] DEBUG Detected obs-dim: obs, detected time-dim: time.
2024-09-02 11:26:10 fv-az1272-597 trajan.trajectory_accessor[3440] DEBUG Detected un-structured (2D) trajectory dataset
2024-09-02 11:26:10 fv-az1272-597 trajan.traj2d[3440] DEBUG Condensing 844 observations.
2024-09-02 11:26:10 fv-az1272-597 trajan.traj2d[3440] DEBUG Condensed observations from: 844 to 844
2024-09-02 11:26:10 fv-az1272-597 trajan.trajectory_accessor[3440] DEBUG Detecting time-dimension for "obs"..
2024-09-02 11:26:10 fv-az1272-597 trajan.trajectory_accessor[3440] DEBUG Detected obs-dim: obs, detected time-dim: time.
2024-09-02 11:26:10 fv-az1272-597 trajan.trajectory_accessor[3440] DEBUG Detected un-structured (2D) trajectory dataset
2024-09-02 11:26:10 fv-az1272-597 trajan.traj[3440] DEBUG No grid-mapping specified, checking if coordinates are lon/lat..
2024-09-02 11:26:10 fv-az1272-597 trajan.traj[3440] DEBUG No grid-mapping specified, checking if coordinates are lon/lat..
2024-09-02 11:26:10 fv-az1272-597 trajan.traj[3440] DEBUG No grid-mapping specified, checking if coordinates are lon/lat..
2024-09-02 11:26:10 fv-az1272-597 trajan.traj[3440] DEBUG No grid-mapping specified, checking if coordinates are lon/lat..
<xarray.Dataset> Size: 27kB
Dimensions: (trajectory: 1, obs: 844)
Coordinates:
* trajectory (trajectory) <U18 72B 'dev864475040536665'
* obs (obs) int64 7kB 0 1 2 3 4 5 6 7 ... 837 838 839 840 841 842 843
Data variables:
time (trajectory, obs) datetime64[ns] 7kB 2022-05-10T12:56:16.5000...
lon (trajectory, obs) float64 7kB 5.332 5.332 5.332 ... 5.25 5.25
lat (trajectory, obs) float64 7kB 60.38 60.38 60.38 ... 60.45 60.45
Attributes:
Conventions: CF-1.10
featureType: trajectory
geospatial_lat_min: 60.38380600000001
geospatial_lat_max: 60.45141000000001
geospatial_lon_min: 5.244493666666667
geospatial_lon_max: 5.33242
time_coverage_start: 2020-01-01T00:00:07.250000
time_coverage_end: 2022-05-20T11:03:46.120000
Open the Norkyst800 model
nk = xr.open_dataset(
'https://thredds.met.no/thredds/dodsC/sea/norkyst800m/1h/aggregate_be')
Use cf-xarray to get the CRS
gm = nk.cf['grid_mapping']
nk_crs = pyproj.CRS.from_cf(gm.attrs)
print(nk_crs)
2024-09-02 11:26:21 fv-az1272-597 pyproj[3440] DEBUG PROJ_ERROR: proj_create: several objects matching this name: Krovak (Greenwich), Equal Earth Greenwich, Laborde Grid (Greenwich), Modified Krovak (Greenwich), Krovak East North (Greenwich), Modified Krovak East North (Greenwich), ...
{"$schema": "https://proj.org/schemas/v0.2/projjson.schema.json", "type": "ProjectedCRS", "name": "undefined", "base_crs": {"$schema": "https://proj.org/schemas/v0.7/projjson.schema.json", "type": "GeographicCRS", "name": "undefined", "datum": {"type": "GeodeticReferenceFrame", "name": "undefined", "ellipsoid": {"name": "undefined", "semi_major_axis": 6378137, "semi_minor_axis": 6356752.3142}}, "coordinate_system": {"subtype": "ellipsoidal", "axis": [{"name": "Longitude", "abbreviation": "lon", "direction": "east", "unit": "degree"}, {"name": "Latitude", "abbreviation": "lat", "direction": "north", "unit": "degree"}]}}, "conversion": {"$schema": "https://proj.org/schemas/v0.7/projjson.schema.json", "type": "Conversion", "name": "unknown", "method": {"name": "Polar Stereographic (variant B)", "id": {"authority": "EPSG", "code": 9829}}, "parameters": [{"name": "Latitude of standard parallel", "value": 60, "unit": "degree", "id": {"authority": "EPSG", "code": 8832}}, {"name": "Longitude of origin", "value": 70, "unit": "degree", "id": {"authority": "EPSG", "code": 8833}}, {"name": "False easting", "value": 3192800, "unit": "metre", "id": {"authority": "EPSG", "code": 8806}}, {"name": "False northing", "value": 1784000, "unit": "metre", "id": {"authority": "EPSG", "code": 8807}}]}, "coordinate_system": {"$schema": "https://proj.org/schemas/v0.7/projjson.schema.json", "type": "CoordinateSystem", "subtype": "Cartesian", "axis": [{"name": "Easting", "abbreviation": "E", "direction": "east", "unit": "metre"}, {"name": "Northing", "abbreviation": "N", "direction": "north", "unit": "metre"}]}}
Grid the drifter dataset to the timesteps of the model.
times = nk.sel(time=slice('2022-05-10', '2022-05-20')).time.values
ds = ds.traj.gridtime(times)
ds = ds.dropna('time', how='all')
print(ds)
2024-09-02 11:26:21 fv-az1272-597 trajan.trajectory_accessor[3440] DEBUG Detecting time-dimension for "obs"..
2024-09-02 11:26:21 fv-az1272-597 trajan.trajectory_accessor[3440] DEBUG Detected obs-dim: obs, detected time-dim: time.
2024-09-02 11:26:21 fv-az1272-597 trajan.trajectory_accessor[3440] DEBUG Detected un-structured (2D) trajectory dataset
<xarray.Dataset> Size: 6kB
Dimensions: (trajectory: 1, time: 239)
Coordinates:
* time (time) datetime64[ns] 2kB 2022-05-10T13:00:00 ... 2022-05-20T...
* trajectory (trajectory) <U18 72B 'dev864475040536665'
Data variables:
lon (trajectory, time) float64 2kB 5.332 5.332 5.332 ... 5.249 5.25
lat (trajectory, time) float64 2kB 60.38 60.38 60.38 ... 60.44 60.45
Attributes:
Conventions: CF-1.10
featureType: trajectory
geospatial_lat_min: 60.38380600000001
geospatial_lat_max: 60.45141000000001
geospatial_lon_min: 5.244493666666667
geospatial_lon_max: 5.33242
time_coverage_start: 2020-01-01T00:00:07.250000
time_coverage_end: 2022-05-20T11:03:46.120000
Transform the drifter dataset to the CRS of the model
tx, ty = ds.traj.transform(nk_crs, ds.traj.tx, ds.traj.ty)
2024-09-02 11:26:21 fv-az1272-597 trajan.trajectory_accessor[3440] DEBUG Detected obs-dim: time, detected time-dim: time.
2024-09-02 11:26:21 fv-az1272-597 trajan.trajectory_accessor[3440] DEBUG Detected structured (1D) trajectory dataset
2024-09-02 11:26:21 fv-az1272-597 trajan.traj[3440] DEBUG No grid-mapping specified, checking if coordinates are lon/lat..
By making sure the coordinates has defined dimensions xarray can select along the dimensions, and does not return slices along all coordinates. See xarray slicing for more details.
tx = xr.DataArray(tx, dims=['trajectory', 'time'])
ty = xr.DataArray(ty, dims=['trajectory', 'time'])
Extract the values of a variable for the trajectory
temp = nk.isel(depth=0).sel(time=ds.time,
X=tx,
Y=ty,
method='nearest').temperature
print(temp)
<xarray.DataArray 'temperature' (time: 239, trajectory: 1)> Size: 956B
[239 values with dtype=float32]
Coordinates:
X (trajectory, time) float64 2kB 3.416e+05 3.416e+05 ... 3.464e+05
Y (trajectory, time) float64 2kB 4.344e+05 4.344e+05 ... 4.416e+05
depth float64 8B 0.0
* time (time) datetime64[ns] 2kB 2022-05-10T13:00:00 ... 2022-05-20T11:...
lat (trajectory, time) float64 2kB 60.38 60.38 60.38 ... 60.44 60.45
lon (trajectory, time) float64 2kB 5.33 5.33 5.33 ... 5.243 5.243 5.249
Dimensions without coordinates: trajectory
Attributes:
units: Celsius
time: ocean_time
grid: grid
location: face
field: temperature, scalar, series
grid_mapping: projection_stere
long_name: Sea water potential temperature
standard_name: sea_water_potential_temperature
_ChunkSizes: [ 1 1 17 2602]
Notice that the lat and lon variables from Norkyst match the original lat and lon from the dataset.
plt.figure()
temp.plot()
plt.show()
![depth = 0.0 [m]](../_images/sphx_glr_example_model_along_trajectory_001.png)
Total running time of the script: (1 minutes 25.720 seconds)