feat: 统计一公里内高干扰小区

main.py

@@ -6,6 +6,7 @@ from decimal import Decimal, InvalidOperation
 import hashlib
 import io
 import json
+import math
 import os
 from dataclasses import dataclass
 from datetime import datetime, timedelta, timezone
@@ -62,6 +63,9 @@ MIN_INTERFERENCE_BY_NETWORK_TYPE = {
     "4G": Decimal("-110"),
 }
 
+EARTH_RADIUS_KM = 6371.0088
+NEARBY_RADIUS_KM = 1.0
+
 DEFAULT_STORAGE_ID = "stg_4d9a910d72"
 DEFAULT_ROOT = "/网优日常优化数据文档/(勿删)干扰定时小时指标"
 CELL_DATA_DIRECTORIES = (
@@ -178,6 +182,7 @@ SUMMARY_HEADER = (
     "longitude",
     "latitude",
     "azimuth",
+    "nearby_count",
 )
 
 class ProcessingError(RuntimeError):
@@ -290,6 +295,7 @@ class ApiSummaryStore:
                     longitude DECIMAL(10,6) NULL,
                     latitude DECIMAL(10,6) NULL,
                     azimuth DECIMAL(6,2) NOT NULL DEFAULT 0,
+                    nearby_count INT NOT NULL DEFAULT 0,
                     PRIMARY KEY (metric_time, cgi)
                 ) ENGINE=InnoDB DEFAULT CHARSET=utf8mb4
             """,
@@ -307,6 +313,7 @@ class ApiSummaryStore:
             "longitude": "DECIMAL(10,6) NULL",
             "latitude": "DECIMAL(10,6) NULL",
             "azimuth": "DECIMAL(6,2) NOT NULL DEFAULT 0",
+            "nearby_count": "INT NOT NULL DEFAULT 0",
         }
         for column, definition in migrations.items():
             if column not in existing_columns:
@@ -337,7 +344,7 @@ class ApiSummaryStore:
                     START TRANSACTION;
                     DELETE FROM `{self.table}` WHERE metric_time = {metric_literal};
                     INSERT INTO `{self.table}`
-                        (metric_time, network_type, cgi, cell_name, interference_dbm, longitude, latitude, azimuth)
+                        (metric_time, network_type, cgi, cell_name, interference_dbm, longitude, latitude, azimuth, nearby_count)
                     VALUES
                         {values};
                     DELETE FROM `{self.table}` WHERE metric_time <> {metric_literal};
@@ -392,6 +399,7 @@ class ApiSummaryStore:
                 sql_decimal_literal(row["longitude"], "longitude", nullable=True),
                 sql_decimal_literal(row["latitude"], "latitude", nullable=True),
                 sql_decimal_literal(row["azimuth"], "azimuth"),
+                str(int(row["nearby_count"])),
             )
         ) + ")"
 
@@ -670,6 +678,7 @@ def process(
                 }
             )
 
+        populate_nearby_counts(summary_rows)
         summary_rows.sort(key=lambda item: (item["cgi"], item["cell_name"]))
         summary_name = f"interference_summary_{window}.csv"
         write_dict_csv(temp_dir / summary_name, SUMMARY_HEADER, summary_rows)
@@ -760,6 +769,7 @@ def summary_record(
         "longitude": longitude,
         "latitude": latitude,
         "azimuth": azimuth,
+        "nearby_count": "0",
     }
 
 
@@ -774,6 +784,60 @@ def passes_interference_threshold(record: dict[str, str]) -> bool:
     return interference >= MIN_INTERFERENCE_BY_NETWORK_TYPE[network_type]
 
 
+def populate_nearby_counts(rows: list[dict[str, str]]) -> None:
+    counts = [0] * len(rows)
+    spatial_index: dict[tuple[int, int, int], list[tuple[int, float, float]]] = {}
+    for index, row in enumerate(rows):
+        row["nearby_count"] = "0"
+        if not row["longitude"] or not row["latitude"]:
+            continue
+        try:
+            longitude = float(row["longitude"])
+            latitude = float(row["latitude"])
+        except ValueError as exc:
+            raise ProcessingError(f"Invalid coordinates for CGI {row['cgi']}") from exc
+        if not math.isfinite(longitude) or not math.isfinite(latitude) or not -180 <= longitude <= 180 or not -90 <= latitude <= 90:
+            raise ProcessingError(f"Invalid coordinates for CGI {row['cgi']}")
+
+        bucket = spatial_bucket(longitude, latitude)
+        for x_offset in (-1, 0, 1):
+            for y_offset in (-1, 0, 1):
+                for z_offset in (-1, 0, 1):
+                    nearby_bucket = (bucket[0] + x_offset, bucket[1] + y_offset, bucket[2] + z_offset)
+                    for other_index, other_longitude, other_latitude in spatial_index.get(nearby_bucket, ()):
+                        if haversine_distance_km(longitude, latitude, other_longitude, other_latitude) <= NEARBY_RADIUS_KM:
+                            counts[index] += 1
+                            counts[other_index] += 1
+        spatial_index.setdefault(bucket, []).append((index, longitude, latitude))
+
+    for index, count in enumerate(counts):
+        rows[index]["nearby_count"] = str(count)
+
+
+def spatial_bucket(longitude: float, latitude: float) -> tuple[int, int, int]:
+    longitude = math.radians(longitude)
+    latitude = math.radians(latitude)
+    radius_at_latitude = EARTH_RADIUS_KM * math.cos(latitude)
+    return (
+        math.floor(radius_at_latitude * math.cos(longitude) / NEARBY_RADIUS_KM),
+        math.floor(radius_at_latitude * math.sin(longitude) / NEARBY_RADIUS_KM),
+        math.floor(EARTH_RADIUS_KM * math.sin(latitude) / NEARBY_RADIUS_KM),
+    )
+
+
+def haversine_distance_km(longitude1: float, latitude1: float, longitude2: float, latitude2: float) -> float:
+    longitude1, latitude1, longitude2, latitude2 = map(
+        math.radians,
+        (longitude1, latitude1, longitude2, latitude2),
+    )
+    longitude_delta = longitude2 - longitude1
+    latitude_delta = latitude2 - latitude1
+    haversine = math.sin(latitude_delta / 2) ** 2 + (
+        math.cos(latitude1) * math.cos(latitude2) * math.sin(longitude_delta / 2) ** 2
+    )
+    return 2 * EARTH_RADIUS_KM * math.asin(math.sqrt(min(1.0, haversine)))
+
+
 def build_cgi(prefix: str, record: dict[str, object], columns: tuple[str, ...], path: str, row: int) -> str:
     parts = [required(record, column, path, row) for column in columns]
     return "-".join(([prefix] if prefix else []) + parts)