Read variables from a model along the trajectories of drifters.

import xarray as xr
import cf_xarray as _
import pyproj
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

Open the Norkyst800 model

nk = xr.open_dataset(
    'https://thredds.met.no/thredds/dodsC/sea/norkyst800m/1h/aggregate_be')
print(nk)

<xarray.Dataset> Size: 89TB
Dimensions:                  (s_rho: 35, s_w: 36, Y: 902, X: 2602, time: 71509,
                              depth: 16)
Coordinates:
  * X                        (X) float64 21kB 0.0 800.0 ... 2.08e+06 2.081e+06
  * Y                        (Y) float64 7kB 0.0 800.0 ... 7.2e+05 7.208e+05
  * depth                    (depth) float64 128B 0.0 3.0 10.0 ... 2e+03 3e+03
  * time                     (time) datetime64[ns] 572kB 2017-02-20 ... 2025-...
    lat                      (Y, X) float64 19MB ...
    lon                      (Y, X) float64 19MB ...
Dimensions without coordinates: s_rho, s_w
Data variables: (12/20)
    Cs_r                     (s_rho) float64 280B ...
    Cs_w                     (s_w) float64 288B ...
    forecast_reference_time  datetime64[ns] 8B ...
    hc                       float64 8B ...
    projection_stere         int32 4B ...
    angle                    (Y, X) float64 19MB ...
    ...                       ...
    ubar                     (time, Y, X) float32 671GB ...
    v                        (time, depth, Y, X) float32 11TB ...
    v_northward              (time, depth, Y, X) float32 11TB ...
    vbar                     (time, Y, X) float32 671GB ...
    w                        (time, depth, Y, X) float32 11TB ...
    zeta                     (time, Y, X) float32 671GB ...
Attributes: (12/68)
    file:                       /home/metno_op/run/norkyst-800m_2017/ocean_hi...
    type:                       ROMS/TOMS history file
    var_info:                   /home/metno_op/sea/ROMS/metroms_apps/norkyst-...
    rst_file:                   /home/metno_op/run/norkyst-800m_2017/ocean_rs...
    his_file:                   /home/metno_op/run/norkyst-800m_2017/ocean_hi...
    avg_file:                   /home/metno_op/run/norkyst-800m_2017/ocean_av...
    ...                         ...
    publisher_name:             Norwegian Meteorological Institute
    publisher_institution:      Norwegian Meteorological Institute
    publisher_email:            data-management-group@met.no
    publisher_url:              https://www.met.no
    project:                    Ocean and Ice - Research to Operation (HI-R2O)
    license:                    https://spdx.org/licenses/CC-BY-4.0 (CC-BY-4.0)

Use cf-xarray to get the CRS

gm = nk.cf['grid_mapping']
nk_crs = pyproj.CRS.from_cf(gm.attrs)
print(nk_crs)

# Converting to Cartopy is not possible to do automatically, unfortunately, so we define it manually:
ncrs = ccrs.Stereographic(true_scale_latitude=60,
                          central_latitude=90,
                          central_longitude=70,
                          false_easting=3192800,
                          false_northing=1784000)

{"$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"}]}}

Plot northward surface current

plt.axes(projection=ncrs)
# plt.axes(projection=ccrs.Mercator())
nk.isel(time=-1, depth=0).v_northward.plot.imshow(
    transform=ncrs)  # TODO: SUPPLY CORRECT TRANSFORM HERE!

plt.scatter(5.32,
            60.4,
            10,
            transform=ccrs.PlateCarree(),
            marker='*',
            label='Bergen')

plt.legend()
plt.show()