1+ """
2+ build_global_trajectories.py (v2)
3+ ====================================
4+ 新增字段：
5+ - day : 代表性日期 (1-31)
6+ - time_of_day : 代表性时刻，单位分钟，取整到10分钟 (0-1430)
7+ """
8+ import json , os , math
9+ import pandas as pd
10+ import numpy as np
11+
12+ CLEANED_FILE = r"E:\datavisul\ais_output\cleaned_data\ais_cleaned_with_vesseltype.parquet"
13+ OUT_FILE = os .path .join (os .path .dirname (os .path .abspath (__file__ )), "vessel_trajectories_global.json" )
14+
15+ MAX_POINTS_PER_SEGMENT = 20
16+ MIN_POINTS = 3
17+ SAMPLE_RG_FRAC = 0.20
18+
19+ VESSEL_TYPE_COLORS = {
20+ "Cargo" : "#2196F3" , "Cargo - Hazard A" : "#1565C0" , "Cargo - Hazard B" : "#1976D2" ,
21+ "Cargo - Hazard C" : "#1E88E5" , "Cargo - Hazard D" : "#42A5F5" ,
22+ "Tanker" : "#FF5722" , "Tanker - Hazard A" : "#E64A19" , "Tanker - Hazard B" : "#FF7043" ,
23+ "Tanker - Hazard C" : "#FF8A65" , "Tanker - Hazard D" : "#FFAB91" ,
24+ "Fishing" : "#4CAF50" , "Tug" : "#FF9800" , "Passenger" : "#9C27B0" ,
25+ "Passenger - Hazard B" : "#AB47BC" , "Sailing" : "#00BCD4" ,
26+ "Pleasure craft" : "#8BC34A" , "Search and rescue" : "#F44336" ,
27+ "Dredging or underwater ops" : "#795548" , "Towing" : "#FF6F00" ,
28+ "Towing astern" : "#FFA000" , "Military ops" : "#607D8B" ,
29+ "Law enforcement" : "#455A64" , "Pilot vessel" : "#26C6DA" ,
30+ "High-speed craft" : "#00E5FF" , "Anti-pollution" : "#69F0AE" ,
31+ "Diving ops" : "#40C4FF" , "Port tender" : "#B0BEC5" ,
32+ "Spare local vessel" : "#78909C" , "Noncombatant ship" : "#546E7A" ,
33+ "Medical transport" : "#EF9A9A" , "Wing in ground" : "#CE93D8" ,
34+ "Other" : "#9E9E9E" , "Unknown" : "#607D8B" ,
35+ }
36+
37+ def classify_pattern (mean_sog , total_km ):
38+ if mean_sog < 0.3 :
39+ return "anchored"
40+ if total_km < 20 or mean_sog < 4 :
41+ return "coastal"
42+ return "inbound"
43+
44+ def haversine_km (lat1 , lon1 , lat2 , lon2 ):
45+ R = 6371
46+ r = math .pi / 180
47+ dlat = (lat2 - lat1 ) * r
48+ dlon = (lon2 - lon1 ) * r
49+ a = math .sin (dlat / 2 )** 2 + math .cos (lat1 * r )* math .cos (lat2 * r )* math .sin (dlon / 2 )** 2
50+ return 2 * R * math .asin (math .sqrt (max (0 , a )))
51+
52+ def extract_time_fields (dt_series ):
53+ """
54+ 从 datetime Series 提取代表性的 month, day, time_of_day。
55+ time_of_day: 小时*60 + 分钟，取整到10分钟，范围 0~1430。
56+ 返回 (month, day, time_of_day)，失败时返回 (1, 1, 0)。
57+ """
58+ valid = dt_series .dropna ()
59+ if len (valid ) == 0 :
60+ return 1 , 1 , 0
61+ month = int (valid .dt .month .mode ().iloc [0 ])
62+ day = int (valid .dt .day .mode ().iloc [0 ])
63+ hour = int (valid .dt .hour .mode ().iloc [0 ])
64+ minute = int (valid .dt .minute .mode ().iloc [0 ])
65+ minute = (minute // 10 ) * 10
66+ time_of_day = hour * 60 + minute
67+ return month , day , time_of_day
68+
69+ print (f"Reading { CLEANED_FILE } )
70+
71+ import pyarrow .parquet as _pq
72+ import gc as _gc
73+
74+ _NEED = {"MMSI" , "LAT" , "LON" , "SOG" , "BaseDateTime" , "VesselType" , "VesselTypeLabel" , "TrackSegmentID" }
75+ _pf = _pq .ParquetFile (CLEANED_FILE )
76+ n_rg = _pf .metadata .num_row_groups
77+ _available = set (_pf .schema_arrow .names )
78+ _read_cols = sorted (_NEED & _available )
79+
80+ rng = np .random .default_rng (42 )
81+ selected_rgs = sorted (rng .choice (n_rg , max (1 , int (round (n_rg * SAMPLE_RG_FRAC ))), replace = False ))
82+ print (f" { n_rg } { len (selected_rgs )} { SAMPLE_RG_FRAC * 100 :.0f} )
83+ print (f" Columns: { _read_cols } )
84+
85+ chunks = []
86+ for i , rg_idx in enumerate (selected_rgs ):
87+ if i % 20 == 0 :
88+ print (f" Reading row-group batch { i + 1 } { len (selected_rgs )} , flush = True )
89+ chunks .append (_pf .read_row_group (rg_idx , columns = _read_cols ).to_pandas ())
90+
91+ df = pd .concat (chunks , ignore_index = True )
92+ del chunks , _pf
93+ _gc .collect ()
94+
95+ required = {"MMSI" , "LAT" , "LON" , "SOG" , "BaseDateTime" }
96+ missing = required - set (df .columns )
97+ if missing :
98+ raise ValueError (f"Missing columns: { missing } )
99+
100+ if "VesselType" in df .columns and "VesselTypeLabel" in df .columns :
101+ df .loc [df ["VesselType" ] == 36 , "VesselTypeLabel" ] = "Sailing"
102+ df .loc [df ["VesselType" ] == 36.0 , "VesselTypeLabel" ] = "Sailing"
103+ df .loc [df ["VesselType" ] == 37 , "VesselTypeLabel" ] = "Pleasure craft"
104+ df .loc [df ["VesselType" ] == 37.0 , "VesselTypeLabel" ] = "Pleasure craft"
105+
106+ if "VesselTypeLabel" in df .columns :
107+ df ["VesselTypeLabel" ] = df ["VesselTypeLabel" ].astype ("category" )
108+
109+ has_segment_id = "TrackSegmentID" in df .columns
110+ print (f"TrackSegmentID available: { has_segment_id } )
111+ print (f"Loaded { len (df ):,} { df ['MMSI' ].nunique ():,} )
112+
113+ df ["BaseDateTime" ] = pd .to_datetime (df ["BaseDateTime" ], utc = True , errors = "coerce" )
114+ df = df .dropna (subset = ["LAT" , "LON" ])
115+
116+ vtype_col = None
117+ for c in ["VesselTypeLabel" , "vessel_type" , "VesselType" ]:
118+ if c in df .columns :
119+ vtype_col = c
120+ break
121+
122+ trajectories = []
123+
124+ if has_segment_id :
125+ group_keys = ["MMSI" , "TrackSegmentID" ]
126+ groups = df .groupby (group_keys , sort = False )
127+ total_groups = len (groups )
128+ print (f"Total (MMSI, TrackSegmentID) groups: { total_groups :,} )
129+
130+ for i , ((mmsi , seg_id ), grp ) in enumerate (groups ):
131+ if i % 2000 == 0 :
132+ print (f" Segment { i + 1 } { total_groups } )
133+
134+ if len (grp ) < MIN_POINTS :
135+ continue
136+ if len (grp ) > MAX_POINTS_PER_SEGMENT :
137+ idx = np .round (np .linspace (0 , len (grp ) - 1 , MAX_POINTS_PER_SEGMENT )).astype (int )
138+ grp = grp .iloc [idx ]
139+
140+ vtype = "Unknown"
141+ if vtype_col :
142+ mode_vals = grp [vtype_col ].dropna ()
143+ if len (mode_vals ):
144+ vtype = str (mode_vals .mode ().iloc [0 ])
145+
146+ color = VESSEL_TYPE_COLORS .get (vtype , "#9E9E9E" )
147+ mean_sog = float (grp ["SOG" ].mean ()) if grp ["SOG" ].notna ().any () else 0.0
148+ first , last = grp .iloc [0 ], grp .iloc [- 1 ]
149+ total_km = haversine_km (float (first ["LAT" ]), float (first ["LON" ]),
150+ float (last ["LAT" ]), float (last ["LON" ]))
151+ pattern = classify_pattern (mean_sog , total_km )
152+
153+ month , day , time_of_day = extract_time_fields (grp ["BaseDateTime" ])
154+
155+ pts = [
156+ {"lat" : round (float (r ["LAT" ]), 5 ),
157+ "lon" : round (float (r ["LON" ]), 5 ),
158+ "sog" : round (float (r ["SOG" ]) if pd .notna (r ["SOG" ]) else 0.0 , 1 )}
159+ for _ , r in grp .iterrows ()
160+ ]
161+
162+ if len (pts ) >= 2 :
163+ trajectories .append ({
164+ "mmsi" : str (mmsi ),
165+ "vessel_type" : vtype ,
166+ "color" : color ,
167+ "pattern" : pattern ,
168+ "month" : month ,
169+ "day" : day , # 新增
170+ "time_of_day" : time_of_day , # 新增
171+ "points" : pts ,
172+ })
173+
174+ else :
175+ print ("TrackSegmentID not found — reconstructing segments from time gaps ..." )
176+ groups = df .groupby ("MMSI" , sort = False )
177+ total_vessels = len (groups )
178+
179+ for i , (mmsi , grp ) in enumerate (groups ):
180+ if i % 500 == 0 :
181+ print (f" Vessel { i + 1 } { total_vessels } )
182+
183+ grp = grp .sort_values ("BaseDateTime" ).reset_index (drop = True )
184+ if len (grp ) < MIN_POINTS :
185+ continue
186+
187+ dt_h = grp ["BaseDateTime" ].diff ().dt .total_seconds ().fillna (0 ) / 3600
188+ lats , lons = grp ["LAT" ].values , grp ["LON" ].values
189+ dkm = np .array ([0.0 ] + [
190+ haversine_km (lats [j - 1 ], lons [j - 1 ], lats [j ], lons [j ])
191+ for j in range (1 , len (grp ))
192+ ])
193+ implied_spd = np .where (dt_h > 0 , dkm / dt_h , 0 )
194+
195+ breaks = np .where ((dt_h > 6 ) | (implied_spd > 100 ))[0 ]
196+ seg_starts = [0 ] + list (breaks )
197+ seg_ends = list (breaks ) + [len (grp )]
198+
199+ for s , e in zip (seg_starts , seg_ends ):
200+ seg = grp .iloc [s :e ]
201+ if len (seg ) < MIN_POINTS :
202+ continue
203+ if len (seg ) > MAX_POINTS_PER_SEGMENT :
204+ idx = np .round (np .linspace (0 , len (seg ) - 1 , MAX_POINTS_PER_SEGMENT )).astype (int )
205+ seg = seg .iloc [idx ]
206+
207+ vtype = "Unknown"
208+ if vtype_col :
209+ mode_vals = seg [vtype_col ].dropna ()
210+ if len (mode_vals ):
211+ vtype = str (mode_vals .mode ().iloc [0 ])
212+
213+ color = VESSEL_TYPE_COLORS .get (vtype , "#9E9E9E" )
214+ mean_sog = float (seg ["SOG" ].mean ()) if seg ["SOG" ].notna ().any () else 0.0
215+ first , last = seg .iloc [0 ], seg .iloc [- 1 ]
216+ total_km = haversine_km (float (first ["LAT" ]), float (first ["LON" ]),
217+ float (last ["LAT" ]), float (last ["LON" ]))
218+ pattern = classify_pattern (mean_sog , total_km )
219+
220+ month , day , time_of_day = extract_time_fields (seg ["BaseDateTime" ])
221+
222+ pts = [
223+ {"lat" : round (float (r ["LAT" ]), 5 ),
224+ "lon" : round (float (r ["LON" ]), 5 ),
225+ "sog" : round (float (r ["SOG" ]) if pd .notna (r ["SOG" ]) else 0.0 , 1 )}
226+ for _ , r in seg .iterrows ()
227+ ]
228+
229+ if len (pts ) >= 2 :
230+ trajectories .append ({
231+ "mmsi" : str (mmsi ),
232+ "vessel_type" : vtype ,
233+ "color" : color ,
234+ "pattern" : pattern ,
235+ "month" : month ,
236+ "day" : day , # 新增
237+ "time_of_day" : time_of_day , # 新增
238+ "points" : pts ,
239+ })
240+
241+ print (f"\n Total trajectories: { len (trajectories ):,} )
242+
243+ MAX_TOTAL_TRAJ = 60000
244+ if len (trajectories ) > MAX_TOTAL_TRAJ :
245+ rng2 = np .random .default_rng (123 )
246+ keep = sorted (rng2 .choice (len (trajectories ), MAX_TOTAL_TRAJ , replace = False ).tolist ())
247+ trajectories = [trajectories [i ] for i in keep ]
248+ print (f"Capped to { MAX_TOTAL_TRAJ :,} )
249+
250+ with open (OUT_FILE , "w" ) as f :
251+ json .dump (trajectories , f , separators = ("," , ":" ))
252+
253+ size_mb = os .path .getsize (OUT_FILE ) / 1e6
254+ print (f"Written to { OUT_FILE } { size_mb :.1f} )
255+ print ("Done." )
0 commit comments