inst_eddies.py.in

#!/usr/bin/env python3

"""This script creates the directories SHPC/Slice_{my_slice}/Cyclones
and SHPC/Slice_{my_slice}/Anticyclones, containing a time slice of
instantaneous eddies, and the file perf_report.csv. If
SHPC/Slice_{my_slice}/{orientation} or perf_report.csv already exist
then they are replaced.

"""

import shutil
import nco
import os
import subprocess
import shapefile
import sys
import time
import csv

def loop_inst_eddies(files, bbox, d, my_slice):
    """files is a sequence of file paths. bbox is a tuple of four real
    numbers. d is a date index (integer value) for the first
    date. my_slice is the slice index (integer value) in the SHPC.

    """

    inst_eddies_exe = "@CMAKE_CURRENT_BINARY_DIR@/inst_eddies"

    if not os.access(inst_eddies_exe, os.X_OK):
        sys.exit(inst_eddies_exe + " not found or not executable")

    if os.access("config_nml.txt", os.R_OK):
        print("inst_eddies.loop_inst_eddies: Will use config_nml.txt.")
    else:
        sys.exit('"config_nml.txt" not found')

    nco_instance = nco.Nco()

    for orient in ["Cyclones", "Anticyclones"]:
        my_dir = f"SHPC/Slice_{my_slice}/{orient}"
        if os.access(my_dir, os.F_OK): shutil.rmtree(my_dir)
        os.makedirs(my_dir)

    perf_report = open("perf_report.csv", "w", newline = '')
    writer = csv.writer(perf_report, lineterminator = "\n")
    writer.writerow(["", "s", "s"]) # units
    writer.writerow(["date", "elapsed NCO", "elapsed Fortran"]) # long names
    writer.writerow(["date", "elapsed_NCO", "elapsed_Fortran"])

    for nc_file in files:
        print("inst_eddies.loop_inst_eddies: Date:", d)
        print("inst_eddies.loop_inst_eddies: Input file:", nc_file,
              flush = True)

        if os.access(nc_file, os.F_OK):
            my_pc = time.perf_counter()

            if bbox:
                xmin, xmax, ymin, ymax = 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"])
            os.symlink("unpacked.nc", "h.nc")
            os.symlink("unpacked.nc", "uv.nc")
            elapsed_NCO = time.perf_counter() - my_pc
            my_pc = time.perf_counter()
            subprocess.run(inst_eddies_exe, check = True, input = "&dates_nml "
                           f"date = {d}, slice = {my_slice}/\n", text = True)
            elapsed_Fortran = time.perf_counter() - my_pc
            os.remove("h.nc")
            os.remove("uv.nc")
            writer.writerow([d, elapsed_NCO, elapsed_Fortran])
        else:
            print("inst_eddies.loop_inst_eddies: Missing file:", nc_file)

            for orient in ["Cyclones", "Anticyclones"]:
                with open(f"SHPC/Slice_{my_slice}/{orient}/ishape_last.txt",
                          "a") as ishape_last, \
                    shapefile.Reader(f"SHPC/Slice_{my_slice}/{orient}/"
                                     "extremum") as shpr:
                    n_shapes = len(shpr)
                    ishape_last.write(str(n_shapes - 1) + "\n")

        d += 1

    perf_report.close()
    os.remove("unpacked.nc")
    if bbox: os.remove("cropped.nc")

    if my_slice == 0 and not os.access("SHPC/n_slices.txt", os.F_OK):
        with open("SHPC/n_slices.txt", "w") as f: f.write("1\n")

if __name__ == "__main__":
    import argparse
    
    argparser = argparse.ArgumentParser()
    argparser.add_argument("file", nargs = "+",
                           help = "NetCDF file containing SSH and velocity "
                           "at a single date")
    argparser.add_argument("-b", "--bbox", nargs=4, type = float,
                           metavar = ("xmin", "xmax", "ymin", "ymax"))
    args = argparser.parse_args()
    d = input("Enter first date index (integer value): ")
    d = int(d)
    my_slice = input("Enter slice index (integer value): ")
    my_slice = int(my_slice)
    loop_inst_eddies(args.file, args.bbox, d, my_slice)