Import script written by Mevludin ISIC

gt_segmented_cf.py

+#!/usr/bin/env python3
+
+# A script that takes a segmented graph in the gt format and performs the 
+# non-local cost function calculation where each each will have a cost 
+# function assigned to it. 
+# Inputs (to be changed within the script):
+    # Orientation
+    # node_id_param_file location and name
+    # gt graph name and location
+    # shapefiles names and location
+# Output: a gt graph with cost functions and an xml graph with cost functions
+
+import graph_tool
+import time
+import json
+import math
+from os import path
+import shapefile
+import datetime
+from numpy import loadtxt
+
+######################
+# Utility functions: #
+######################
+
+def node_to_date_eddy(node_index, e_overestim, d_init):
+    """Convert from node identification in graph to eddy identification in
+    shapefiles.
+
+    """
+    k = (node_index - 1) // e_overestim
+    d = d_init + k
+    return {"eddy_index": node_index - k * e_overestim, "date_index": k,
+            "days_1950": d}
+
+def comp_ishape(d_init, ishape_last, days_1950, eddy_index): 
+    k = days_1950 - d_init
+    
+    if k == 0:
+        ishape = eddy_index - 1
+    else:
+        ishape = ishape_last[k - 1] + eddy_index
+        
+    return ishape 
+
+#######################
+# Necessary functions #
+#######################
+
+def calculate_radii_and_rossby(start, end, inc, segment, e_overestim, d_init, ext_shp_rcd, max_rcd, ishape_last):
+    
+    radii = 0 #[m]
+    rossby = 0 #[1/s]
+    
+    days_modifier = 0
+    
+    for i in range(start, end, inc):
+        current_eddy = node_to_date_eddy(segment[i], e_overestim, d_init)
+        date = datetime.date(1950,1,1) + datetime.timedelta(current_eddy['days_1950'])
+        year = date.year
+        
+        # calculate the location in the shapefile
+        d_init_1 = ext_shp_rcd[year][0].record[1]
+        location = comp_ishape(d_init_1, ishape_last[year], 
+                               current_eddy['days_1950'], current_eddy['eddy_index'])
+        
+        # now that we have the location in the shapefiles, we need to get the radius and the rossby number
+        
+        # rossby:
+        lat_in_deg = ext_shp_rcd[year][location].shape.points[0][1] #[deg]
+        f = 2*2*math.pi/(24*3600)*math.sin(math.radians(lat_in_deg)) # [1/s]
+        
+        V_max = ext_shp_rcd[year][location].record[4] #[m/s]
+        R_Vmax = max_rcd[year][location]['r_eq_area'] * 1000 #[m]
+        
+        
+        if (V_max < 100):
+            # calculate Ro and Delta_Ro
+            Ro = V_max / (f * R_Vmax) #[]
+        else:
+            Ro = 0
+            days_modifier += 1
+        
+        ####### RADII #######
+        radii += R_Vmax # [m]
+        ####### ROSSBY ######
+        rossby += Ro # []
+        
+    return {"radii": radii, "rossby": rossby, "days_modifier": days_modifier}
+
+###############################
+# grab e_overestim and d_init #
+###############################
+
+node_id_param_file = "/data/misic/Convert_Global/node_id_param.json"
+with open(node_id_param_file) as f: node_id_param = json.load(f)
+# assign attributes to the whole graphs
+e_overestim = node_id_param["e_overestim"]
+d_init = node_id_param["d_init"]
+
+################################################
+# set some values needed for the cost function #
+################################################
+
+delta_cent_mean = 3.8481 # [km]
+delta_cent_std = 8.0388
+
+delta_ro_mean = -0.0025965 # []
+delta_ro_std = 5.2168
+
+delta_r_mean = -0.0094709 * 1000 #[m] 
+delta_r_std = 8.6953 * 1000
+
+orientation = 'anti' 
+
+##########################
+# load the edgelist file #
+##########################
+
+t0 = time.perf_counter()
+
+print('Loading gt file...')
+g = graph_tool.Graph()
+g.load(f'global_segmented_{orientation}.gt')
+t1 = time.perf_counter()
+print(f'Loading done: {g}')
+print(t1 - t0)
+
+g.set_fast_edge_removal()
+
+##################################
+# Calculating the Cost Function: #
+##################################
+
+#######################
+# Load the shapefiles #
+#######################
+
+# load the shapefiles and the ishape_last files
+
+print('Loading extrema and ishape_last files.')
+    
+# setup the dicts for the shapefiles and the ishape_last files
+extremum = {}
+max_speed = {}
+ishape_last = {}
+year_bool = {}
+
+shp_tr_dir = '/data/guez/Oceanic_eddies/Global_Matlab'
+    
+for year in range(1993,2019):
+    shp_dir = path.join(shp_tr_dir, f'SHPC_{orientation}_{year}')
+    extremum[year] = shapefile.Reader(path.join(shp_dir, 'extremum.shp'))
+    max_speed[year] = shapefile.Reader(path.join(shp_dir, 'max_speed_contour.shp'))
+    ishape_last[year] = loadtxt(path.join(shp_dir, 'ishape_last.txt'), 
+                                dtype = int, delimiter = "\n", unpack = False) 
+    year_bool[year] = False
+
+print('Loading of the shapefiles and ishape_last files done.')
+
+count = 1 
+
+t0 = time.perf_counter()
+# change if there is a change over the number of days to average
+num_of_days_to_avg = 7
+
+# we iterate on the years because the segments are non sequential in time
+for yr in range(1993, 2019):
+    
+    print(f'Current year: {yr}') 
+
+    # grab the shapefiles
+
+    ext_shp_rcd = {}
+    max_rcd = {}
+    
+    ext_shp_rcd[yr] = extremum[yr].shapeRecords()
+    max_rcd[yr] = max_speed[yr].records()
+        
+    # check if we are reaching the last two years or the last year
+    if (yr <= 2017):
+        ext_shp_rcd[yr+1] = extremum[yr + 1].shapeRecords()
+        max_rcd[yr+1] = max_speed[yr + 1].records()
+        
+        if (yr <= 2016):
+            ext_shp_rcd[yr+2] = extremum[yr + 2].shapeRecords()
+            max_rcd[yr+2] = max_speed[yr + 2].records()    
+    
+    # iterate on the vertices 
+    for n in g.vertices():
+        # Get the segment and the number of days
+        segment = g.vp.segment[n]      
+        
+        # calculate the indexes and dates
+        first = node_to_date_eddy(segment[0], e_overestim, d_init)
+        
+        # get days since 1950
+        first_days = first['days_1950']
+        
+        # get year
+        first_date = datetime.date(1950,1,1) + datetime.timedelta(first['days_1950'])
+        first_year = first_date.year
+        
+        #print(f'Segment size: {len(segment)}, first: {first_year}, last: {last_year}, yr: {yr}')        
+        
+        # chech if we are in the current year
+        
+        if (first_year == yr):
+            num_of_days = len(segment)
+            
+            # start processing
+            
+            first_pos = []
+            last_pos = []
+            
+            last = node_to_date_eddy(segment[-1], e_overestim, d_init)
+            last_days = last['days_1950']
+            
+            last_date = datetime.date(1950,1,1) + datetime.timedelta(last['days_1950'])
+            last_year = last_date.year
+            
+            print(f'Segment size: {len(segment)}, first: {first_year}, last: {last_year}, yr: {yr}')
+        
+            # calculate the location in the shapefile
+            
+            d_init_1 = ext_shp_rcd[first_year][0].record[1]
+            d_init_2 = ext_shp_rcd[last_year][0].record[1]
+            
+            first_loc = comp_ishape(d_init_1, ishape_last[first_year], first['days_1950'], first['eddy_index'])
+            last_loc = comp_ishape(d_init_2, ishape_last[last_year], last['days_1950'], last['eddy_index'])
+            
+            # grab the centers
+            
+            first_pos = ext_shp_rcd[first_year][first_loc].shape.points[0]
+            last_pos = ext_shp_rcd[last_year][last_loc].shape.points[0]
+            
+            ##### STORE POSITIONS IN THE VPS ######
+            g.vp.pos_first[n] = first_pos # [deg, deg]
+            g.vp.pos_last[n] = last_pos # [deg, deg]
+        
+            # if the segments are longer than the # of days over which to avg
+        
+            if (num_of_days > num_of_days_to_avg):
+                first_radii = 0 # [m]
+                last_radii = 0 # [m]
+                
+                first_rossby = 0 # []
+                last_rossby = 0 # []
+                
+                # First 7 days calculation
+                first_res = calculate_radii_and_rossby(0, num_of_days_to_avg, 1, segment, e_overestim, 
+                                                 d_init, ext_shp_rcd, max_rcd, ishape_last)      
+                
+                # average and assign radii
+                first_radii = first_res['radii'] / num_of_days_to_avg
+                g.vp.first_av_rad[n] = first_radii        
+                
+                # grab the days modifier
+                modifier = first_res['days_modifier']
+                
+                if (num_of_days_to_avg - modifier > 0):
+                    # Average and assign the rossbies:
+                    first_rossby = first_res['rossby'] / (num_of_days_to_avg - modifier)
+                    g.vp.first_av_ros[n] = first_rossby
+                else:
+                    # there is division by zero, average rossby is undefinied
+                    pass
+                    
+                # Last 7 days calculation
+                last_res = calculate_radii_and_rossby(len(segment) - 1, 
+                                                      len(segment) - (num_of_days_to_avg + 1), -1,
+                                                      segment, e_overestim, d_init, ext_shp_rcd, 
+                                                      max_rcd, ishape_last)                  
+                    
+                # Average and assign the last radii
+                last_radii = last_res['radii'] / num_of_days_to_avg
+                g.vp.last_av_rad[n] = last_radii        
+                
+                
+                # grab the days modifier
+                modifier = last_res['days_modifier']
+                
+                if (num_of_days_to_avg - modifier > 0):
+                    # Average and assign the rossbies:
+                    last_rossby = last_res['rossby'] / (num_of_days_to_avg - modifier)
+                    g.vp.last_av_ros[n] = last_rossby
+                else:
+                    # there is division by zero, average rossby is undefinied
+                    pass
+            # else, the number of eddies in a segment is lower than the # of 
+            # days over which to average, the values will be the same except
+            # for the positions
+            else:
+                res = calculate_radii_and_rossby(0, num_of_days, 1, segment, e_overestim, 
+                                                 d_init, ext_shp_rcd, max_rcd, ishape_last)              
+                
+                # grab the days modifier
+                modifier = res['days_modifier']
+                
+                if (num_of_days - modifier > 0):
+                    # Average and assign the rossbies:
+                    rossby = res['rossby'] / (num_of_days - modifier)
+                    g.vp.first_av_ros[n] = rossby
+                    g.vp.last_av_ros[n] = rossby
+                else:
+                    # there is division by zero, average rossby is undefinied
+                    pass
+                    
+                # Average and assign the radii
+                    
+                radii = res['radii'] / num_of_days
+                g.vp.first_av_rad[n] = radii
+                g.vp.last_av_rad[n] = radii
+            
+###############################
+# Calculate the cost function #
+###############################
+
+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
+
+    # calculate Delta_cent: numbers used for conversion obtained from:
+        # https://stackoverflow.com/questions/1253499/simple-calculations-for-working-with-lat-lon-and-km-distance
+    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:
+        print("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.nl_cost_function[edge] = cf    
+
+
+################################
+# Writing to an GT Graph File #
+################################
+
+#g.save("segments_cyclo_gl_cf.xml")
+#print("Done writing xml.")
+g.save(f'segmented_{orientation}_cf.gt')
+print("Done writing gt.")
+g.save(f'segmented_{orientation}_cf.xml')
+print("Done writing xml.")
+#g.save(f'{orientation}_segmented_cf.gv', fmt = "dot")
+#print("Done writing dot file.")
+
+print('All done')