In script plot_snapshot.py, bug fix for type of argument scale. Add

the possibility to restrict the plotted domain because the script
takes a very long time for a global domain.

Tests/plot_snapshot.py

@@ -9,6 +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.
+
 """
 
 import argparse
@@ -19,11 +21,12 @@ import netCDF4
 import math
 import itertools
 import cartopy.crs as ccrs
+import sys
 
 parser = argparse.ArgumentParser()
 parser.add_argument("-v", "--velocity", help = "plot velocity field",
                     action = "store_true")
-parser.add_argument("-s", "--scale", default = 20,
+parser.add_argument("-s", "--scale", default = 20, type = float,
                     help = "scale of arrows for the velocity field")
 args = parser.parse_args()
 
@@ -31,11 +34,20 @@ with netCDF4.Dataset("h.nc") as f:
     longitude = f.variables["lon"][:]
     latitude = f.variables["lat"][:]
     
-lon_2d, lat_2d = np.meshgrid(longitude, latitude)
-# llcrnrlon = math.floor(longitude[0])
-# llcrnrlat = math.floor(latitude[0])
-# urcrnrlon = math.ceil(longitude[-1])
-# urcrnrlat = math.ceil(latitude[-1])
+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")
+
+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])
 
 plt.figure()
 src_crs = ccrs.PlateCarree()
@@ -61,51 +73,56 @@ 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]
-    
-    if shape_rec_extr.record[4] == 0:
-            # Anti-cyclone
-        color = "red"
-    else:
-        color = "blue"
-        
-    lines = ax.plot(points[0], points[1], markersize = 10, color = color,
-                    fillstyle = "none", transform = src_crs)
-
-    if shape_rec_extr.record[5] == 0:
-        # 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 shape_outer.shapeType != shapefile.NULL:
-        points = np.array(shape_outer.points)
-        lines = ax.plot(points[:, 0], points[:, 1], color = color,
-                        transform = src_crs)
+    points[0] += math.ceil((llcrnrlon - points[0]) / 360) * 360
+    # (in [llcrnrlon, llcrnrlon + 2 pi[)
+
+    if points[0] <= urcrnrlon and llcrnrlat <= points[1] <= urcrnrlat:
+        if shape_rec_extr.record[4] == 0:
+                # Anti-cyclone
+            color = "red"
+        else:
+            color = "blue"
+
+        lines = ax.plot(points[0], points[1], markersize = 10, color = color,
+                        fillstyle = "none", 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 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)
+        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 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 args.velocity:
     with netCDF4.Dataset("uv.nc") as f:
-        quiver_return = ax.quiver(longitude, latitude, f["u"][0], f["v"][0],
+        quiver_return = ax.quiver(longitude[lon_mask], latitude[lat_mask],
+                                  f["u"][0, lat_mask, lon_mask],
+                                  f["v"][0, lat_mask, lon_mask],
                                   scale = args.scale, scale_units = "width",
                                   transform = src_crs)
-    plt.quiverkey(quiver_return, 0.75, 0.9, 1, r"1 m s$^{-1}$",
+    plt.quiverkey(quiver_return, 0.9, 0.9, 1, r"1 m s$^{-1}$",
                   coordinates = "figure")
 
 ax.gridlines(draw_labels = True)