From 93dba70e116cb2b5dd1c60968ec82f04cfe1bb1b Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Sat, 13 Mar 2021 00:36:58 +0100
Subject: [PATCH] New script to loop on `inst_eddies`
New script to loop on detection of instantaneous eddies.
---
Inst_eddies/inst_eddies.py | 83 ++++++++++++++++++++++++++++++++++++++
1 file changed, 83 insertions(+)
create mode 100755 Inst_eddies/inst_eddies.py
diff --git a/Inst_eddies/inst_eddies.py b/Inst_eddies/inst_eddies.py
new file mode 100755
index 00000000..d8d3090a
--- /dev/null
+++ b/Inst_eddies/inst_eddies.py
@@ -0,0 +1,83 @@
+#!/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")
--
GitLab