Skip to content
Snippets Groups Projects
Commit a90ebb4c authored by Lionel GUEZ's avatar Lionel GUEZ
Browse files

Remove script

No longer useful since there are no interpolated eddies any longer.
parent a06cbf67
No related branches found
No related tags found
No related merge requests found
#!/usr/bin/env python3
"""This script renumbers interpolated eddies. It is useful only to
compare the results of runs of eddy_graph with different numbers of
MPI processes, for debugging. Renumbering is done after sorting by
date index, then latitude then longitude. Even for a run with one MPI
process, the interpolated eddies as output by the run may not be
sorted this way.
"""
import numpy as np
import argparse
import fiona
from os import path
import csv
import os
import sys
import shutil
parser = argparse.ArgumentParser()
parser.add_argument("input_dir",
help = "containing the SHPC of interpolated eddies and "
"edge list")
args = parser.parse_args()
if path.samefile(args.input_dir, os.getcwd()):
sys.exit("The input directory cannot be the current directory.")
my_input = path.join(args.input_dir, "number_eddies.csv")
n_vis, n_interp = np.loadtxt(my_input, skiprows=2, dtype=int, usecols=(1, 2),
unpack = True)
with open(my_input) as f:
f.readline()
f.readline()
k1 = f.readline().split()[0]
k1 = int(k1)
lon = []
lat = []
date_index = []
eddy_index_old = []
my_input = path.join(args.input_dir, "SHPC")
with fiona.open(my_input, layer = "extremum") as extremum:
for record in extremum:
lon.append(record["geometry"]["coordinates"][0])
lat.append(record["geometry"]["coordinates"][1])
date_index.append(record["properties"]["date_index"])
eddy_index_old.append(record["properties"]["eddy_index"])
lon = np.array(lon)
lat = np.array(lat)
date_index = np.array(date_index)
eddy_index_old = np.array(eddy_index_old)
ind = np.lexsort((lon, lat, date_index))
previous_date = None
eddy_index_map = {}
for ishape in ind:
if date_index[ishape] == previous_date:
eddy_index_new += 1
else:
eddy_index_new = int(n_vis[date_index[ishape] - k1]) + 1
previous_date = date_index[ishape]
if eddy_index_new != eddy_index_old[ishape]:
eddy_index_map[(date_index[ishape], eddy_index_old[ishape])] \
= eddy_index_new
for layer in fiona.listlayers(my_input):
with fiona.open(my_input, layer = layer) as src:
with fiona.open("SHPC", "w", layer = layer, **src.meta) as sink:
for ishape in ind:
record = src[int(ishape)]
try:
record["properties"]["eddy_index"] \
= eddy_index_map[(date_index[ishape],
eddy_index_old[ishape])]
except KeyError:
pass
sink.write(record)
ishape_last = np.cumsum(n_interp) - 1
i = np.searchsorted(ishape_last, 0)
# (Ensure that the first line corresponds to the first date present in
# interpolated eddies.)
np.savetxt("SHPC/ishape_last.txt", ishape_last[i:], fmt = "%d")
my_input = path.join(args.input_dir, "SHPC/grid_nml.txt")
shutil.copy(my_input, "SHPC")
for orientation in ["cyclo", "anti"]:
my_input = path.join(args.input_dir, f"edgelist_{orientation}.csv")
edgelist = np.genfromtxt(my_input, skip_header = 1, names = True,
dtype = "i8," * 4 + "f8")
for row in edgelist:
for k, i in [("k1", "i1"), ("k2", "i2")]:
if row[i] > n_vis[row[k] - k1]:
try:
row[i] = eddy_index_map[(row[k], row[i])]
except KeyError:
pass
edgelist.sort(order = ["k1", "i1"])
with open(f"edgelist_{orientation}.csv", "w", newline='') as sink:
# Title lines:
sink.write('"predecessor date subscript" "predecessor eddy subscript" '
'"successor date subscript" "successor eddy subscript"\n')
sink.write("k1 i1 k2 i2 weight\n")
writer = csv.writer(sink, delimiter = " ", lineterminator = "\n")
writer.writerows(edgelist)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment