Merge branch '0.4.x' into main

.github/pull_request_template.md

@@ -22,7 +22,7 @@ This includes a collection of test files which should be run according to which
 >
 > The files which have been changed during this PR can be listed using the command
 
-    git diff --name-only 0.3.x
+    git diff --name-only 0.4.x
 
 - [ ] My changes require one or more test files to be updated for all regression tests to pass.   
 
@@ -38,6 +38,6 @@ This includes a collection of test files which should be run according to which
 - [ ] I have commented my code, particularly in hard-to-understand areas.  
 - [ ] I have updated the documentation of the codebase where required.  
 - [ ] My changes generate no new warnings.   
-- [ ] My PR has been made to the `0.3.x` branch (**DO NOT SUBMIT A PR TO MAIN**)  
+- [ ] My PR has been made to the `0.4.x` branch (**DO NOT SUBMIT A PR TO MAIN**)  
 
 

.gitignore

@@ -27,3 +27,5 @@ docs_builder/*
 testing.ipynb
 docs/build/*
 dev_venv/*
+*.ipynb
+*.log

polar_route/__init__.py

@@ -1,4 +1,4 @@
-__version__ = "0.3.12"
+__version__ = "0.4.0"
 __description__ = "PolarRoute: Long-distance maritime polar route planning taking into account complex changing environmental conditions"
 __license__ = "MIT"
 __author__ = "Autonomous Marine Operations Planning (AMOP) Team, AI Lab, British Antarctic Survey"

polar_route/cli.py

@@ -152,7 +152,7 @@ def optimise_routes_cli():
     if args.dijkstra:
         # Form a unique name for the dijkstra output
         dijkstra_output_file_strs = output_file_strs
-        dijkstra_output_file_strs[0] += '_dijkstra'
+        dijkstra_output_file_strs[-2] += '_dijkstra'
 
         dijkstra_output_file = '.'.join(dijkstra_output_file_strs)
         logging.info(f"\tOutputting dijkstra path to {dijkstra_output_file}")

polar_route/crossing_smoothing.py

@@ -1047,8 +1047,8 @@ def blocked(self,new_cell,cell_a,cell_b):
         end   = cell_b['SIC']
         max_new = new_cell['SIC']
 
-        percentage_diff1  = (max_new-start)
-        percentage_diff2  = (max_new-end)
+        percentage_diff1  = (max_new-start)*100
+        percentage_diff2  = (max_new-end)*100
 
         if (percentage_diff1 <= self.blocked_sic*start) or (percentage_diff2 <= self.blocked_sic*end) or (max_new<=self.blocked_sic):
             return False

polar_route/route_calc.py

@@ -3,7 +3,7 @@
 import numpy as np
 import pandas as pd
 import geopandas as gpd
-from shapely import wkt
+from shapely import wkt, distance
 from shapely.geometry import Point, LineString, MultiLineString, Polygon
 from polar_route.utils import gpx_route_import
 
@@ -76,7 +76,12 @@ def traveltime_distance(cellbox, wp, cp, speed='speed', vector_x='uC', vector_y=
     su = cellbox[vector_x]
     sv = cellbox[vector_y]
     ssp = cellbox[speed][idx] * (1000 / (60 * 60))
-    traveltime, distance = traveltime_in_cell(x, y, su, sv, ssp)
+    try:
+        traveltime, distance = traveltime_in_cell(x, y, su, sv, ssp)
+    except:
+        traveltime=0
+        distance=0
+
     return traveltime, distance
 
 
@@ -269,29 +274,27 @@ def order_track(df, track_points):
 
     # Loop through crossing points to order them into a track along the route
     while pathing:
-        try:
-            start_point_segment = track_points['startPoints'].iloc[track_id]
-            end_point_segment = track_points['endPoints'].iloc[track_id]
-            path_point.append(start_point_segment)
-            cell_ids.append(track_points['cellID'].iloc[track_id])
-
-            if len(track_points['midPoints'].iloc[track_id]) != 0:
-                for midpnt in track_points['midPoints'].iloc[track_id]:
-                    path_point.append(midpnt)
-                    cell_ids.append(track_points['cellID'].iloc[track_id])
-
-            if end_point_segment == end_point:
-                pathing = False
-            track_id = np.where(track_points['startPoints'] == end_point_segment)[0][0]
-        except IndexError:
+        start_point_segment = track_points['startPoints'].iloc[track_id]
+        end_point_segment = track_points['endPoints'].iloc[track_id]
+        path_point.append(start_point_segment)
+        cell_ids.append(track_points['cellID'].iloc[track_id])
+
+        if len(track_points['midPoints'].iloc[track_id]) != 0:
+            for midpnt in track_points['midPoints'].iloc[track_id]:
+                path_point.append(midpnt)
+                cell_ids.append(track_points['cellID'].iloc[track_id])
+
+        if  distance(end_point_segment,end_point) < 0.05:
             pathing = False
-            path_point.append(end_point_segment)
-            cell_ids.append('NaN')
+        else:
+            track_id     = np.argmin([distance(entry,end_point_segment) for entry in track_points['startPoints']])
+            track_misfit = min([distance(entry,end_point_segment) for entry in track_points['startPoints']])
+            if track_misfit >= 0.05:
+                raise Exception('Path Segmentment not adding - ID={},Misfit={},distance from end={}'.format(track_id,track_misfit,distance(end_point_segment,end_point)))
 
     user_track = pd.DataFrame({'Point': path_point, 'CellID': cell_ids})
     return user_track
 
-
 def route_calc(route_file, mesh_file):
     """
         Function to calculate the fuel/time cost of a user defined route in a given mesh

polar_route/route_planner.py

@@ -10,14 +10,17 @@
 from shapely import wkt, Point, LineString, STRtree, Polygon
 import geopandas as gpd
 import logging
+from io import StringIO  
 
-from pandas.core.common import SettingWithCopyWarning
-warnings.simplefilter(action="ignore", category=SettingWithCopyWarning)
+# Squelching SettingWithCopyWarning 
+pd.options.mode.chained_assignment = None
 
 from polar_route.crossing import NewtonianDistance
 from polar_route.crossing_smoothing import Smoothing,PathValues,find_edge
 from polar_route.config_validation.config_validator import validate_route_config
 from polar_route.config_validation.config_validator import validate_waypoints
+from meshiphi import Boundary
+from meshiphi.utils import longitude_domain
 
 def _flattenCases(id,mesh):
     neighbour_case = []
@@ -134,11 +137,10 @@ def _mesh_boundary_polygon(mesh):
     lat_max = mesh['config']['mesh_info']['region']['lat_max']+tiny_value
     long_min = mesh['config']['mesh_info']['region']['long_min']-tiny_value
     long_max = mesh['config']['mesh_info']['region']['long_max']+tiny_value
-    p1 = Point([long_min, lat_min])
-    p2 = Point([long_min, lat_max])
-    p3 = Point([long_max, lat_max])
-    p4 = Point([long_max, lat_min])
-    return Polygon([p1,p2,p3,p4])
+
+    bounds = Boundary([lat_min, lat_max], [long_min, long_max])
+
+    return bounds.to_polygon()
 
 def _adjust_waypoints(point, cellboxes, max_distance=5):
     '''
@@ -311,8 +313,8 @@ def __init__(self, mesh, config, waypoints, cost_func=NewtonianDistance):
 
             adjusted_point = _adjust_waypoints(point, self.mesh['cellboxes'])
 
-            waypoints_df['Long'][idx] = adjusted_point.x
-            waypoints_df['Lat'][idx] = adjusted_point.y
+            waypoints_df.loc[idx, 'Long'] = adjusted_point.x
+            waypoints_df.loc[idx, 'Lat'] = adjusted_point.y
 
         source_waypoints_df   = waypoints_df[waypoints_df['Source'] == "X"]
         des_waypoints_df      = waypoints_df[waypoints_df['Destination'] == "X"]
@@ -328,7 +330,15 @@ def __init__(self, mesh, config, waypoints, cost_func=NewtonianDistance):
         self.smoothed_paths = None
         self.dijkstra_info = {}
 
+
+
+
+
+
         # ====== Loading Mesh & Neighbour Graph ======
+        # Zeroing currents if vectors names are not defined or zero_currents is defined
+        self.mesh = self._zero_currents(self.mesh)
+
         # Formatting the Mesh and Neighbour Graph to the right form
         self.neighbour_graph = pd.DataFrame(self.mesh['cellboxes']).set_index('id')
         self.neighbour_graph['geometry'] = self.neighbour_graph['geometry'].apply(wkt.loads)
@@ -355,6 +365,7 @@ def __init__(self, mesh, config, waypoints, cost_func=NewtonianDistance):
         if 'SIC' not in self.neighbour_graph:
             self.neighbour_graph['SIC'] = 0.0
 
+
         # ====== Objective Function Information ======
         #  Checking if objective function is in the dijkstra            
         if self.config['objective_function'] != 'traveltime':
@@ -406,14 +417,14 @@ def __init__(self, mesh, config, waypoints, cost_func=NewtonianDistance):
             if len(indices) > 1:
                 raise Exception('Waypoint lies in multiple cell boxes. Please check mesh ! ')
             else:
-                wpts['index'].loc[idx] = int(indices[0])
+                wpts.loc[idx, 'index'] = int(indices[0])
 
         self.mesh['waypoints'] = wpts[~wpts['index'].isnull()]
         self.mesh['waypoints']['index'] = self.mesh['waypoints']['index'].astype(int)
         self.mesh['waypoints'] =  self.mesh['waypoints'].to_json()
 
         # ==== Printing Configuration and Information
-        self.mesh['waypoints'] =  pd.read_json(self.mesh['waypoints'])
+        self.mesh['waypoints'] =  pd.read_json(StringIO(self.mesh['waypoints']))
 
         # # ===== Running the route planner for the given information
         # if ("dijkstra_only" in self.config) and self.config['dijkstra_only']:
@@ -431,6 +442,35 @@ def __init__(self, mesh, config, waypoints, cost_func=NewtonianDistance):
         # else:
         #     self.output = output
 
+    def _zero_currents(self,mesh):
+        '''
+            Applying zero currents to mesh
+
+            Input 
+                mesh (JSON) - MeshiPhi Mesh input
+            Output:
+                mesh (JSON) - MeshiPhi Mesh Corrected
+        '''
+
+        # Zeroing currents if both vectors are defined and zeroed
+        if ('zero_currents' in self.config) and ("vector_names" in self.config):
+            if self.config['zero_currents']:
+                logging.info('Zero Currents for Mesh !')
+                for idx,cell in enumerate(mesh['cellboxes']):
+                    cell[self.config['vector_names'][0]] = 0.0
+                    cell[self.config['vector_names'][1]] = 0.0
+                    mesh['cellboxes'][idx] = cell
+
+        # If no vectors are defined then add zero currents to mesh
+        if 'vector_names' not in self.config:
+            self.config['vector_names'] = ['Vector_x','Vector_y']
+            logging.info('No vector_names defined in config. Zeroing currents in mesh !')
+            for idx,cell in enumerate(mesh['cellboxes']):
+                cell[self.config['vector_names'][0]] = 0.0
+                cell[self.config['vector_names'][1]] = 0.0
+                mesh['cellboxes'][idx] = cell    
+
+        return mesh
 
 
     def to_json(self):
@@ -541,6 +581,8 @@ def _dijkstra_paths(self, start_waypoints, end_waypoints):
                         path['geometry'] = {}
                         path['geometry']['type'] = "LineString"
                         path_points = (np.array(wpt_a_loc+list(np.array(graph['pathPoints'].loc[wpt_b_index])[:-1, :])+wpt_b_loc))
+                        # Ensure all coordinates are in domain -180:180
+                        path_points[:,0] = longitude_domain(path_points[:,0])
                         path['geometry']['coordinates'] = path_points.tolist()
 
                         path['properties'] = {}
@@ -671,7 +713,7 @@ def _dijkstra(self, wpt_name):
         wpts = self.mesh['waypoints'][self.mesh['waypoints']['Name'].isin(self.end_waypoints)]
 
         # Initialising zero traveltime at the source location
-        source_index = int(self.mesh['waypoints'][self.mesh['waypoints']['Name'] == wpt_name]['index'])
+        source_index = int(self.mesh['waypoints'][self.mesh['waypoints']['Name'] == wpt_name]['index'].iloc[0])
 
         for vrbl in self.config['path_variables']:
             self.dijkstra_info[wpt_name].loc[source_index, 'shortest_{}'.format(vrbl)] = 0.0
@@ -785,6 +827,8 @@ def compute_smoothed_routes(self,blocked_metric='SIC',debugging=False):
             # Given a smoothed route path now determine the along path parameters.
             pv             = PathValues()
             path_info      = pv.objective_function(sf.aps,sf.start_waypoint,sf.end_waypoint)
+            # Ensure all coordinates are in domain -180:180
+            path_info['path'][:,0] = longitude_domain(path_info['path'][:,0])
             variables      = path_info['variables']
             TravelTimeLegs = variables['traveltime']['path_values']
             DistanceLegs   = variables['distance']['path_values'] 
@@ -811,6 +855,8 @@ def compute_smoothed_routes(self,blocked_metric='SIC',debugging=False):
             SmoothedPath['properties']['speed']      = list(SpeedLegs)
             SmoothedPaths += [SmoothedPath]
 
+            logging.info('{} iterations'.format(sf.jj))
+
         geojson['type'] = "FeatureCollection"
         geojson['features'] = SmoothedPaths
         self.smoothed_paths = geojson

polar_route/vessel_performance/vessel_performance_modeller.py

@@ -2,6 +2,7 @@
 import logging
 
 from meshiphi.mesh_generation.environment_mesh import EnvironmentMesh
+from meshiphi.mesh_generation.direction import Direction
 from polar_route.vessel_performance.vessel_factory import VesselFactory
 from polar_route.config_validation.config_validator import validate_vessel_config
 from polar_route.utils import timed_call
@@ -41,6 +42,9 @@ def model_accessibility(self):
             self.env_mesh.update_cellbox(i, access_values)
         inaccessible_nodes = [c.id for c in self.env_mesh.agg_cellboxes if c.agg_data['inaccessible']]
         logging.info(f"Found {len(inaccessible_nodes)} inaccessible cells in the mesh")
+        # Split any cells that neighbour inaccessible cells to match their size
+        self.split_neighbouring_cells(inaccessible_nodes)
+        # Remove inaccessible cells from graph
         for in_node in inaccessible_nodes:
             self.env_mesh.neighbour_graph.remove_node_and_update_neighbours(in_node)
 
@@ -65,8 +69,9 @@ def to_json(self):
             Returns:
                 j_mesh (dict): a dictionary representation of the modified mesh.
         """
-        self.env_mesh.config['vessel_info'] = self.config
+
         j_mesh = self.env_mesh.to_json()
+        j_mesh['config']['vessel_info'] = self.config
         return j_mesh
 
     def filter_nans(self):
@@ -76,4 +81,46 @@ def filter_nans(self):
         for i, cellbox in enumerate(self.env_mesh.agg_cellboxes):
             if any(np.isnan(val) for val in cellbox.agg_data.values() if type(val) == float):
                 filtered_data = {k: 0. if np.isnan(v) else v for k, v in cellbox.agg_data.items() if type(v) == float}
-                self.env_mesh.update_cellbox(i, filtered_data)
+                self.env_mesh.update_cellbox(i, filtered_data)
+
+    def split_neighbouring_cells(self, inaccessible_nodes):
+        """
+        Method to split any accessible cells that neighbour inaccessible cells until their sizes match
+
+        Args:
+            inaccessible_nodes (list): List of inaccessible nodes to split around
+
+        """
+        for in_node in inaccessible_nodes:
+            in_cb = self.env_mesh.get_cellbox(in_node)
+            neighbour_nodes = self.get_all_neighbours(in_node)
+            neighbours = [self.env_mesh.get_cellbox(nn) for nn in neighbour_nodes]
+            # Only interested in splitting accessible neighbours
+            acc_neighbours = [n for n in neighbours if not n.agg_data['inaccessible']]
+            # Split neighbouring cells until size matches the inaccessible cell
+            while any(neighbour.boundary.get_width() > in_cb.boundary.get_width() or
+                      neighbour.boundary.get_height() > in_cb.boundary.get_height() for neighbour in acc_neighbours):
+                # Split all larger neighbours
+                for neighbour in acc_neighbours:
+                    if neighbour.boundary.get_width() > in_cb.boundary.get_width():
+                        self.env_mesh.split_and_replace(str(neighbour.id))
+                # Extract new neighbours after splitting
+                neighbour_nodes = self.get_all_neighbours(in_node)
+                neighbours = [self.env_mesh.get_cellbox(nn) for nn in neighbour_nodes]
+                acc_neighbours = [n for n in neighbours if not n.agg_data['inaccessible']]
+
+    def get_all_neighbours(self, cell_id):
+        """
+        Method to get a list of all neighbouring cell ids for a particular cell within the environmental mesh
+
+        Args:
+            cell_id (str): Cell ID to find the neighbours around
+        Returns:
+            neighbours (list): List of IDs of neighbouring cells
+        """
+        neighbours = []
+        direction_obj = Direction()
+        for direction in direction_obj.__dict__.values():
+            neighbours.extend(self.env_mesh.neighbour_graph.get_neighbours(cell_id, str(direction)))
+
+        return neighbours

requirements.txt

@@ -5,7 +5,7 @@ jsonschema
 matplotlib
 netcdf4
 numpy>=1.21.6
-pandas>=1.3.5,<1.5
+pandas>=1.3.5
 pyproj
 pytest
 shapely