In script plot_snapshot.py, add arguments to limit plot window and to

plot the grid. Incorporate the functionality of script
plot_snapshot_light.py into script plot_snapshot.py, with option
--light, to avoid duplication of source code.

Tests/plot_snapshot.py

@@ -9,7 +9,8 @@ contour. Crosses for a valid outermost contour but with no distinct
 max-speed contour. Filled circles for a valid outermost contour with a
 distinct max-speed contour.
 
-This scripts takes about 30 s of CPU for a 90° by 90° window.
+This scripts takes about 30 s of CPU for a 90° by 90° window, about 4
+mn for a 180° by 180° window.
 
 """
 
@@ -28,33 +29,44 @@ 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("-g", "--grid", help = "plot grid", action = "store_true")
+parser.add_argument("-w", "--window", help = "choose a limited plot window",
+                    action = "store_true")
+parser.add_argument("-l", "--light", help = "lighter plot",
+                    action = "store_true")
 args = parser.parse_args()
 
 with netCDF4.Dataset("h.nc") as f:
     longitude = f.variables["lon"][:]
     latitude = f.variables["lat"][:]
     
-llcrnrlon = float(input("llcrnrlon = ? "))
-llcrnrlat = float(input("llcrnrlat = ? "))
-urcrnrlon = float(input("urcrnrlon = ? "))
-urcrnrlat = float(input("urcrnrlat = ? "))
+if args.window:
+    llcrnrlon = float(input("llcrnrlon = ? "))
+    llcrnrlat = float(input("llcrnrlat = ? "))
+    urcrnrlon = float(input("urcrnrlon = ? "))
+    urcrnrlat = float(input("urcrnrlat = ? "))
 
-if urcrnrlon - llcrnrlon > 360:
-    sys.exit("bad values of urcrnrlon and llcrnrlon")
+    if urcrnrlon - llcrnrlon > 360:
+        sys.exit("bad values of urcrnrlon and llcrnrlon")
 
-longitude += np.ceil((llcrnrlon - longitude) / 360) * 360
-# (in [llcrnrlon, llcrnrlon + 2 pi[)
+    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)
-lon_2d, lat_2d = np.meshgrid(longitude[lon_mask], latitude[lat_mask])
+    lon_mask = longitude <= urcrnrlon
+    lat_mask = np.logical_and(latitude >= llcrnrlat, latitude <= urcrnrlat)
+else:
+    lon_mask = True
+    lat_mask = True
 
 plt.figure()
 src_crs = ccrs.PlateCarree()
 projection = ccrs.PlateCarree()
 ax = plt.axes(projection = projection)
-ax.plot(lon_2d.reshape(-1), lat_2d.reshape(-1), transform = src_crs,
-        marker = "+", color = "gray", linestyle = "None")
+
+if args.grid:
+    lon_2d, lat_2d = np.meshgrid(longitude[lon_mask], latitude[lat_mask])
+    ax.plot(lon_2d.reshape(-1), lat_2d.reshape(-1), transform = src_crs,
+            marker = "+", color = "gray", linestyle = "None")
 
 reader_extr = shapefile.Reader("extremum_1")
 reader_outer = shapefile.Reader("outermost_contour_1")
@@ -73,10 +85,13 @@ for shape_rec_extr, shape_outer, shape_m_s \
     in zip(reader_extr.iterShapeRecords(), reader_outer.iterShapes(),
            m_s_iterShapes):
     points = shape_rec_extr.shape.points[0]
-    points[0] += math.ceil((llcrnrlon - points[0]) / 360) * 360
-    # (in [llcrnrlon, llcrnrlon + 2 pi[)
 
-    if points[0] <= urcrnrlon and llcrnrlat <= points[1] <= urcrnrlat:
+    if args.window:
+        points[0] += math.ceil((llcrnrlon - points[0]) / 360) * 360
+        # (in [llcrnrlon, llcrnrlon + 2 pi[)
+
+    if not args.window or (points[0] <= urcrnrlon 
+                          and llcrnrlat <= points[1] <= urcrnrlat):
         if shape_rec_extr.record[4] == 0:
                 # Anti-cyclone
             color = "red"
@@ -86,34 +101,48 @@ for shape_rec_extr, shape_outer, shape_m_s \
         lines = ax.plot(points[0], points[1], markersize = 10, color = color,
                         fillstyle = "none", transform = src_crs)
 
-        if shape_rec_extr.record[5] == 0:
+        if shape_rec_extr.record[5] != 0:
+            if args.light:
+                lines[0].set_marker("+")
+            else:
+                lines[0].set_marker("o")
+        elif not args.light:
             # Invalid outermost contour
             lines[0].set_marker("s")
-        else:
-            lines[0].set_marker("o")
 
-        ax.annotate(str(shape_rec_extr.record[2]),
-                    projection.transform_point(points[0], points[1], src_crs),
-                    xytext = (3, 3), textcoords = "offset points")
+        if not args.light:
+            ax.annotate(str(shape_rec_extr.record[2]),
+                        projection.transform_point(points[0], points[1],
+                                                   src_crs),
+                        xytext = (3, 3), textcoords = "offset points")
 
         if shape_outer.shapeType != shapefile.NULL:
             points = np.array(shape_outer.points)
             lines = ax.plot(points[:, 0], points[:, 1], color = color,
                             transform = src_crs)
 
-            if shape_rec_extr.record[5] == 0:
-                # Invalid outermost contour
-                lines[0].set_marker("s")
-                lines[0].set_fillstyle("none")
-            elif shape_m_s.shapeType == shapefile.NULL:
-                lines[0].set_marker("x")
-            else:
-                lines[0].set_marker("o")
+            if not args.light:
+                if shape_rec_extr.record[5] == 0:
+                    # Invalid outermost contour
+                    lines[0].set_marker("s")
+                    lines[0].set_fillstyle("none")
+                elif shape_m_s.shapeType == shapefile.NULL:
+                    lines[0].set_marker("x")
+                else:
+                    lines[0].set_marker("o")
 
             if shape_m_s.shapeType != shapefile.NULL:
                 points = np.array(shape_m_s.points)
-                ax.plot(points[:, 0], points[:, 1], color = color, marker = "o",
-                        transform = src_crs)
+                
+                if shape_rec_extr.record[4] == 0:
+                    # Anti-cyclone
+                    color = "magenta"
+                else:
+                    color = "cyan"
+
+                lines = ax.plot(points[:, 0], points[:, 1], color = color,
+                                transform = src_crs)
+                if not args.light: lines[0].set_marker("o")
 
 if args.velocity:
     with netCDF4.Dataset("uv.nc") as f:

Tests/plot_snapshot_light.py

-#!/usr/bin/env python3
-
-"""Plots outermost contours and max-speed contours.
-
-Red or magenta for anti-cyclones, blue or cyan for cyclones. Only
-extremums with outermost contour of sufficient amplitude.
-
-"""
-
-import shapefile
-import numpy as np
-import matplotlib.pyplot as plt
-import cartopy.crs as ccrs
-
-reader = shapefile.Reader("extremum_1")
-reader_outer = shapefile.Reader("outermost_contour_1")
-reader_m_s = shapefile.Reader("max_speed_contour_1")
-
-plt.figure()
-src_crs = ccrs.PlateCarree()
-projection = ccrs.PlateCarree()
-ax = plt.axes(projection = projection)
-##ax = plt.axes()
-
-for shape_rec_extr, shape_outer, shape_m_s in zip(reader.iterShapeRecords(),
-                                                  reader_outer.iterShapes(),
-                                                  reader_m_s.iterShapes()):
-    if shape_rec_extr.record[5] == 1:
-        # Valid outermost contour
-        points = shape_rec_extr.shape.points[0]
-
-        if shape_rec_extr.record[4] == 0:
-            # Anti-cyclone
-            color = "red"
-        else:
-            color = "blue"
-        
-        ax.plot(points[0], points[1], "+", markersize = 10, color = color,
-                fillstyle = "none", transform = src_crs)
-
-        points = np.array(shape_outer.points)
-        ax.plot(points[:, 0], points[:, 1], color = color, transform = src_crs)
-
-        if shape_rec_extr.record[4] == 0:
-            # Anti-cyclone
-            color = "magenta"
-        else:
-            color = "cyan"
-        
-        if shape_m_s.shapeType != shapefile.NULL:
-            points = np.array(shape_m_s.points)
-            ax.plot(points[:, 0], points[:, 1], color = color,
-                    transform = src_crs)
-
-ax.set_xlabel("longitude (degrees east)")
-ax.set_ylabel("latitude (degrees north)")
-ax.gridlines(draw_labels = True)
-ax.coastlines()
-
-plt.show()