New script to loop on `inst_eddies`

New script to loop on detection of instantaneous eddies.

Inst_eddies/inst_eddies.py

+#!/usr/bin/env python3
+
+import argparse
+import datetime
+import shutil
+import nco
+import os
+from os import path
+import subprocess
+import shapefile
+
+parser = argparse.ArgumentParser()
+parser.add_argument("first_file")
+parser.add_argument("-l", "--last_date", help = "%%Y-%%m-%%d")
+parser.add_argument("-b", "--bbox", nargs=4, type = float,
+                    metavar = ("xmin", "xmax", "ymin", "ymax"))
+args = parser.parse_args()
+
+if os.access("SHPC_all_dates", os.F_OK): shutil.rmtree("SHPC_all_dates")
+dirname, basename = path.split(args.first_file)
+my_date = datetime.datetime.strptime(basename,
+                                     "dt_global_allsat_phy_l4_%Y%m%d"
+                                     "_20190101.nc").date()
+nc_file = args.first_file
+nco_instance = nco.Nco()
+
+if args.last_date:
+    last_date = datetime.datetime.strptime(args.last_date, "%Y-%m-%d").date()
+else:
+    last_date = my_date
+
+while True:
+    if os.access(nc_file, os.F_OK):
+        if args.bbox:
+            xmin, xmax, ymin, ymax = args.bbox
+            nco_instance.ncks(nc_file, output = "cropped.nc",
+                              options = [f"--dimension=longitude,{xmin},{xmax}",
+                                         f"--dimension=latitude,{ymin},{ymax}"])
+            nc_file = "cropped.nc"
+            
+        nco_instance.ncpdq(nc_file, output = "unpacked.nc",
+                           options = ["--unpack"])
+        nc_file = "unpacked.nc"
+        os.symlink(nc_file, "h.nc")
+        os.symlink(nc_file, "uv.nc")
+        os.mkdir("SHP_triplet")
+        
+        with open("main_nml.txt") as nml_file:
+            subprocess.run("inst_eddies", stdin = nml_file, check = True)
+
+        os.remove("h.nc")
+        os.remove("uv.nc")
+    else:
+        print("Missing file:", nc_file)
+
+    if os.access("SHPC_all_dates", os.F_OK):
+        if os.access("SHP_triplet", os.F_OK):
+            for filename in ["extremum", "max_speed_contour",
+                             "outermost_contour"]:
+                from_file = path.join("SHP_triplet", filename)
+                to_file = path.join("SHPC_all_dates", filename)
+
+                for command in ["dbfcat", "shpcat"]:
+                    subprocess.run([command, from_file, to_file], check = True)
+        
+            shutil.rmtree("SHP_triplet")
+            
+        with open("SHPC_all_dates/ishape_last.txt", "a") as \
+             ishape_last, shapefile.Reader("SHPC_all_dates/extremum") \
+             as shpr:
+            n_shapes = len(shpr)
+            ishape_last.write(str(n_shapes - 1) + "\n")
+    else:
+        if os.access("SHP_triplet", os.F_OK):
+            os.rename("SHP_triplet", "SHPC_all_dates")
+
+    my_date += datetime.timedelta(1)
+    if my_date > last_date: break
+    basename = my_date.strftime("dt_global_allsat_phy_l4_%Y%m%d_20190101.nc")
+    nc_file = path.join(dirname, basename)
+
+os.remove("unpacked.nc")
+os.remove("cropped.nc")