extract_route_between.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from tqdm import tqdm\n",
    "from multiprocessing import Pool\n",
    "from functools import partial\n",
    "import utils.data_preprocess_utils as dpu\n",
    "import geopandas as gpd\n",
    "from utils.data_preprocess_utils import get_config, get_all_files\n",
    "from easydict import EasyDict\n",
    "cfg = EasyDict(get_config('./cfg/route_extract_cfg.yaml'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "files = get_all_files(cfg.raw_data_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "import s2sphere\n",
    "from datetime import datetime, timedelta\n",
    "from concurrent.futures import ThreadPoolExecutor\n",
    "from tqdm import tqdm\n",
    "\n",
    "# Step 2: 创建S2 Cell分层器\n",
    "def get_s2_cells(latitudes, longitudes, level=9):\n",
    "    cell_ids = [s2sphere.CellId.from_lat_lng(s2sphere.LatLng.from_degrees(lat, lon)) for lat, lon in zip(latitudes, longitudes)]\n",
    "    return [cell_id.parent(level).to_token() for cell_id in cell_ids]\n",
    "\n",
    "import pandas as pd\n",
    "import s2sphere\n",
    "from concurrent.futures import ThreadPoolExecutor\n",
    "\n",
    "def process_ais_file(file_path):\n",
    "    # 读取单个AIS数据文件\n",
    "    ais_data = pd.read_csv(file_path, usecols=cfg.columns_to_extract, compression='gzip')\n",
    "    columns_to_check = ['longitude', 'latitude', 'speed', 'heading']\n",
    "    ais_data = ais_data.dropna(subset=columns_to_check).reset_index(drop=True)\n",
    "    ais_data = dpu.convert_time_format(ais_data, time_column_ori='timestamp', time_column_new='timestamp')\n",
    "    ais_data.drop_duplicates(subset=['mmsi', 'timestamp'], inplace=True)\n",
    "    latitudes = ais_data['latitude'].values\n",
    "    longitudes = ais_data['longitude'].values\n",
    "    speeds = ais_data['speed'].values  # 使用数据中的速度信息\n",
    "    mmsi_values = ais_data['mmsi'].values\n",
    "\n",
    "    # 标识停止点\n",
    "    is_stop_point = speeds < 0.2\n",
    "    # 关联S2单元\n",
    "    s2_cell_ids = get_s2_cells(latitudes, longitudes)\n",
    "    # 统计每个cell中停止的船只数量\n",
    "    stopped_ships_per_cell = {}\n",
    "    for i, mmsi in enumerate(mmsi_values):\n",
    "        if is_stop_point[i]:\n",
    "            s2_cell_id = s2_cell_ids[i]\n",
    "            stopped_ships_per_cell.setdefault(s2_cell_id, set())\n",
    "            stopped_ships_per_cell[s2_cell_id].add(mmsi)\n",
    "\n",
    "    return stopped_ships_per_cell\n",
    "\n",
    "# 多线程处理所有文件的停止的船只数量\n",
    "files = get_all_files(cfg.raw_data_dir)\n",
    "files = files\n",
    "total_stopped_ships_per_cell = {}\n",
    "\n",
    "# 使用tqdm显示进度\n",
    "with tqdm(total=len(files), desc=\"Processing files\", unit=\"file\") as pbar:\n",
    "    def process_file(file_path):\n",
    "        stopped_ships_per_cell = process_ais_file(file_path)\n",
    "        \n",
    "        # 合并当前文件的结果到总的结果中\n",
    "        for cell_id, mmsi_set in stopped_ships_per_cell.items():\n",
    "            total_stopped_ships_per_cell.setdefault(cell_id, set())\n",
    "            total_stopped_ships_per_cell[cell_id].update(mmsi_set)\n",
    "\n",
    "        pbar.update(1)  # 更新进度条\n",
    "\n",
    "    # 使用ThreadPoolExecutor进行多线程处理\n",
    "    with ThreadPoolExecutor(max_workers=12) as executor:\n",
    "        files = [file for file in files if file.endswith(\"csv.gz\")]\n",
    "        executor.map(process_file, files)\n",
    "\n",
    "# 打印每个cell中停止的船只数量\n",
    "for cell_id, mmsi_set in total_stopped_ships_per_cell.items():\n",
    "    print(f\"S2 Cell {cell_id}: {len(mmsi_set)} 停止的船只数量\")\n",
    "\n",
    "# 保存结果到文件\n",
    "result_df = pd.DataFrame([(cell_id, len(mmsi_set)) for cell_id, mmsi_set in total_stopped_ships_per_cell.items()], columns=['S2_Cell_ID', 'Stopped_Ships_Count'])\n",
    "result_df.to_csv('stopped_ships_counts_per_cell.csv', index=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "from shapely.geometry import Polygon\n",
    "import s2sphere\n",
    "\n",
    "result_df = pd.read_csv('stopped_ships_counts_per_cell.csv')\n",
    "\n",
    "# Create an empty GeoDataFrame\n",
    "gdf = gpd.GeoDataFrame()\n",
    "\n",
    "# Add geometries based on 'S2_Cell_ID'\n",
    "geometry_list = []\n",
    "for _, row in result_df.iterrows():\n",
    "    cell_id = s2sphere.CellId.from_token(row['S2_Cell_ID'])\n",
    "    cell = s2sphere.Cell(cell_id)\n",
    "    vertices = [s2sphere.LatLng.from_point(cell.get_vertex(i)) for i in range(4)]\n",
    "    polygon = Polygon([(v.lng().degrees, v.lat().degrees) for v in vertices])\n",
    "    geometry_list.append(polygon)\n",
    "\n",
    "# Add the 'geometry' column to the GeoDataFrame\n",
    "gdf['geometry'] = geometry_list\n",
    "\n",
    "# Combine the original DataFrame with the GeoDataFrame\n",
    "gdf = pd.concat([result_df, gdf], axis=1)\n",
    "gdf = gpd.GeoDataFrame(gdf)\n",
    "# Now 'gdf' should have the correct geometries\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save to file\n",
    "gdf.to_file('archorage_dict.geojson', driver='GeoJSON')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# filter anchorage with too few ships\n",
    "min_ships = 10\n",
    "total_stopped_ships_per_cell = {k: v for k, v in total_stopped_ships_per_cell.items() if len(v) > min_ships}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "archorage_dict = {}\n",
    "for k in total_stopped_ships_per_cell.keys():\n",
    "    archorage_dict[k] = {'ship_count': len(total_stopped_ships_per_cell[k])}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.cluster import DBSCAN\n",
    "from sklearn.metrics.pairwise import haversine_distances\n",
    "from math import radians\n",
    "import numpy as np\n",
    "def calculate_geodesic_distance_matrix(matrix1:np.ndarray, matrix2:np.ndarray):\n",
    "    lat1, lon1 = matrix1.T[0], matrix1.T[1]\n",
    "    lat2, lon2 = matrix2.T[0], matrix2.T[1]\n",
    "    \n",
    "    lat1 = np.radians(lat1)\n",
    "    lon1 = np.radians(lon1)\n",
    "    lat2 = np.radians(lat2)\n",
    "    lon2 = np.radians(lon2)\n",
    "    \n",
    "    # Earth radius in nautical miles\n",
    "    earth_radius = 3440.065\n",
    "    \n",
    "    dlat = lat2[:, np.newaxis] - lat1\n",
    "    dlon = lon2[:, np.newaxis] - lon1\n",
    "    \n",
    "    a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2[:, np.newaxis]) * np.sin(dlon/2)**2\n",
    "    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))\n",
    "    distance_matrix = earth_radius * c\n",
    "    \n",
    "    return distance_matrix\n",
    "# 假设 anchorage_dict 是一个字典，包含每个 anchorage 的 s2_cell_id 作为 key，船只数目作为 value\n",
    "anchorage_coordinates = []  # 用于存储 anchorage 的中心点坐标\n",
    "anchorage_ids = []  # 用于存储 anchorage 的 s2_cell_id\n",
    "\n",
    "for s2_cell_id, ship_count in total_stopped_ships_per_cell.items():\n",
    "    # 假设你有 anchorage 的位置信息，可以通过计算中心点坐标\n",
    "    cell_id = s2sphere.CellId.from_token(s2_cell_id)\n",
    "    cell = s2sphere.Cell(cell_id)\n",
    "    vertices = []\n",
    "    for i in range(4):\n",
    "        vertex = s2sphere.LatLng.from_point(cell.get_vertex(i))\n",
    "        vertices.append([vertex.lat().degrees, vertex.lng().degrees])\n",
    "\n",
    "    # 计算 anchorage 的中心点坐标\n",
    "    center_latitude = sum(coord[0] for coord in vertices) / 4\n",
    "    center_longitude = sum(coord[1] for coord in vertices) / 4\n",
    "\n",
    "    anchorage_coordinates.append([center_latitude, center_longitude])\n",
    "    anchorage_ids.append(s2_cell_id)\n",
    "\n",
    "# 转换成 numpy 数组\n",
    "anchorage_coordinates = np.array(anchorage_coordinates)\n",
    "\n",
    "# 计算 haversine 距离矩阵\n",
    "haversine_matrix = calculate_geodesic_distance_matrix(anchorage_coordinates, anchorage_coordinates)\n",
    "\n",
    "\n",
    "\n",
    "# 使用 DBSCAN 聚类\n",
    "dbscan = DBSCAN(eps=15, min_samples=1, metric='precomputed')\n",
    "cluster_labels = dbscan.fit_predict(haversine_matrix)\n",
    "\n",
    "# 将结果添加到 total_stopped_ships_per_cell 中\n",
    "for s2_cell_id, cluster_label in zip(anchorage_ids, cluster_labels):\n",
    "    archorage_dict[s2_cell_id]['cluster_label'] = cluster_label\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "files = get_all_files(cfg.raw_data_dir)\n",
    "\n",
    "for file_path in files:\n",
    "    ais_data = pd.read_csv(file_path, usecols=cfg.columns_to_extract, compression='gzip')\n",
    "    columns_to_check = ['longitude', 'latitude', 'speed', 'heading']\n",
    "    ais_data = ais_data.dropna(subset=columns_to_check).reset_index(drop=True)\n",
    "    ais_data = dpu.convert_time_format(ais_data, time_column_ori='timestamp', time_column_new='timestamp')\n",
    "    ais_data.drop_duplicates(subset=['mmsi', 'timestamp'], inplace=True)\n",
    "    break\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import s2sphere\n",
    "\n",
    "# 假设 ais_data 是包含船舶轨迹数据的 DataFrame，包括 'Latitude', 'Longitude', 'Timestamp' 等列\n",
    "# archorage_dict 是包含 archorage 信息的字典，key 是 s2_cell_id，value 是 cluster\n",
    "\n",
    "\n",
    "from shapely.geometry import Polygon, Point\n",
    "\n",
    "def get_archorage_status(row, archorage_dict, current_archorage):\n",
    "    lat = row['latitude']\n",
    "    lon = row['longitude']\n",
    "\n",
    "    for s2_cell_id, cluster in archorage_dict.items():\n",
    "        if cluster == current_archorage:\n",
    "            continue  # If the current trajectory is already in the same numbered archorage, skip\n",
    "\n",
    "        cell_id = s2sphere.CellId.from_token(s2_cell_id)\n",
    "        cell = s2sphere.Cell(cell_id)\n",
    "        vertices = [s2sphere.LatLng.from_point(cell.get_vertex(i)) for i in range(4)]\n",
    "\n",
    "        # Create a Polygon using the vertices\n",
    "        polygon = Polygon([(v.lng().degrees, v.lat().degrees) for v in vertices])\n",
    "\n",
    "        # Check if the point is inside the polygon\n",
    "        point = Point(lon, lat)\n",
    "        if point.within(polygon):\n",
    "            return cluster\n",
    "\n",
    "    return None\n",
    "\n",
    "\n",
    "def process_ais_trajectory(ais_data, archorage_dict):\n",
    "    # 初始化状态列\n",
    "    ais_data['Status'] = None\n",
    "    current_archorage = None\n",
    "    \n",
    "    # 遍历轨迹点，标记状态\n",
    "    for i in range(1, len(ais_data)):\n",
    "        current_status = ais_data.loc[i, 'Status']\n",
    "        \n",
    "        if current_status is not None:\n",
    "            continue  # 如果状态已经标记过，跳过计算\n",
    "        \n",
    "        archorage_status = get_archorage_status(ais_data.iloc[i], archorage_dict,current_archorage)\n",
    "        \n",
    "        ais_data.at[i, 'Status'] = 'Within Archorage' if archorage_status is not None else 'Outside Archorage'\n",
    "        if ais_data.at[i, 'Status'] == 'Within Archorage' and ais_data.at[i, 'speed'] < 0.2:\n",
    "            ais_data.at[i, 'Status'] = 'Stop'\n",
    "    return ais_data\n",
    "\n",
    "# 调用函数进行处理\n",
    "processed_ais_data = process_ais_trajectory(single_vessel_data, total_stopped_ships_per_cell)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import s2sphere\n",
    "\n",
    "# 假设 ais_data 是包含船舶轨迹数据的 DataFrame，包括 'Latitude', 'Longitude', 'Timestamp' 等列\n",
    "# archorage_dict 是包含 archorage 信息的字典，key 是 s2_cell_id，value 是 cluster\n",
    "\n",
    "\n",
    "from shapely.geometry import Polygon, Point\n",
    "def precompute_polygons(archorage_dict):\n",
    "    precomputed_polygons = {}\n",
    "    for s2_cell_id, value in archorage_dict.items():\n",
    "        cell_id = s2sphere.CellId.from_token(s2_cell_id)\n",
    "        cell = s2sphere.Cell(cell_id)\n",
    "        vertices = [s2sphere.LatLng.from_point(cell.get_vertex(i)) for i in range(4)]\n",
    "        polygon = Polygon([(v.lng().degrees, v.lat().degrees) for v in vertices])\n",
    "        precomputed_polygons[value['cluster_label']] = polygon\n",
    "    return precomputed_polygons\n",
    "\n",
    "def get_archorage_status_vectorized(row, precomputed_polygons, current_archorage):\n",
    "    lat = row['latitude']\n",
    "    lon = row['longitude']\n",
    "\n",
    "    if current_archorage is not None:\n",
    "        return current_archorage\n",
    "\n",
    "    point = Point(lon, lat)\n",
    "    for cluster, polygon in precomputed_polygons.items():\n",
    "        if point.within(polygon):\n",
    "            return cluster\n",
    "    return None\n",
    "\n",
    "def process_ais_trajectory_vectorized(ais_data, precomputed_polygons):\n",
    "    ais_data['Status'] = None\n",
    "    current_archorage = None\n",
    "\n",
    "    def process_row(row):\n",
    "        nonlocal current_archorage\n",
    "        current_status = row['Status']\n",
    "\n",
    "        if current_status is not None:\n",
    "            return row\n",
    "\n",
    "        archorage_status = get_archorage_status_vectorized(row, precomputed_polygons, current_archorage)\n",
    "        row['Status'] = 'Within Archorage' if archorage_status is not None else 'Underway'\n",
    "        if row['speed'] < 0.2: row['Status'] = 'Stop'\n",
    "\n",
    "        return row\n",
    "\n",
    "    ais_data = ais_data.apply(process_row, axis=1)\n",
    "    return ais_data\n",
    "\n",
    "# Precompute polygons outside the loop\n",
    "precomputed_polygons = precompute_polygons(archorage_dict)\n",
    "\n",
    "# Call the vectorized function for processing\n",
    "processed_ais_data = process_ais_trajectory_vectorized(ais_data, precomputed_polygons)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "processed_ais_data.to_csv('processed_ais_data.csv', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set(style=\"whitegrid\")\n",
    "plt.figure(figsize=(10, 5))\n",
    "import matplotlib.pyplot as plt\n",
    "single_vessel_data = single_vessel_data[single_vessel_data['timestamp'] > '2021-03-31']\n",
    "single_vessel_data = single_vessel_data[single_vessel_data['timestamp'] < '2021-04-3']\n",
    "with_anchorage = single_vessel_data[processed_ais_data['Status'] == 'Within Archorage']\n",
    "stop = single_vessel_data[processed_ais_data['Status'] == 'Stop']\n",
    "underway = single_vessel_data[processed_ais_data['Status'] == 'Outside Archorage']\n",
    "plt.scatter(underway['timestamp'], underway['speed'], s=1, c='blue')\n",
    "plt.scatter(stop['timestamp'], stop['speed'], s=1, c='green')\n",
    "plt.scatter(with_anchorage['timestamp'], with_anchorage['speed'], s=1, c='red')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "processed_ais_data['Status'].value_counts()\n",
    "m = folium.Map(location=[single_vessel_data['latitude'].mean(), single_vessel_data['longitude'].mean()], zoom_start=10)\n",
    "# visualize the trajs using polyline\n",
    "for mmsi, group in processed_ais_data.groupby('mmsi'):\n",
    "    for status, group_status in group.groupby('Status'):\n",
    "        # draw scatter plot for each status\n",
    "        if status == 'Stop':\n",
    "            for row in group_status.itertuples():\n",
    "                folium.CircleMarker([row.latitude, row.longitude], popup=status, radius=2, color='green', fill_opacity=0.6, opacity=0.6).add_to(m)\n",
    "        elif status == 'Within Archorage':\n",
    "            for row in group_status.itertuples():\n",
    "                folium.CircleMarker([row.latitude, row.longitude], popup=status, radius=2, color='red', fill_opacity=0.6, opacity=0.6).add_to(m)\n",
    "        else:\n",
    "            for row in group_status.itertuples():\n",
    "                folium.CircleMarker([row.latitude, row.longitude], popup=status, radius=2, color='blue', fill_opacity=0.6, opacity=0.6).add_to(m)\n",
    "            \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import folium\n",
    "import s2sphere\n",
    "import utils.visualize_utils as vu\n",
    "# Sample S2 Cell IDs for visualization\n",
    "sample_s2_cell_ids = list(total_stopped_ships_per_cell.keys())\n",
    "\n",
    "def visualize_s2_cells(s2_cell_ids, m=None):\n",
    "    # Create a Folium map centered on a specific location (e.g., Beijing)\n",
    "    m = m\n",
    "\n",
    "    for s2_cell_id in s2_cell_ids:\n",
    "        cell_id = s2sphere.CellId.from_token(s2_cell_id)\n",
    "        cell = s2sphere.Cell(cell_id)\n",
    "        vertices = []\n",
    "        for i in range(4):\n",
    "            vertex = s2sphere.LatLng.from_point(cell.get_vertex(i))\n",
    "            vertices.append([vertex.lat().degrees, vertex.lng().degrees])\n",
    "\n",
    "        # Draw a polygon on the map representing the S2 cell\n",
    "        folium.Polygon(\n",
    "            locations=vertices,\n",
    "            # different cluster different color\n",
    "\n",
    "            color=vu.get_color(total_stopped_ships_per_cell[s2_cell_id]['cluster']),\n",
    "            fill=True,\n",
    "            fill_color=vu.get_color(total_stopped_ships_per_cell[s2_cell_id]['cluster']),\n",
    "            fill_opacity=0.4,\n",
    "            tooltip=s2_cell_id,\n",
    "            popup= f'Stop Count: {len(total_stopped_ships_per_cell[s2_cell_id])}'\n",
    "        ).add_to(m)\n",
    "\n",
    "    return m\n",
    "\n",
    "# Visualize the sample S2 Cell IDs\n",
    "visualize_s2_cells(sample_s2_cell_ids, m)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "single_vessel_data_speed = single_vessel_data['speed']\n",
    "print(single_vessel_data_speed.shape)\n",
    "single_vessel_data_speed_diff_1 = single_vessel_data_speed.diff(1)\n",
    "single_vessel_data_speed_diff_2 = single_vessel_data_speed.diff(2)\n",
    "# drop the data that speed difference is too large\n",
    "single_vessel_data_speed_1 = single_vessel_data_speed[single_vessel_data_speed_diff_1.abs() < 5]\n",
    "single_vessel_data_speed_2 = single_vessel_data_speed[single_vessel_data_speed_diff_2.abs() < 10]\n",
    "# calculate the union of the two data\n",
    "single_vessel_data_speed = single_vessel_data_speed_1[single_vessel_data_speed_1.index.isin(single_vessel_data_speed_2.index)]\n",
    "print(single_vessel_data_speed.shape)\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set(style=\"whitegrid\")\n",
    "plt.figure(figsize=(30, 6))\n",
    "ax = sns.lineplot(data=single_vessel_data_speed, markers=True, dashes=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import folium\n",
    "\n",
    "\n",
    "single_vessel_data_stop = single_vessel_data[(single_vessel_data['status'] == 1.0)]\n",
    "map = folium.Map(location=[single_vessel_data_stop['latitude'].mean(), single_vessel_data_stop['longitude'].mean()], zoom_start=10)\n",
    "for index, row in single_vessel_data_stop.iterrows():\n",
    "    folium.CircleMarker(location=[row['latitude'], row['longitude']], popup=row['speed']).add_to(map)\n",
    "map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.figure(figsize=(30, 10))\n",
    "underway = single_vessel_data[single_vessel_data['behavior'] == 'underway']\n",
    "enter_port = single_vessel_data[single_vessel_data['behavior'] == 'enter_port']\n",
    "exit_port = single_vessel_data[single_vessel_data['behavior'] == 'exit_port']\n",
    "stop = single_vessel_data[single_vessel_data['behavior'] == 'stop']\n",
    "plt.scatter(underway['timestamp'], underway['speed'], s=1, c='blue')\n",
    "plt.scatter(enter_port['timestamp'], enter_port['speed'], s=1, c='green')\n",
    "plt.scatter(exit_port['timestamp'], exit_port['speed'], s=1, c='red')\n",
    "plt.scatter(stop['timestamp'], stop['speed'], s=1, c='yellow')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.signal import find_peaks\n",
    "\n",
    "def detect_speed_pattern(speed_sequence, threshold=0.1):\n",
    "    # 计算速度的一阶差分\n",
    "    velocity_diff = np.diff(speed_sequence)\n",
    "\n",
    "    # 找到速度下降点\n",
    "    neg_peaks, _ = find_peaks(-velocity_diff, height=threshold)\n",
    "\n",
    "    # 找到速度上升点\n",
    "    pos_peaks, _ = find_peaks(velocity_diff, height=threshold)\n",
    "\n",
    "    # 筛选保持接近于0的速度变化区间\n",
    "    valid_ranges = []\n",
    "    for neg_peak in neg_peaks:\n",
    "        closest_pos_peak = pos_peaks[(pos_peaks > neg_peak) & (pos_peaks < neg_peak + 10)]  # 10是一个可以调整的窗口大小\n",
    "        if closest_pos_peak.size > 0:\n",
    "            valid_ranges.append((neg_peak, closest_pos_peak[0]))\n",
    "\n",
    "    return valid_ranges\n",
    "\n",
    "# 示例\n",
    "speed_sequence = np.array([0.1, 0.2, -0.05, -0.02, -0.01, 0.05, 0.08, 0.1, 0.2, 0.15, 0.12, 0.2, 0.25, 0.3])\n",
    "pattern_ranges = detect_speed_pattern(speed_sequence, threshold=0.1)\n",
    "\n",
    "print(\"Detected pattern ranges:\", pattern_ranges)\n",
    "plt.plot(speed_sequence)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import folium\n",
    "\n",
    "\n",
    "m = folium.Map(location=[df_wp['latitude'].mean(), df_wp['longitude'].mean()], zoom_start=10)\n",
    "for i in tqdm(range(0, df_wp.shape[0])):\n",
    "    folium.Marker([df_wp['latitude'].iloc[i], df_wp['longitude'].iloc[i]]).add_to(m)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "df = pd.read_csv('./JPtoKSA.csv', parse_dates=['timestamp'])\n",
    "df.rename(columns={'timestamp': 'time', 'id':'MMSI'}, inplace=True)\n",
    "df_wp = df[(df['speed'] <=0.4) | (df['status'] == 1.0)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import utils.geohash_utils as ghu\n",
    "precision_max = 6\n",
    "cfg = EasyDict(get_config('./cfg/route_extract_cfg.yaml'))\n",
    "bounds = cfg.bounds\n",
    "tree = ghu.QuadTree(bounds, precision_max)\n",
    "tree.init_subgrids(tree.root)\n",
    "tree.init_ship_dic_df(df_wp_clustered)\n",
    "# tree.init_ship_cnt()\n",
    "# tree.filter_nodes(1)\n",
    "leaf_nodes = tree.get_leaf_nodes()\n",
    "leaf_nodes = ghu.assign_label(leaf_nodes)\n",
    "filtered_grids_map = ghu.visualize_node_list(df_wp_clustered, leaf_nodes, color='red')\n",
    "filtered_grids_map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_wp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import utils.waypoint_extract_utils as weu\n",
    "cfg = EasyDict(get_config('./cfg/route_extract_cfg.yaml'))\n",
    "df_wp_clustered = df_wp.drop_duplicates(subset=['MMSI', 'time']).reset_index(drop=True)\n",
    "df_wp_points_labels = weu.clustering_hdbscan(df_wp_clustered, cfg)\n",
    "df_wp_clustered['cluster_label'] = df_wp_points_labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import folium\n",
    "tiles= 'https://wprd01.is.autonavi.com/appmaptile?x={x}&y={y}&z={z}&lang=zh_cn&size=1&scl=1&style=7'\n",
    "# visualize the node list\n",
    "m = folium.Map(location=[55, 12], zoom_start=6, tiles=tiles,attr='高德-常规图',)\n",
    "\n",
    "for label, group in df_wp_clustered.groupby('cluster_label'):\n",
    "    r = (label * 37) % 256\n",
    "    g = (label * 73) % 256\n",
    "    b = (label * 101) % 256\n",
    "    color = 'rgb(%d, %d, %d)' % (r, g, b)\n",
    "    for idx, row in group.iterrows():\n",
    "        \n",
    "        folium.CircleMarker(location=[row['latitude'], row['longitude']], radius=1, color=color).add_to(m)\n",
    "m.save('./JPtoKSA.html')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import utils.geohash_utils as ghu\n",
    "precision_max = 8\n",
    "bounds = cfg.bounds\n",
    "tree = ghu.QuadTree(bounds, precision_max)\n",
    "tree.init_subgrids(tree.root)\n",
    "tree.init_ship_dic_df(df_wp_clustered)\n",
    "# tree.init_ship_cnt()\n",
    "# tree.filter_nodes(1)\n",
    "leaf_nodes = tree.get_leaf_nodes()\n",
    "leaf_nodes = ghu.assign_label(leaf_nodes)\n",
    "filtered_grids_map = ghu.visualize_node_list(df_wp_clustered, leaf_nodes, color='red')\n",
    "filtered_grids_map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "waypoints = ghu.leaf_nodes2waypoints(leaf_nodes)\n",
    "convex_hulls = ghu.get_convex_hulls(waypoints)\n",
    "convex_hulls.to_feather(cfg.convex_hull_save_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def process_file(file_path, cfg):\n",
    "    if file_path.endswith('csv.gz'):\n",
    "        df = pd.read_csv(file_path, usecols=cfg.columns_to_extract, compression='gzip')\n",
    "        columns_to_check = ['longitude', 'latitude', 'speed', 'heading']\n",
    "        df = df.dropna(subset=columns_to_check).reset_index(drop=True)\n",
    "        df = dpu.convert_time_format(df, time_column_ori='timestamp', time_column_new='timestamp')\n",
    "        df.drop_duplicates(subset=['mmsi', 'timestamp'], inplace=True)\n",
    "        location_1_polygon = gpd.read_file(cfg.location_1)\n",
    "        location_2_polygon = gpd.read_file(cfg.location_2)\n",
    "        df_to, df_from = dpu.filter_trajs_between(df, location_1_polygon, location_2_polygon, group_column='mmsi', time_column='timestamp')\n",
    "\n",
    "        return (df_to, df_from)\n",
    "\n",
    "\n",
    "def process_file_wrapper(args):\n",
    "    return process_file(*args)\n",
    "\n",
    "\n",
    "cfg = EasyDict(get_config('./cfg/route_extract_cfg.yaml'))\n",
    "file_paths = get_all_files(cfg.raw_data_dir)\n",
    "\n",
    "location_1_name = cfg.location_1.split('/')[-1].split('.')[0]\n",
    "location_2_name = cfg.location_2.split('/')[-1].split('.')[0]\n",
    "traj_from = pd.DataFrame(columns=cfg.columns_to_extract)\n",
    "traj_to = pd.DataFrame(columns=cfg.columns_to_extract)\n",
    "\n",
    "with Pool(processes=16) as pool:\n",
    "    # Use tqdm to display progress\n",
    "    with tqdm(total=len(file_paths), desc=\"Processing files\", unit=\"file\") as pbar:\n",
    "        # Define the partial function with fixed arguments\n",
    "        partial_process_file = partial(process_file, cfg=cfg)\n",
    "        \n",
    "        # Use imap_unordered to process files in parallel and get results asynchronously\n",
    "        results = pool.imap_unordered(partial_process_file, file_paths)\n",
    "\n",
    "        # Iterate over results and update progress bar\n",
    "        for df_to, df_from in results:\n",
    "            traj_from = pd.concat([traj_from, df_from], ignore_index=True)\n",
    "            traj_to = pd.concat([traj_to, df_to], ignore_index=True)\n",
    "            pbar.update(1)\n",
    "\n",
    "traj_from.to_csv(f'{cfg.save_dir}{location_1_name}to{location_2_name}.csv.gz', index=False, compression='gzip')\n",
    "traj_to.to_csv(f'{cfg.save_dir}{location_2_name}to{location_1_name}.csv.gz', index=False, compression='gzip')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import folium\n",
    "import utils.visualize_utils as vis\n",
    "cfg = EasyDict(get_config('./cfg/route_extract_cfg.yaml'))\n",
    "location_1_name = cfg.location_1.split('/')[-1].split('.')[0]\n",
    "location_2_name = cfg.location_2.split('/')[-1].split('.')[0]\n",
    "traj_from = pd.read_csv(f'{cfg.save_dir}{location_1_name}to{location_2_name}.csv.gz', compression='gzip')\n",
    "traj_to = pd.read_csv(f'{cfg.save_dir}{location_2_name}to{location_1_name}.csv.gz', compression='gzip')\n",
    "m_to = vis.visulize_trajs(traj_to, mode='lines', heatmap=False, save=False)\n",
    "m_from = vis.visulize_trajs(traj_to, mode='lines', heatmap=False, save=False)\n",
    "location_1_polygon = gpd.read_file(cfg.location_1)\n",
    "location_2_polygon = gpd.read_file(cfg.location_2)\n",
    "folium.GeoJson(location_1_polygon).add_to(m_to).add_to(m_from)\n",
    "folium.GeoJson(location_2_polygon).add_to(m_to).add_to(m_from)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "rapids-23.10",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}