Write two SHPT for cyclones and anticyclones

To progressively replace the output of Matlab, we modify it by
steps. So, for now at least, we want to keep the Matlab indexing of
eddies, which is separate for cyclones and anticyclones. So we have to
create two triplets of shapefiles, one for cyclones and one for anticyclones.

Analysis/mat_to_shapefiles.py

@@ -8,15 +8,19 @@
 
 import shapefile
 import numpy as np
+from os import path
 
-def write(matlab_data):
-    """matlab_data is a dictionary."""
+def write(matlab_data, eddy_orientation, SHP_triplet):
+    """matlab_data is a dictionary. matlab_data[eddy_orientation] is a
+    numpy.ndarray with shape (number of eddies, 21).
 
-    with shapefile.Writer("Snapshot/extremum", 
+    """
+
+    with shapefile.Writer(path.join(SHP_triplet, "extremum"),
                           shapeType = shapefile.POINT) as w_extr, \
-      shapefile.Writer("Snapshot/outermost_contour", 
+      shapefile.Writer(path.join(SHP_triplet, "outermost_contour"),
                        shapeType = shapefile.POLYGON) as w_outer, \
-      shapefile.Writer("Snapshot/max_speed_contour", 
+      shapefile.Writer(path.join(SHP_triplet, "max_speed_contour"),
                        shapeType = shapefile.POLYGON) as w_max_speed:
         w_extr.field("ssh", "N", 13, 6)
         w_extr.field("date_index", "N", 6)
@@ -39,49 +43,46 @@ def write(matlab_data):
 
         eddy_index = 0
 
-        for eddy_orientation in ["Anticyclonic_Cell", "Cyclonic_Cell"]:
-            # matlab_data[eddy_orientation] is a numpy.ndarray with
-            # shape (number of eddies, 21).
-            cyclone = eddy_orientation == "Cyclonic_Cell"
-
-            for eddy in matlab_data[eddy_orientation]:
-                # eddy is a numpy.ndarray with shape (21,). The
-                # description of each element of eddy is given by
-                # matlab_data["Fields"]. Note that amplitudes in
-                # eddy[7] and eddy[11] are positive regardless of eddy
-                # orientation.
-                eddy_index += 1
-
-                w_extr.point(eddy[0], eddy[1])
-                i = np.argwhere(matlab_data["X"] == eddy[0]).item()
-                j = np.argwhere(matlab_data["Y"] == eddy[1]).item()
-                speed = 1e4 if np.isnan(eddy[12]) else eddy[12]
-                w_extr.record(ssh = matlab_data["ADT"][i, j], 
-                              date_index = matlab_data["date_num"],
-                              eddy_index = eddy_index, interpolat = 0, 
-                              cyclone = cyclone, valid = 1, speed = speed)
-
-                if cyclone:
-                    ssh = matlab_data["ADT"][i, j] + eddy[7]
-                else:
-                    ssh =  matlab_data["ADT"][i, j] - eddy[7]
-
-                w_outer.record(r_eq_area = eddy[6], ssh = ssh, 
+        cyclone = eddy_orientation == "Cyclonic_Cell"
+
+        for eddy in matlab_data[eddy_orientation]:
+            # eddy is a numpy.ndarray with shape (21,). The
+            # description of each element of eddy is given by
+            # matlab_data["Fields"]. Note that amplitudes in
+            # eddy[7] and eddy[11] are positive regardless of eddy
+            # orientation.
+            eddy_index += 1
+
+            w_extr.point(eddy[0], eddy[1])
+            i = np.argwhere(matlab_data["X"] == eddy[0]).item()
+            j = np.argwhere(matlab_data["Y"] == eddy[1]).item()
+            speed = 1e4 if np.isnan(eddy[12]) else eddy[12]
+            w_extr.record(ssh = matlab_data["ADT"][i, j],
+                          date_index = matlab_data["date_num"],
+                          eddy_index = eddy_index, interpolat = 0,
+                          cyclone = cyclone, valid = 1, speed = speed)
+
+            if cyclone:
+                ssh = matlab_data["ADT"][i, j] + eddy[7]
+            else:
+                ssh =  matlab_data["ADT"][i, j] - eddy[7]
+
+            w_outer.record(r_eq_area = eddy[6], ssh = ssh,
+                           date_index = matlab_data["date_num"],
+                           eddy_index = eddy_index, radius4 = - 1)
+            polyline = np.stack((eddy[4], eddy[5]), axis = 1)
+            w_outer.poly([polyline])
+
+            if cyclone:
+                ssh = matlab_data["ADT"][i, j] + eddy[11]
+            else:
+                ssh =  matlab_data["ADT"][i, j] - eddy[11]
+
+            w_max_speed.record(r_eq_area = eddy[10], ssh = ssh,
                                date_index = matlab_data["date_num"], 
-                               eddy_index = eddy_index, radius4 = - 1)
-                polyline = np.stack((eddy[4], eddy[5]), axis = 1)
-                w_outer.poly([polyline])
-
-                if cyclone:
-                    ssh = matlab_data["ADT"][i, j] + eddy[11]
-                else:
-                    ssh =  matlab_data["ADT"][i, j] - eddy[11]
-
-                w_max_speed.record(r_eq_area = eddy[10], ssh = ssh, 
-                                   date_index = matlab_data["date_num"], 
-                                   eddy_index = eddy_index)
-                polyline = np.stack((eddy[8], eddy[9]), axis = 1)
-                w_max_speed.poly([polyline])
+                               eddy_index = eddy_index)
+            polyline = np.stack((eddy[8], eddy[9]), axis = 1)
+            w_max_speed.poly([polyline])
 
 if __name__ == "__main__":
     import sys
@@ -89,4 +90,5 @@ if __name__ == "__main__":
     
     if len(sys.argv) != 2: sys.exit("Required argument: .mat file")
     matlab_data = sio.loadmat(sys.argv[1], squeeze_me = True)
-    write(matlab_data)
+    write(matlab_data, "Anticyclonic_Cell", "Snapshot_anti")
+    write(matlab_data, "Cyclonic_Cell", "Snapshot_cyclo")