How to plot Ocean Currents with Cartopy

Question

I am trying to plot a netCDF4 file containing ocean currents from a NASA database for a project, but I keep getting errors such as "x and y coordinates are not compatible with the shape of the vector components".

I have tried changing the streamplot to a contourf (when I did it said that it needed to be a 2d array) which I tried to change but I could not get it to work.

import os
import matplotlib.pyplot as plt
from netCDF4 import Dataset as netcdf_dataset
import numpy as np
import cartopy.crs as ccrs

fname = "oscar_vel2019.nc.gz.nc"

data=netcdf_dataset(fname)
v = data.variables['v'][0, :, :, :]
vf = data.variables['vm'][0, :, :, :]
u = data.variables['u'][0, :, :, :]
uf = data.variables['um'][0, :, :, :]
lats = data.variables['latitude'][:]
lons = data.variables['longitude'][:]
ax = plt.axes(projection=ccrs.PlateCarree())

mymap=plt.streamplot(lons, lats, u, v, 60, transform=ccrs.PlateCarree())

ax.coastlines()

plt.show()

I would like it to work such that the ocean currents are visible on the plot and to show the movement of particles in the currents through an animation. I really don't have much knowledge with this which is why I am asking. Here is the link from which I got the file: https://podaac-opendap.jpl.nasa.gov/opendap/hyrax/allData/oscar/preview/L4/oscar_third_deg/oscar_vel2019.nc.gz.html

Answer 1

OK, I downloaded the data. The problem is that u and v are 4-dimensional, so you need to squeeze out the "depth" dimension. Cartopy also doesn't accept longitudes greater than 180, and you probably won't get away with stream plotting the whole thing. Also, density=60 will take forever...

This is ugly, but gives you the idea.

import xarray as xr
import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

with xr.open_dataset('/Users/jklymak/downloads/oscar_vel2019.nc.gz.nc') as ds:
    print(ds)

    ax = plt.axes(projection=ccrs.PlateCarree())

    dec = 10
    lon = ds.longitude.values[::dec]
    lon[lon>180] = lon[lon>180] - 360
    mymap=plt.streamplot(lon, ds.latitude.values[::dec], ds.u.values[0, 0, ::dec, ::dec], ds.v.values[0, 0, ::dec, ::dec], 6, transform=ccrs.PlateCarree())
    ax.coastlines()
    plt.show()