Use cartopy instead of basemap

Inst_eddies/Tests/test_mean_speed.py

 #!/usr/bin/env python3
 
 import netCDF4
-from mpl_toolkits import basemap
+import cartopy.crs as ccrs
 from scipy import interpolate
 import shapefile
 import matplotlib.pyplot as plt
 import numpy as np
 import sys
 import math
+from numpy import ma
 
 with netCDF4.Dataset("uv.nc") as f:
     u = f.variables["u"][:]
@@ -44,12 +45,9 @@ lat_0 = float(input("latitude extremum = ? "))
 lon_0_rad = lon_0 * np.pi/180
 lat_0_rad = lat_0 * np.pi/180
 
-m = basemap.Basemap(projection="stere", lon_0 = lon_0, lat_0 = lat_0,
-                    llcrnrlon = lon[0], llcrnrlat = lat[0],
-                    urcrnrlon = lon[- 1], urcrnrlat= lat[- 1])
-u_rot, v_rot = m.rotate_vector(u, v, lon, lat)
-ui_rot, vi_rot = m.rotate_vector(ui, vi, lon_interp[:- 1], lat_interp[:- 1])
-
+src_crs = ccrs.PlateCarree()
+projection = ccrs.Stereographic(central_latitude = lat_0,
+                                central_longitude = lon_0)
 x = np.cos(lat_0_rad) * (lon_interp_rad - lon_0_rad)
 y = lat_interp_rad  - lat_0_rad
 
@@ -57,23 +55,32 @@ v_azim = (x * vi - y * ui) / np.sqrt(x**2 + y**2)
 print("mean azimuthal speed = ", np.mean(v_azim), "m s-1")
 
 plt.figure()
-m.plot(lon_interp, lat_interp, latlon=True)
-m.drawparallels(range(math.floor(lat[0]), math.ceil(lat[- 1])),
-                labels = [1, 0, 0, 0])
-m.drawmeridians(range(math.floor(lon[0]), math.ceil(lon[- 1])),
-                labels = [0, 0, 0, 1])
-m.quiver(lon_interp[:- 1], lat_interp[:- 1], ui_rot, vi_rot, latlon = True)
-m.scatter(lon_0, lat_0, latlon = True)
+ax = plt.axes(projection = projection)
+ax.plot(lon_interp, lat_interp, transform = src_crs)
+ax.gridlines(ylocs = range(math.floor(lat[0]), math.ceil(lat[- 1])),
+             xlocs = range(math.floor(lon[0]), math.ceil(lon[- 1])))
+u_src_crs = ui / np.cos(lat_interp[:- 1, np.newaxis] / 180 * np.pi)
+v_src_crs = vi
+magnitude = ma.sqrt(ui**2 + vi**2)
+magn_src_crs = ma.sqrt(u_src_crs**2 + v_src_crs**2)
+ax.quiver(lon_interp[:- 1], lat_interp[:- 1],
+          u_src_crs * magnitude / magn_src_crs,
+          v_src_crs * magnitude / magn_src_crs, transform = src_crs)
+ax.scatter(lon_0, lat_0, transform = src_crs)
 plt.title("interpolated velocity")
 
 plt.figure()
-m.plot(lon_interp, lat_interp, latlon=True)
-m.drawparallels(range(math.floor(lat[0]), math.ceil(lat[- 1])),
-                labels = [1, 0, 0, 0])
-m.drawmeridians(range(math.floor(lon[0]), math.ceil(lon[- 1])),
-                labels = [0, 0, 0, 1])
-m.quiver(lon2d, lat2d, u_rot, v_rot, latlon = True)
-m.scatter(lon_0, lat_0, latlon = True)
+ax.plot(lon_interp, lat_interp, transform = src_crs)
+ax.gridlines(ylocs = range(math.floor(lat[0]), math.ceil(lat[- 1])),
+             xlocs = range(math.floor(lon[0]), math.ceil(lon[- 1])))
+ax = plt.axes(projection = projection)
+u_src_crs = u / np.cos(lat[:, np.newaxis] / 180 * np.pi)
+v_src_crs = v
+magnitude = ma.sqrt(u**2 + v**2)
+magn_src_crs = ma.sqrt(u_src_crs**2 + v_src_crs**2)
+ax.quiver(lon2d, lat2d, u_src_crs * magnitude / magn_src_crs,
+          v_src_crs * magnitude / magn_src_crs, transform = src_crs)
+ax.scatter(lon_0, lat_0, transform = src_crs)
 plt.title("velocity at grid points")
 
 plt.show()