New component twice of the derived type eddy, for diagnostic of

program behavior.

In plot_snapshot.py, no need to accumulate all content of shapefile
extremum_1, just process it record by record, as for other shapefiles.

Sources/derived_types.f

@@ -21,6 +21,7 @@ module derived_types
                       ! amplitude of ssh between extremum and
                       ! outermost contour
      logical interpolated
+     logical twice ! outermost contour computed twice
      
      integer:: delta_in = huge(0)
      ! Maximum difference in time index where there is a direct

Sources/dispatch_snapshot.f

@@ -22,7 +22,7 @@ contains
 
     TYPE(shpfileobject), intent(inout):: hshp_outermost
     ! shapefile outermost_contour_$m. We assume that the fields in the
-    ! DBF file are, in that order: area, ssh, date index, eddy index.
+    ! DBF file are, in that order: area, ssh, date index, eddy index, twice.
 
     TYPE(shpfileobject), intent(inout):: hshp_max_speed
     ! shapefile x_speed_contour_$m. We assume that the fields in the

Sources/get_snapshot.f

@@ -147,6 +147,8 @@ contains
             noise_around(i) &
             = any(s%extr_map(llc(1, i):urc(1, i), llc(2, i):urc(2, i)) < 0)
 
+       s%list_vis%twice = flat_extr .and. s%list_vis%suff_amp .and. noise_around
+
        do i = 1, s%number_vis_eddies
           if (s%list_vis(i)%suff_amp .and. noise_around(i) &
                .or. .not. flat_extr(i)) then

Sources/write_eddy.f

@@ -23,7 +23,7 @@ contains
 
     TYPE(shpfileobject), intent(inout):: hshp_outermost
     ! shapefile outermost_contour_$m. We assume that the fields in the
-    ! DBF file are, in that order: area, ssh, date index, eddy index.
+    ! DBF file are, in that order: area, ssh, date index, eddy index, twice.
 
     TYPE(shpfileobject), intent(inout):: hshp_max_speed
     ! shapefile x_speed_contour_$m. We assume that the fields in the
@@ -72,6 +72,7 @@ contains
          e%outermost_contour%ssh)
     call dbf_write_attribute_03(hshp_outermost, ishape, 2, k)
     call dbf_write_attribute_03(hshp_outermost, ishape, 3, i)
+    call dbf_write_attribute_03(hshp_outermost, ishape, 4, merge(1, 0, e%twice))
 
     call dbf_write_attribute_03(hshp_max_speed, ishape, 0, &
          e%max_speed_contour%area)

Tests/plot_snapshot.py

@@ -14,7 +14,6 @@ import shapefile
 import numpy as np
 import matplotlib.pyplot as plt
 import netCDF4
-import matplotlib.markers
 
 parser = argparse.ArgumentParser()
 parser.add_argument("-v", "--velocity", help = "plot velocity field",
@@ -30,40 +29,32 @@ lon_2d, lat_2d = np.meshgrid(longitude, latitude)
 plt.figure()
 plt.plot(lon_2d.reshape(-1), lat_2d.reshape(-1), "+", color="gray")
 
-# Extrema:
-
 reader = shapefile.Reader("extremum_1")
-points = []
-cyclone = []
-suff_amp_list = []
-
-for shape_rec in reader.iterShapeRecords():
-    points.append(shape_rec.shape.points[0])
-    cyclone.append(shape_rec.record[4])
-    suff_amp_list.append(shape_rec.record[5])
-
-points = np.array(points)
-colors = np.choose(cyclone, ("red", "blue"))
-
-plt.scatter(points[:, 0], points[:, 1], s = 100, edgecolors = colors,
-            facecolors = "none")
-
-for i, e in enumerate(points, start = 1):
-    plt.annotate(str(i), e, xytext = (3, 3), textcoords = "offset points")
-
-# Outermost and max-speed contours:
-
 reader_outer = shapefile.Reader("outermost_contour_1")
 reader_m_s = shapefile.Reader("max_speed_contour_1")
 
-for shape_outer, shape_m_s, my_color, suff_amp in zip(reader_outer.iterShapes(),
-                                                      reader_m_s.iterShapes(),
-                                                      colors, suff_amp_list):
+for shape_rec_extr, shape_outer, shape_m_s in zip(reader.iterShapeRecords(),
+                                                  reader_outer.iterShapes(),
+                                                  reader_m_s.iterShapes()):
+    points = shape_rec_extr.shape.points[0]
+    
+    if shape_rec_extr.record[4] == 0:
+            # Anti-cyclone
+        color = "red"
+    else:
+        color = "blue"
+        
+    plt.plot(points[0], points[1], "o", markersize = 10, color = color,
+             fillstyle = "none")
+    plt.annotate(str(shape_rec_extr.record[2]), points, xytext = (3, 3),
+                 textcoords = "offset points")
+
     if shape_outer.shapeType != shapefile.NULL:
         points = np.array(shape_outer.points)
-        lines = plt.plot(points[:, 0], points[:, 1], color = my_color)
+        lines = plt.plot(points[:, 0], points[:, 1], color = color)
 
-        if suff_amp == 0:
+        if shape_rec_extr.record[5] == 0:
+            # Insufficient amplitude
             lines[0].set_marker("s")
             lines[0].set_fillstyle("none")
         elif shape_m_s.shapeType == shapefile.NULL:
@@ -73,7 +64,7 @@ for shape_outer, shape_m_s, my_color, suff_amp in zip(reader_outer.iterShapes(),
 
         if shape_m_s.shapeType != shapefile.NULL:
             points = np.array(shape_m_s.points)
-            plt.plot(points[:, 0], points[:, 1], color = my_color, marker = "o")
+            plt.plot(points[:, 0], points[:, 1], color = color, 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
+
+reader = shapefile.Reader("extremum_1")
+reader_outer = shapefile.Reader("outermost_contour_1")
+reader_m_s = shapefile.Reader("max_speed_contour_1")
+
+plt.figure()
+
+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:
+        # Sufficient amplitude
+        points = shape_rec_extr.shape.points[0]
+
+        if shape_rec_extr.record[4] == 0:
+            # Anti-cyclone
+            color = "red"
+        else:
+            color = "blue"
+        
+        plt.plot(points[0], points[1], "+", markersize = 10, color = color,
+                 fillstyle = "none")
+
+        points = np.array(shape_outer.points)
+        lines = plt.plot(points[:, 0], points[:, 1], color = color)
+
+        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)
+            plt.plot(points[:, 0], points[:, 1], color = color)
+
+plt.xlabel("longitude (degrees east)")
+plt.ylabel("latitude (degrees north)")
+
+plt.show()

Tests/test_get_snapshot.f

@@ -24,7 +24,7 @@ program test_get_snapshot
 
   TYPE(shpfileobject) hshp_outermost
   ! shapefile outermost_contour_$m. The fields in the DBF file are,
-  ! in that order: area, ssh, date index, eddy index.
+  ! in that order: area, ssh, date index, eddy index, twice.
 
   TYPE(shpfileobject) hshp_max_speed
   ! shapefile x_speed_contour_$m. The fields in the DBF file are, in
@@ -85,6 +85,8 @@ program test_get_snapshot
        nwidth = 4, ndecimals = 0)
   call dbf_add_field_03(ifield, hshp_outermost, 'eddy_index', ftinteger, &
        nwidth = 5, ndecimals = 0)
+  call dbf_add_field_03(ifield, hshp_outermost, 'twice', ftinteger, &
+       nwidth = 1, ndecimals = 0)
 
   call shp_create_03("max_speed_contour_1", shpt_polygon, hshp_max_speed)
   call dbf_add_field_03(ifield, hshp_max_speed, 'area', ftdouble, nwidth = 20, &