-
Lionel GUEZ authored
The vertex property `avg_Rossby` was never set after `n_valid_Rossby == 0`. It appears that it had then a default value of 0. Better to explicitly set it to 0.
Lionel GUEZ authoredThe vertex property `avg_Rossby` was never set after `n_valid_Rossby == 0`. It appears that it had then a default value of 0. Better to explicitly set it to 0.
cost_function.py 8.87 KiB
#!/usr/bin/env python3
"""A script that takes the graph of segments without cost functions
and computes the non-local cost functions applied to edges.
Input:
- "node_id_param.json", expected to be in the current directory
- the graph of segments without cost functions, "segments.gt" or
"segments.graphml", expected to be in the current directory
- shapefiles, specified as command line arguments
Output: the graph of segments with cost functions,
'segments_cost_functions.gt'
"""
import graph_tool
import time
import json
import math
import report_graph
import util_eddies
import bisect
import argparse
Omega = 2 * math.pi / 86164.
def calculate_radii_rossby(list_eddies, e_overestim, handlers, array_d_init):
"""Compute average on list_eddies of Rossby number and radius of
maximum speed contour.
"""
avg_rad = 0 # in m
avg_Rossby = 0
n_valid_Rossby = 0
for n in list_eddies:
date_index, eddy_index = report_graph.node_to_date_eddy(n, e_overestim)
i_SHPC = get_SHPC(array_d_init, date_index)
ishape = util_eddies.comp_ishape(handlers[i_SHPC], date_index,
eddy_index)
# Now that we have the location in the shapefiles, we need to
# get the radius and the rossby number:
shapeRec = handlers[i_SHPC]["readers"]["extremum"].shapeRecord(ishape)
lat_in_deg = shapeRec.shape.points[0][1] # in degrees
f = 2 * Omega * math.sin(math.radians(lat_in_deg)) # in s-1
V_max = shapeRec.record[4] # in m/s
radius = handlers[i_SHPC]["readers"]["max_speed_contour"]\
.record(ishape)['r_eq_area'] * 1000 # in m
if abs(V_max) < 100:
avg_Rossby += V_max / (f * radius)
n_valid_Rossby += 1
avg_rad += radius # in m
avg_rad /= len(list_eddies)
if n_valid_Rossby != 0: avg_Rossby /= n_valid_Rossby
return {"avg_rad": avg_rad, "avg_Rossby": avg_Rossby}
def get_SHPC(array_d_ini, date_index):
i_SHPC = bisect.bisect(array_d_init, date_index) - 1
assert i_SHPC >= 0
return i_SHPC
t0 = time.perf_counter()
timings = open("timings.txt", "w")
parser = argparse.ArgumentParser()
parser.add_argument("SHPC_dir", nargs='+')
parser.add_argument("--graphml", action = "store_true",
help = "save to graphml format")
args = parser.parse_args()
# Grab e_overestim:
with open("node_id_param.json") as f: node_id_param = json.load(f)
e_overestim = node_id_param["e_overestim"]
# Set some values needed for the cost function:
delta_cent_mean = 3.8481 # [in km]
delta_cent_std = 8.0388
delta_ro_mean = -0.0025965
delta_ro_std = 5.2168
delta_r_mean = -0.0094709 * 1000 # [in m]
delta_r_std = 8.6953 * 1000
# Load the graph_tool file:
print('Loading graph...')
try:
g = graph_tool.load_graph('segments.gt')
except FileNotFoundError:
g = graph_tool.load_graph('segments.graphml')
print('Loading done...')
print("Input graph:")
print("Number of vertices:", g.num_vertices())
print("Number of edges:", g.num_edges())
print("Internal properties:")
g.list_properties()
t1 = time.perf_counter()
timings.write(f"loading: {t1 - t0:.0f} s\n")
t0 = t1
g.vp['pos_first'] = g.new_vp('object')
g.vp['pos_last'] = g.new_vp('object')
g.vp['first_av_rad'] = g.new_vp('float')
g.vp['first_av_ros'] = g.new_vp('float')
g.vp['last_av_rad'] = g.new_vp('float')
g.vp['last_av_ros'] = g.new_vp('float')
g.ep['cost_function'] = g.new_ep('float')
# Set up the list of SHPC:
handlers = [util_eddies.open_shpc(shpc_dir) for shpc_dir in args.SHPC_dir]
array_d_init = [handler["d_init"] for handler in handlers]
# (create the list once and for all)
num_of_days_to_avg = 7 # number of days to average
print("Iterating on vertices...")
for n in g.vertices():
if n.in_degree() != 0 or n.out_degree() != 0:
segment = g.vp.inst_eddies[n]
num_of_days = len(segment)
# Calculate the date index, the eddy index and the SHPC index of
# the first instantaneous eddy in the segment, then grab the
# position of the extremum and store it in a vertex property:
date_index, eddy_index = report_graph.node_to_date_eddy(segment[0],
e_overestim)
i_SHPC = get_SHPC(array_d_init, date_index)
ishape = util_eddies.comp_ishape(handlers[i_SHPC], date_index,
eddy_index)
g.vp.pos_first[n] = handlers[i_SHPC]["readers"]["extremum"]\
.shape(ishape).points[0] # [in degrees]
# Same for last instantaneous eddy in the segment:
date_index, eddy_index = report_graph.node_to_date_eddy(segment[-1],
e_overestim)
i_SHPC = get_SHPC(array_d_init, date_index)
ishape = util_eddies.comp_ishape(handlers[i_SHPC], date_index,
eddy_index)
g.vp.pos_last[n] = handlers[i_SHPC]["readers"]["extremum"]\
.shape(ishape).points[0] # [in degrees]
if (num_of_days > num_of_days_to_avg):
# The segment is longer than the number of days over which
# to average
# First 7 days calculation
first_res = calculate_radii_rossby(segment[:num_of_days_to_avg],
e_overestim, handlers,
array_d_init)
# Average and assign the first radii:
g.vp.first_av_rad[n] = first_res['avg_rad']
# Average and assign the rossbies:
g.vp.first_av_ros[n] = first_res['avg_Rossby']
# Last 7 days calculation:
last_res = calculate_radii_rossby(segment[- num_of_days_to_avg:],
e_overestim, handlers,
array_d_init)
# Average and assign the last radii
g.vp.last_av_rad[n] = last_res['avg_rad']
# Average and assign the rossbies:
g.vp.last_av_ros[n] = last_res['avg_Rossby']
else:
# The number of eddies in a segment is lower than the number
# of days over which to average. The values will be the same
# except for the positions.
res = calculate_radii_rossby(segment, e_overestim, handlers,
array_d_init)
# Average and assign the rossbies:
avg_Rossby = res['avg_Rossby']
g.vp.first_av_ros[n] = avg_Rossby
g.vp.last_av_ros[n] = avg_Rossby
# Average and assign the radii
avg_rad = res['avg_rad']
g.vp.first_av_rad[n] = avg_rad
g.vp.last_av_rad[n] = avg_rad
t1 = time.perf_counter()
timings.write(f"iterating on vertices: {t1 - t0:.0f} s\n")
t0 = t1
print("Iterating on edges...")
for edge in g.edges():
source_node = edge.source()
target_node = edge.target()
cf = -10000
lat_for_conv = (g.vp.pos_last[source_node][1] +
g.vp.pos_first[target_node][1]) / 2
# (latitude needed for conversion of degrees to kilometers)
lat_for_conv = math.radians(lat_for_conv) # need to convert to radians
# because of the wrapping issue (360° wrapping incorrectly to 0°),
# we check for that here
lon_diff = abs(g.vp.pos_last[source_node][0] \
- g.vp.pos_first[target_node][0])
if (lon_diff > 300):
lon_diff = 360 - lon_diff
Delta_Cent = math.sqrt((lon_diff * 111.32 * math.cos(lat_for_conv))**2
+ ((g.vp.pos_last[source_node][1]
- g.vp.pos_first[target_node][1]) * 110.574)**2)
# calculate the first term
first_term = ((Delta_Cent - delta_cent_mean)/delta_cent_std) ** 2
# Rossbies:
if (g.vp.first_av_ros[target_node] and g.vp.last_av_ros[source_node]):
Delta_Ro = g.vp.last_av_ros[source_node] \
- g.vp.first_av_ros[target_node]
else:
# At least one of the rossbies is invalid.
# Delta_Ro = delta_ro_mean
Delta_Ro = 0
# Calculate the second term
second_term = ((Delta_Ro - delta_ro_mean)/delta_ro_std ) ** 2
# R_Vmax 1 and 2 already exist, just get the delta
Delta_R_Vmax = g.vp.last_av_rad[source_node] \
- g.vp.first_av_rad[target_node]
# Calculate the third term
third_term = ((Delta_R_Vmax - delta_r_mean)/delta_r_std) ** 2
#############################
# calculate the cost function
#############################
cf = math.sqrt(first_term + second_term + third_term)
# assign as weight to the edge
g.ep.cost_function[edge] = cf
t1 = time.perf_counter()
timings.write(f"iterating on edges: {t1 - t0:.0f} s\n")
t0 = t1
print("Saving...")
if args.graphml:
g.save('segments_cost_functions.graphml')
else:
g.save('segments_cost_functions.gt')
print('All done')
t1 = time.perf_counter()
timings.write(f"saving: {t1 - t0:.0f} s\n")
timings.close()