1+ rule create_maritime_network :
2+ input :
3+ nodes = "{OUTPUT_DIR}/input/networks/maritime/nodes.gpq" ,
4+ edges_no_geom = "{OUTPUT_DIR}/input/networks/maritime/edges_by_cargo/maritime_base_network_general_cargo.pq" ,
5+ edges_visualisation = "{OUTPUT_DIR}/input/networks/maritime/edges.gpq" ,
6+ output :
7+ nodes = "{OUTPUT_DIR}/maritime_network/nodes.gpq" ,
8+ edges = "{OUTPUT_DIR}/maritime_network/edges.gpq" ,
9+ run :
10+ import igraph as ig
11+ import geopandas as gpd
12+ from shapely .geometry import Point
13+ from shapely .ops import linemerge
14+ from tqdm import tqdm
15+
16+ from open_gira .network_creation import preprocess_maritime_network
17+
18+ # possible cargo types = ("container", "dry_bulk", "general_cargo", "roro", "tanker")
19+ # for now, just use 'general_cargo'
20+ maritime_nodes , maritime_edges_no_geom = preprocess_maritime_network (
21+ input .nodes ,
22+ input .edges_no_geom
23+ )
24+
25+ if config ["study_country_iso_a3" ] == "THA" :
26+ # put Bangkok port in the right place...
27+ maritime_nodes .loc [maritime_nodes .name == "Bangkok_Thailand" , "geometry" ] = Point ((100.5753 , 13.7037 ))
28+
29+ # %%
30+ # Jasper's maritime edges in 'edges_by_cargo' do not contain geometry
31+ # this is because the AIS data that they were derived from only contain origin and destination port, not route
32+ # this is a pain for visualisation, so we will create a geometry for each from `maritime_vis_edges`
33+
34+ maritime_vis_edges = gpd .read_parquet (input .edges_visualisation )
35+ vis_graph = ig .Graph .DataFrame (maritime_vis_edges , directed = True , use_vids = False )
36+
37+ maritime_edges = maritime_edges_no_geom .copy ()
38+ change_of_port_mask = maritime_edges_no_geom .from_port != maritime_edges_no_geom .to_port
39+ port_pairs_to_generate_geom_for = maritime_edges_no_geom [change_of_port_mask ]
40+ for index , row in tqdm (port_pairs_to_generate_geom_for .iterrows (), total = len (port_pairs_to_generate_geom_for )):
41+ edge_list = vis_graph .get_shortest_path (row .from_port , row .to_port , weights = "distance" , output = "epath" )
42+ route_edges = maritime_vis_edges .iloc [edge_list ]
43+ route_linestring = linemerge (list (route_edges .geometry ))
44+ maritime_edges .loc [index , "geometry" ] = route_linestring
45+
46+ maritime_edges = gpd .GeoDataFrame (maritime_edges ).set_crs (epsg = 4326 )
47+
48+ maritime_nodes .to_parquet (output .nodes )
49+ maritime_edges .to_parquet (output .edges )
50+
51+
52+ rule plot_maritime_network :
53+ input :
54+ nodes = "{OUTPUT_DIR}/maritime_network/nodes.gpq" ,
55+ edges = "{OUTPUT_DIR}/maritime_network/edges.gpq" ,
56+ output :
57+ edges_plot = "{OUTPUT_DIR}/maritime_network/edges.png" ,
58+ run :
59+ import geopandas as gpd
60+ import matplotlib
61+ import matplotlib .pyplot as plt
62+ import numpy as np
63+
64+ from open_gira .plot .utils import chop_at_antimeridian
65+
66+ matplotlib .use ("Agg" )
67+ plt .style .use ("bmh" )
68+
69+ world = gpd .read_file (gpd .datasets .get_path ('naturalearth_lowres' ))
70+ world .geometry = world .geometry .boundary
71+
72+ maritime_nodes = gpd .read_parquet (input .nodes )
73+ maritime_edges = gpd .read_parquet (input .edges )
74+
75+ # whole network
76+ f , ax = plt .subplots (figsize = (16 , 7 ))
77+ chop_at_antimeridian (maritime_edges , drop_null_geometry = True ).plot (
78+ ax = ax ,
79+ linewidth = 0.5 ,
80+ alpha = 0.8
81+ )
82+ world .plot (ax = ax , lw = 0.5 , alpha = 0.2 )
83+ ax .set_xticks (np .linspace (- 180 , 180 , 13 ))
84+ ax .set_yticks ([- 60 , - 30 , 0 , 30 , 60 ])
85+ ax .set_ylim (- 65 , 85 )
86+ ax .set_xlim (- 180 , 180 )
87+ ax .grid (alpha = 0.3 )
88+ ax .set_xlabel ("Longitude [deg]" )
89+ ax .set_ylabel ("Latitude [deg]" )
90+ f .savefig (output .edges_plot )
91+
92+
93+ rule plot_port_connections :
94+ input :
95+ nodes = "{OUTPUT_DIR}/maritime_network/nodes.gpq" ,
96+ edges = "{OUTPUT_DIR}/maritime_network/edges.gpq" ,
97+ output :
98+ port_trade_plots = directory ("{OUTPUT_DIR}/maritime_network/port_trade_plots" ),
99+ run :
100+ import os
101+
102+ import geopandas as gpd
103+ import matplotlib
104+ import matplotlib .pyplot as plt
105+ from tqdm import tqdm
106+
107+ from open_gira .plot .utils import chop_at_antimeridian
108+
109+ matplotlib .use ("Agg" )
110+ plt .style .use ("bmh" )
111+
112+ maritime_nodes = gpd .read_parquet (input .nodes )
113+ maritime_edges = gpd .read_parquet (input .edges )
114+
115+ # disambiguate the global view and plot the routes from each port, one port at a time
116+ world = gpd .read_file (gpd .datasets .get_path ('naturalearth_lowres' ))
117+ world .geometry = world .geometry .boundary
118+
119+ ports = maritime_edges .from_port .unique ()
120+ os .makedirs (output .port_trade_plots )
121+ for port_id in tqdm (ports ):
122+ filepath = os .path .join (output .port_trade_plots , f"{ port_id } )
123+ f , ax = plt .subplots (figsize = (10 ,10 ))
124+ port = maritime_nodes [maritime_nodes .id == f"{ port_id } ]
125+ routes = maritime_edges [maritime_edges .from_port == port_id ]
126+ if all (routes .geometry .isna ()):
127+ continue
128+ try :
129+ chop_at_antimeridian (routes , drop_null_geometry = True ).plot (
130+ column = "to_port" ,
131+ categorical = True ,
132+ ax = ax
133+ )
134+ except ValueError :
135+ print (f"Failed to plot { port_id } )
136+ continue
137+ maritime_nodes [maritime_nodes .id == f"{ port_id } ].plot (
138+ ax = ax ,
139+ markersize = 500 ,
140+ marker = "*" ,
141+ facecolor = "none" ,
142+ color = "r"
143+ )
144+ xmin , xmax = ax .get_xlim ()
145+ ymin , ymax = ax .get_ylim ()
146+ world .plot (ax = ax , linewidth = 0.5 , alpha = 0.4 )
147+ ax .set_xlim (xmin , xmax )
148+ ax .set_ylim (ymin , ymax )
149+ port_name , = port .name
150+ ax .set_title (f"{ port_id } { port_name .replace ('_' , ', ' )} )
151+ ax .get_xaxis ().set_visible (False )
152+ ax .get_yaxis ().set_visible (False )
153+
154+ f .savefig (filepath )
155+ plt .close (f )
0 commit comments