Add scripts to compare with Matlab output

Add scripts to compare with Matlab output and script to concatenate shapefiles.

Analysis/convert_Matlab_v6.sh

+if (($# != 1))
+then
+    echo "Required argument: input \".mat\" file"
+    exit 1
+fi
+
+set -xe
+
+my_base=`basename $1 .mat`
+
+cat >/tmp/convert_v6_$$.m <<EOF
+clear
+cd $PWD
+load('$1')
+save('${my_base}_v6', '-v6')
+exit
+EOF
+
+matlab -nojvm -r "run('/tmp/convert_v6_$$.m')"

Analysis/make_graph_from_Matlab.py

+#!/usr/bin/env python3
+
+import scipy.io as sio
+
+Assoc_eddies_Anti_max = sio.loadmat("Association_eddies_Anti_max_v6.mat",
+                                          squeeze_me=True)
+Assoc_eddies_Cyclo_max = sio.loadmat("Association_eddies_Cyclo_max_v6.mat",
+                                          squeeze_me=True)
+Nanti = sio.loadmat("Parameters_Anticyclonic_Eddies_v6.mat", 
+                    squeeze_me = True, variable_names = ["Nanti"])["Nanti"]
+Ncyclo = sio.loadmat("Parameters_Cyclonic_Eddies_v6.mat", squeeze_me = True,
+                     variable_names = ["Ncyclo"])["Ncyclo"]
+
+def convert_FEI(orientation, eddy_index_Matlab):
+    """Convert from Matlab eddy identification to eddy
+    identification in shapefiles.
+
+    We assume that the shapefiles are created by the script
+    "mat_to_shapefiles.py". This script processes a single snapshot
+    and all the anticyclones are processed before all the cyclones."""
+
+    if orientation == "anticyclone":
+        if eddy_index_Matlab <= Nanti[0]:
+            date = Assoc_eddies_Anti_max["date_num"][0]
+            eddy_index_F = eddy_index_Matlab
+        else:
+            date = Assoc_eddies_Anti_max["date_num"][1]
+            eddy_index_F = eddy_index_Matlab - 100
+    else:
+        if eddy_index_Matlab <= Ncyclo[0]:
+            date = Assoc_eddies_Anti_max["date_num"][0]
+            eddy_index_F = Nanti[0] + eddy_index_Matlab
+        else:
+            date = Assoc_eddies_Anti_max["date_num"][1]
+            eddy_index_F = Nanti[1] + eddy_index_Matlab - 100
+
+    return date, eddy_index_F
+
+def print_edgelist(orientation, total_number, id_child):
+    for i1 in range(total_number):
+        predecessor = (orientation, i1 + 1)
+        try:
+            for eddy_index_Matlab in id_child[i1]:
+               print(*convert_FEI(*predecessor), 
+                     *convert_FEI(orientation, eddy_index_Matlab))
+        except TypeError:
+            print(*convert_FEI(*predecessor), 
+                   *convert_FEI(orientation, id_child[i1]))
+
+print_edgelist("anticyclone", Nanti[0], Assoc_eddies_Anti_max["id_child"][:, 0])
+print_edgelist("cyclone", Ncyclo[0], Assoc_eddies_Cyclo_max["id_child"][:, 0])

Analysis/mat_to_shapefiles.py

+#!/usr/bin/env python3
+
+import shapefile
+import numpy as np
+
+def write(matlab_data):
+    """matlab_data is a dict."""
+
+    with shapefile.Writer("Snapshot/extremum", 
+                          shapeType = shapefile.POINT) as w_extr, \
+      shapefile.Writer("Snapshot/outermost_contour", 
+                       shapeType = shapefile.POLYGON) as w_outer, \
+      shapefile.Writer("Snapshot/max_speed_contour", 
+                       shapeType = shapefile.POLYGON) as w_max_speed:
+        w_extr.field("ssh", "N", 13, 6)
+        w_extr.field("date_index", "N", 6)
+        w_extr.field("eddy_index", "N", 5)
+        w_extr.field("interpolat", "N", 1)
+        w_extr.field("cyclone", "N", 1)
+        w_extr.field("valid", "N", 1)
+        w_extr.field("speed", "N", 13, 6)
+
+        w_outer.field("r_eq_area", "N", 10, 4)
+        w_outer.field("ssh", "N", 13, 6)
+        w_outer.field("date_index", "N", 6)
+        w_outer.field("eddy_index", "N", 5)
+        w_outer.field("radius4", "N", 2)
+
+        w_max_speed.field("r_eq_area", "N", 10, 4)
+        w_max_speed.field("ssh", "N", 13, 6)
+        w_max_speed.field("date_index", "N", 6)
+        w_max_speed.field("eddy_index", "N", 5)
+
+        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).
+
+            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
+                cyclone = eddy_orientation == "Cyclonic_Cell"
+
+                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)
+                polyline = np.stack((eddy[8], eddy[9]), axis = 1)
+                w_max_speed.poly([polyline])
+
+if __name__ == "__main__":
+    import sys
+    import scipy.io as sio
+    
+    if len(sys.argv) != 2: sys.exit("Required argument: .mat file")
+    matlab_data = sio.loadmat(sys.argv[1], squeeze_me = True)
+    write(matlab_data)

loop_shpcat.py

+#!/usr/bin/env python3
+
+import datetime
+import subprocess
+from os import path
+import os
+import time
+import glob
+import shutil
+
+initial_dir = os.getcwd()
+
+for year in range(2011, 2017):
+    my_path = path.join(initial_dir, str(year))
+    os.chdir(my_path)
+
+    for src in glob.iglob(f"{year}_01/{year}_01_01/*"):
+        shutil.copy(src, os.curdir)
+
+    my_date = datetime.date(year, 1, 2)
+    final_date = datetime.date(year, 12, 31)
+    t0 = time.perf_counter()
+
+    with  open("perf_report.csv", "w") as perf_report:
+        perf_report.write("elapsed time, in s\n")
+
+        while my_date <= final_date:
+            my_dir = my_date.strftime("%Y_%m/%Y_%m_%d")
+
+            for filename in ["extremum", "max_speed_contour", 
+                             "outermost_contour"]:
+                from_file = path.join(my_dir, filename)
+
+                for command in ["dbfcat", "shpcat"]:
+                    subprocess.run([command, from_file, filename], check = True)
+
+            t1 = time.perf_counter()
+            perf_report.write(str(t1 - t0) + "\n")
+            t0 = t1
+            my_date += datetime.timedelta(1)