Extract plotting of velocity in separate script

The script `plot_eddy_contours.py` was becoming too complicated.

Inst_eddies/Analysis/plot_eddy_contours.py

@@ -153,10 +153,6 @@ if __name__ == "__main__":
     import netCDF4
 
     parser = argparse.ArgumentParser()
-    parser.add_argument("-v", "--velocity", help = "plot velocity field",
-                        action = "store_true")
-    parser.add_argument("-s", "--scale", default = 20, type = float,
-                        help = "scale of arrows for the velocity field")
     parser.add_argument("-d", "--date", type = int,
                         help = "date, as days since 1950-1-1")
     parser.add_argument("-g", "--grid", help = "plot grid",
@@ -175,7 +171,7 @@ if __name__ == "__main__":
                         help = "Save file to specified format")
     args = parser.parse_args()
 
-    if args.grid or args.velocity:
+    if args.grid:
         with netCDF4.Dataset("h.nc") as f:
             if "lon" in f.variables:
                 lon = "lon"
@@ -193,14 +189,14 @@ if __name__ == "__main__":
         if urcrnrlon - llcrnrlon > 360:
             sys.exit("bad values of urcrnrlon and llcrnrlon")
 
-        if args.grid or args.velocity:
+        if args.grid:
             longitude += np.ceil((llcrnrlon - longitude) / 360) * 360
             # (in [llcrnrlon, llcrnrlon + 2 pi[)
 
             lon_mask = longitude <= urcrnrlon
             lat_mask = np.logical_and(latitude >= llcrnrlat,
                                       latitude <= urcrnrlat)
-    elif args.grid or args.velocity:
+    elif args.grid:
         lon_mask = np.ones(len(longitude), dtype = bool)
         lat_mask = np.ones(len(latitude), dtype = bool)
 
@@ -218,18 +214,6 @@ if __name__ == "__main__":
 
     if args.window is None: plot_grid_bb(args.shpc_dir[0], ax)
 
-    if args.velocity:
-        with netCDF4.Dataset("uv.nc") as f:
-            quiver_return = ax.quiver(longitude[lon_mask],
-                                      latitude[lat_mask],
-                                      f["ugos"][0, lat_mask][:, lon_mask],
-                                      f["vgos"][0, lat_mask][:, lon_mask],
-                                      scale = args.scale,
-                                      scale_units = "width",
-                                      transform = src_crs)
-        plt.quiverkey(quiver_return, 0.9, 0.9, 1, r"1 m s$^{-1}$",
-                      coordinates = "figure")
-
     for shpc_dir in args.shpc_dir:
         handler = util_eddies.open_shpc(shpc_dir)
 

Inst_eddies/Analysis/plot_velocity.py

+#!/usr/bin/env python3
+
+import numpy as np
+import cartopy.crs as ccrs
+import sys
+import matplotlib.pyplot as plt
+import argparse
+import netCDF4
+
+parser = argparse.ArgumentParser()
+parser.add_argument("-s", "--scale", default = 20, type = float,
+                    help = "scale of arrows for the velocity field")
+parser.add_argument("-w", "--window", help = "choose a limited plot window",
+                    type = float, nargs = 4,
+                    metavar = ("llcrnrlon", "llcrnrlat", "urcrnrlon",
+                               "urcrnrlat"))
+parser.add_argument("--save", metavar = "format",
+                    help = "Save file to specified format")
+args = parser.parse_args()
+
+with netCDF4.Dataset("h.nc") as f:
+    if "lon" in f.variables:
+        lon = "lon"
+        lat = "lat"
+    else:
+        lon = "longitude"
+        lat = "latitude"
+
+    longitude = f[lon][:]
+    latitude = f[lat][:]
+
+if args.window is not None:
+    llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat = args.window
+
+    if urcrnrlon - llcrnrlon > 360:
+        sys.exit("bad values of urcrnrlon and llcrnrlon")
+
+    longitude += np.ceil((llcrnrlon - longitude) / 360) * 360
+    # (in [llcrnrlon, llcrnrlon + 2 pi[)
+
+    lon_mask = longitude <= urcrnrlon
+    lat_mask = np.logical_and(latitude >= llcrnrlat,
+                              latitude <= urcrnrlat)
+else:
+    lon_mask = np.ones(len(longitude), dtype = bool)
+    lat_mask = np.ones(len(latitude), dtype = bool)
+
+fig = plt.figure()
+src_crs = ccrs.PlateCarree()
+projection = ccrs.PlateCarree()
+##projection = ccrs.NorthPolarStereo()
+ax = plt.axes(projection = projection)
+
+with netCDF4.Dataset("uv.nc") as f:
+    quiver_return = ax.quiver(longitude[lon_mask],
+                              latitude[lat_mask],
+                              f["ugos"][0, lat_mask][:, lon_mask],
+                              f["vgos"][0, lat_mask][:, lon_mask],
+                              scale = args.scale,
+                              scale_units = "width",
+                              transform = src_crs)
+plt.quiverkey(quiver_return, 0.9, 0.9, 1, r"1 m s$^{-1}$",
+              coordinates = "figure")
+
+ax.gridlines(draw_labels = True)
+ax.coastlines()
+
+if args.save:
+    plt.savefig(f"plot_velocity.{args.save}")
+    print(f'Created "plot_velocity.{args.save}".')
+else:
+    plt.show()

Inst_eddies/Analysis/tests.json

@@ -59,14 +59,9 @@
 	"description": "With window option."
     },
     {
-	"title": "Plot_eddy_contours_vel",
+	"title": "Plot_velocity",
 	"command":
-	[
-	    "$src_dir/Inst_eddies/Analysis/plot_eddy_contours.py",
-	    "$tests_old_dir/Extraction_eddies_region_4/SHPC_anti",
-	    "--velocity", "--save=png"
-	],
-	"description": "With velocity option.",
+	["$src_dir/Inst_eddies/Analysis/plot_velocity.py", "--save=png"],
 	"required":
 	[
 	    ["$src_dir/Inst_eddies/Tests/Input/h_region_4.nc", "h.nc"],