chengkai3
2131 lines
78 KiB
Python
2131 lines
78 KiB
Python
from __future__ import annotations
|
|
|
|
import asyncio
|
|
import math
|
|
import re
|
|
from dataclasses import dataclass
|
|
from datetime import datetime, timedelta
|
|
from typing import Any
|
|
from uuid import uuid4
|
|
|
|
from fastapi import HTTPException, UploadFile, status
|
|
from sqlalchemy import and_, case, func, insert, or_, select
|
|
from sqlalchemy.orm import Session
|
|
|
|
from ..models.base import utcnow
|
|
from ..models.lightning_event import LightningCurrentEvent
|
|
from ..models.lightning_sample import LightningCurrentSample
|
|
from ..models.line_tower import LineTower
|
|
from ..schemas.lightning import (
|
|
LightningCurrentEventListResponse,
|
|
LightningCurrentEventSummary,
|
|
LightningCurrentEventUpdateRequest,
|
|
LightningDistributionEventBrief,
|
|
LightningDistributionGridCell,
|
|
LightningDistributionImportResponse,
|
|
LightningDistributionScatterPoint,
|
|
LightningDistributionStatsResponse,
|
|
LightningDistributionSummary,
|
|
LightningDistributionReportResponse,
|
|
LightningCurrentExceedancePoint,
|
|
LightningCurrentExceedanceResponse,
|
|
LightningCurrentImportResponse,
|
|
LightningPolarityStats,
|
|
LightningCurrentSampleItem,
|
|
LightningCurrentSampleListResponse,
|
|
LightningPolarity,
|
|
LightningSourceStats,
|
|
LightningSyntheticCompareResponse,
|
|
LightningSyntheticDatasetStats,
|
|
LightningTowerBufferEventItem,
|
|
LightningTowerBufferStatsResponse,
|
|
LightningTowerTerrainComputeRequest,
|
|
LightningTowerTerrainComputeResponse,
|
|
LightningTowerTerrainMetrics,
|
|
)
|
|
from .push_service import publish_topic
|
|
|
|
LIGHTNING_TOPIC = "admin.lightning-currents"
|
|
TEXT_ENCODINGS = ("utf-8-sig", "utf-8", "gbk", "latin-1")
|
|
MAX_SAMPLES = 2_000_000
|
|
INSERT_CHUNK_SIZE = 5_000
|
|
MAX_DISTRIBUTION_ROWS = 2_000_000
|
|
STANDARD_WAVE_SHAPES: tuple[tuple[str, float, float], ...] = (
|
|
("10/350", 10.0, 350.0),
|
|
("8/20", 8.0, 20.0),
|
|
("1.2/50", 1.2, 50.0),
|
|
)
|
|
TOKEN_SPLIT_PATTERN = re.compile(r"[,\t; ]+")
|
|
DEGREE_TO_KM = 111.32
|
|
TERRAIN_ALGORITHM_VERSION = "horn_3x3.v1"
|
|
|
|
|
|
@dataclass
|
|
class ParsedSeries:
|
|
currents_ka: list[float]
|
|
time_us: list[float] | None
|
|
warnings: list[str]
|
|
|
|
|
|
@dataclass
|
|
class WaveFeatures:
|
|
peak_current_ka: float | None
|
|
peak_abs_current_ka: float | None
|
|
wavefront_time_t1_us: float | None
|
|
half_value_time_t2_us: float | None
|
|
steepness_ka_per_us: float | None
|
|
action_integral_j_ohm: float | None
|
|
wave_shape: str | None
|
|
polarity: LightningPolarity
|
|
stroke_count: int
|
|
stroke_peaks_json: list[dict[str, Any]]
|
|
|
|
|
|
def serialize_lightning_event(event: LightningCurrentEvent) -> LightningCurrentEventSummary:
|
|
return LightningCurrentEventSummary(
|
|
id=event.id,
|
|
event_id=event.event_id,
|
|
source_file_name=event.source_file_name,
|
|
event_time=event.event_time,
|
|
sample_count=event.sample_count,
|
|
sample_interval_us=event.sample_interval_us,
|
|
sampling_frequency_hz=event.sampling_frequency_hz,
|
|
peak_current_ka=event.peak_current_ka,
|
|
peak_abs_current_ka=event.peak_abs_current_ka,
|
|
wavefront_time_t1_us=event.wavefront_time_t1_us,
|
|
half_value_time_t2_us=event.half_value_time_t2_us,
|
|
steepness_ka_per_us=event.steepness_ka_per_us,
|
|
action_integral_j_ohm=event.action_integral_j_ohm,
|
|
wave_shape=event.wave_shape,
|
|
polarity=event.polarity, # type: ignore[arg-type]
|
|
stroke_count=event.stroke_count,
|
|
stroke_peaks_json=event.stroke_peaks_json or [],
|
|
region_id=event.region_id,
|
|
location_tag=event.location_tag,
|
|
city=event.city,
|
|
longitude=event.longitude,
|
|
latitude=event.latitude,
|
|
altitude_m=event.altitude_m,
|
|
sensor_model=event.sensor_model,
|
|
install_position=event.install_position,
|
|
weather_level=event.weather_level,
|
|
pressure_hpa=event.pressure_hpa,
|
|
humidity_percent=event.humidity_percent,
|
|
is_synthetic=event.is_synthetic,
|
|
feature_json=event.feature_json or {},
|
|
notes=event.notes,
|
|
create_date=event.create_date,
|
|
create_user=event.create_user,
|
|
update_date=event.update_date,
|
|
update_user=event.update_user,
|
|
)
|
|
|
|
|
|
def serialize_lightning_sample(sample: LightningCurrentSample) -> LightningCurrentSampleItem:
|
|
return LightningCurrentSampleItem(
|
|
id=sample.id,
|
|
event_ref_id=sample.event_ref_id,
|
|
seq_no=sample.seq_no,
|
|
time_us=sample.time_us,
|
|
current_ka=sample.current_ka,
|
|
)
|
|
|
|
|
|
def get_lightning_event_by_id(db: Session, event_id: str) -> LightningCurrentEvent | None:
|
|
return db.execute(select(LightningCurrentEvent).where(LightningCurrentEvent.id == event_id)).scalar_one_or_none()
|
|
|
|
|
|
def get_lightning_event_by_event_id(db: Session, event_id: str) -> LightningCurrentEvent | None:
|
|
normalized = _normalize_str(event_id)
|
|
if normalized is None:
|
|
return None
|
|
return db.execute(
|
|
select(LightningCurrentEvent).where(func.lower(LightningCurrentEvent.event_id) == normalized.lower())
|
|
).scalar_one_or_none()
|
|
|
|
|
|
def list_lightning_events(
|
|
db: Session,
|
|
*,
|
|
keyword: str | None,
|
|
region_id: str | None,
|
|
polarity: str | None,
|
|
wave_shape: str | None,
|
|
is_synthetic: bool | None,
|
|
limit: int,
|
|
offset: int,
|
|
) -> LightningCurrentEventListResponse:
|
|
stmt = select(LightningCurrentEvent)
|
|
total_stmt = select(func.count()).select_from(LightningCurrentEvent)
|
|
|
|
normalized_keyword = (keyword or "").strip()
|
|
if normalized_keyword:
|
|
like = f"%{normalized_keyword}%"
|
|
predicate = or_(
|
|
LightningCurrentEvent.event_id.ilike(like),
|
|
LightningCurrentEvent.location_tag.ilike(like),
|
|
LightningCurrentEvent.city.ilike(like),
|
|
LightningCurrentEvent.source_file_name.ilike(like),
|
|
)
|
|
stmt = stmt.where(predicate)
|
|
total_stmt = total_stmt.where(predicate)
|
|
|
|
normalized_region = _normalize_str(region_id)
|
|
if normalized_region:
|
|
stmt = stmt.where(LightningCurrentEvent.region_id == normalized_region)
|
|
total_stmt = total_stmt.where(LightningCurrentEvent.region_id == normalized_region)
|
|
|
|
normalized_polarity = _normalize_str(polarity)
|
|
if normalized_polarity:
|
|
stmt = stmt.where(LightningCurrentEvent.polarity == normalized_polarity.lower())
|
|
total_stmt = total_stmt.where(LightningCurrentEvent.polarity == normalized_polarity.lower())
|
|
|
|
normalized_wave_shape = _normalize_str(wave_shape)
|
|
if normalized_wave_shape:
|
|
stmt = stmt.where(func.lower(LightningCurrentEvent.wave_shape) == normalized_wave_shape.lower())
|
|
total_stmt = total_stmt.where(func.lower(LightningCurrentEvent.wave_shape) == normalized_wave_shape.lower())
|
|
|
|
if is_synthetic is not None:
|
|
stmt = stmt.where(LightningCurrentEvent.is_synthetic == is_synthetic)
|
|
total_stmt = total_stmt.where(LightningCurrentEvent.is_synthetic == is_synthetic)
|
|
|
|
total = int(db.scalar(total_stmt) or 0)
|
|
rows = db.execute(
|
|
stmt.order_by(
|
|
LightningCurrentEvent.event_time.desc(),
|
|
LightningCurrentEvent.update_date.desc(),
|
|
LightningCurrentEvent.id.desc(),
|
|
)
|
|
.offset(offset)
|
|
.limit(limit)
|
|
).scalars().all()
|
|
|
|
return LightningCurrentEventListResponse(
|
|
items=[serialize_lightning_event(item) for item in rows],
|
|
total=total,
|
|
limit=limit,
|
|
offset=offset,
|
|
)
|
|
|
|
|
|
def update_lightning_event(
|
|
db: Session,
|
|
event_id: str,
|
|
payload: LightningCurrentEventUpdateRequest,
|
|
*,
|
|
actor_user_id: str,
|
|
) -> LightningCurrentEventSummary | None:
|
|
event = get_lightning_event_by_id(db, event_id)
|
|
if not event:
|
|
return None
|
|
|
|
data = payload.model_dump(exclude_unset=True)
|
|
|
|
if "event_time" in data:
|
|
event.event_time = data["event_time"]
|
|
if "region_id" in data:
|
|
event.region_id = _normalize_str(data["region_id"])
|
|
if "location_tag" in data:
|
|
event.location_tag = _normalize_str(data["location_tag"])
|
|
if "city" in data:
|
|
event.city = _normalize_str(data["city"])
|
|
if "longitude" in data:
|
|
event.longitude = data["longitude"]
|
|
if "latitude" in data:
|
|
event.latitude = data["latitude"]
|
|
if "altitude_m" in data:
|
|
event.altitude_m = data["altitude_m"]
|
|
if "sensor_model" in data:
|
|
event.sensor_model = _normalize_str(data["sensor_model"])
|
|
if "install_position" in data:
|
|
event.install_position = _normalize_str(data["install_position"])
|
|
if "weather_level" in data:
|
|
event.weather_level = _normalize_str(data["weather_level"])
|
|
if "pressure_hpa" in data:
|
|
event.pressure_hpa = data["pressure_hpa"]
|
|
if "humidity_percent" in data:
|
|
event.humidity_percent = data["humidity_percent"]
|
|
if "is_synthetic" in data and data["is_synthetic"] is not None:
|
|
event.is_synthetic = bool(data["is_synthetic"])
|
|
if "wave_shape" in data:
|
|
event.wave_shape = _normalize_str(data["wave_shape"])
|
|
if "notes" in data:
|
|
event.notes = _normalize_str(data["notes"])
|
|
if "feature_json" in data and data["feature_json"] is not None:
|
|
event.feature_json = dict(data["feature_json"])
|
|
|
|
event.update_user = actor_user_id
|
|
event.update_date = utcnow()
|
|
db.commit()
|
|
|
|
saved = get_lightning_event_by_id(db, event_id)
|
|
if not saved:
|
|
return None
|
|
|
|
_publish_lightning_change("lightning.updated", {"action": "updated", "event_id": saved.id})
|
|
return serialize_lightning_event(saved)
|
|
|
|
|
|
def delete_lightning_event(db: Session, event_id: str) -> bool:
|
|
event = get_lightning_event_by_id(db, event_id)
|
|
if not event:
|
|
return False
|
|
|
|
db.delete(event)
|
|
db.commit()
|
|
_publish_lightning_change("lightning.deleted", {"action": "deleted", "event_id": event_id})
|
|
return True
|
|
|
|
|
|
def list_lightning_samples(
|
|
db: Session,
|
|
*,
|
|
event_id: str,
|
|
limit: int,
|
|
offset: int,
|
|
) -> LightningCurrentSampleListResponse:
|
|
total = int(
|
|
db.scalar(
|
|
select(func.count()).select_from(LightningCurrentSample).where(LightningCurrentSample.event_ref_id == event_id)
|
|
)
|
|
or 0
|
|
)
|
|
rows = db.execute(
|
|
select(LightningCurrentSample)
|
|
.where(LightningCurrentSample.event_ref_id == event_id)
|
|
.order_by(LightningCurrentSample.seq_no.asc())
|
|
.offset(offset)
|
|
.limit(limit)
|
|
).scalars().all()
|
|
|
|
return LightningCurrentSampleListResponse(
|
|
items=[serialize_lightning_sample(item) for item in rows],
|
|
total=total,
|
|
limit=limit,
|
|
offset=offset,
|
|
)
|
|
|
|
|
|
def import_lightning_event_from_file(
|
|
db: Session,
|
|
*,
|
|
file: UploadFile,
|
|
actor_user_id: str,
|
|
event_id: str | None,
|
|
event_time: datetime | None,
|
|
sample_interval_us: float | None,
|
|
region_id: str | None,
|
|
location_tag: str | None,
|
|
city: str | None,
|
|
longitude: float | None,
|
|
latitude: float | None,
|
|
altitude_m: float | None,
|
|
sensor_model: str | None,
|
|
install_position: str | None,
|
|
weather_level: str | None,
|
|
pressure_hpa: float | None,
|
|
humidity_percent: float | None,
|
|
is_synthetic: bool,
|
|
notes: str | None,
|
|
) -> LightningCurrentImportResponse:
|
|
content = file.file.read()
|
|
if not content:
|
|
raise HTTPException(status_code=status.HTTP_400_BAD_REQUEST, detail="上传文件为空")
|
|
|
|
text = _decode_text_bytes(content)
|
|
parsed = _parse_current_series(text)
|
|
if len(parsed.currents_ka) < 2:
|
|
raise HTTPException(status_code=status.HTTP_400_BAD_REQUEST, detail="有效采样点数量不足（至少需要 2 个点）")
|
|
if len(parsed.currents_ka) > MAX_SAMPLES:
|
|
raise HTTPException(
|
|
status_code=status.HTTP_400_BAD_REQUEST,
|
|
detail=f"采样点过多（>{MAX_SAMPLES}），请拆分后再导入",
|
|
)
|
|
|
|
time_axis_us, inferred_interval_us, inferred_fs_hz, time_warnings = _build_time_axis(
|
|
parsed.time_us,
|
|
sample_interval_us=sample_interval_us,
|
|
sample_count=len(parsed.currents_ka),
|
|
)
|
|
warnings = [*parsed.warnings, *time_warnings]
|
|
features = _extract_wave_features(parsed.currents_ka, time_axis_us, inferred_interval_us)
|
|
|
|
normalized_event_id = _normalize_str(event_id) or _generate_event_id()
|
|
if get_lightning_event_by_event_id(db, normalized_event_id):
|
|
raise HTTPException(status_code=status.HTTP_409_CONFLICT, detail=f"事件编号已存在：{normalized_event_id}")
|
|
|
|
now = utcnow()
|
|
event = LightningCurrentEvent(
|
|
event_id=normalized_event_id,
|
|
source_file_name=_normalize_str(file.filename),
|
|
event_time=event_time,
|
|
sample_count=len(parsed.currents_ka),
|
|
sample_interval_us=inferred_interval_us,
|
|
sampling_frequency_hz=inferred_fs_hz,
|
|
peak_current_ka=features.peak_current_ka,
|
|
peak_abs_current_ka=features.peak_abs_current_ka,
|
|
wavefront_time_t1_us=features.wavefront_time_t1_us,
|
|
half_value_time_t2_us=features.half_value_time_t2_us,
|
|
steepness_ka_per_us=features.steepness_ka_per_us,
|
|
action_integral_j_ohm=features.action_integral_j_ohm,
|
|
wave_shape=features.wave_shape,
|
|
polarity=features.polarity,
|
|
stroke_count=features.stroke_count,
|
|
stroke_peaks_json=features.stroke_peaks_json,
|
|
region_id=_normalize_str(region_id),
|
|
location_tag=_normalize_str(location_tag),
|
|
city=_normalize_str(city),
|
|
longitude=longitude,
|
|
latitude=latitude,
|
|
altitude_m=altitude_m,
|
|
sensor_model=_normalize_str(sensor_model),
|
|
install_position=_normalize_str(install_position),
|
|
weather_level=_normalize_str(weather_level),
|
|
pressure_hpa=pressure_hpa,
|
|
humidity_percent=humidity_percent,
|
|
is_synthetic=bool(is_synthetic),
|
|
notes=_normalize_str(notes),
|
|
feature_json={
|
|
"source": "import-file",
|
|
"has_explicit_time_axis": parsed.time_us is not None,
|
|
},
|
|
create_user=actor_user_id,
|
|
update_user=actor_user_id,
|
|
create_date=now,
|
|
update_date=now,
|
|
)
|
|
db.add(event)
|
|
db.flush()
|
|
|
|
sample_insert_stmt = insert(LightningCurrentSample)
|
|
total_samples = len(parsed.currents_ka)
|
|
for start in range(0, total_samples, INSERT_CHUNK_SIZE):
|
|
end = min(start + INSERT_CHUNK_SIZE, total_samples)
|
|
chunk = [
|
|
{
|
|
"event_ref_id": event.id,
|
|
"seq_no": idx + 1,
|
|
"time_us": time_axis_us[idx],
|
|
"current_ka": parsed.currents_ka[idx],
|
|
}
|
|
for idx in range(start, end)
|
|
]
|
|
db.execute(sample_insert_stmt, chunk)
|
|
|
|
db.commit()
|
|
saved = get_lightning_event_by_id(db, event.id)
|
|
if not saved:
|
|
raise HTTPException(status_code=status.HTTP_500_INTERNAL_SERVER_ERROR, detail="雷电流事件保存失败")
|
|
|
|
_publish_lightning_change("lightning.imported", {"action": "imported", "event_id": saved.id})
|
|
return LightningCurrentImportResponse(
|
|
event=serialize_lightning_event(saved),
|
|
warning_count=len(warnings),
|
|
warnings=warnings,
|
|
)
|
|
|
|
|
|
def import_lightning_distribution_from_file(
|
|
db: Session,
|
|
*,
|
|
file: UploadFile,
|
|
actor_user_id: str,
|
|
event_year: int | None,
|
|
region_id: str | None,
|
|
location_tag: str | None,
|
|
city: str | None,
|
|
is_synthetic: bool,
|
|
notes: str | None,
|
|
) -> LightningDistributionImportResponse:
|
|
content = file.file.read()
|
|
if not content:
|
|
raise HTTPException(status_code=status.HTTP_400_BAD_REQUEST, detail="上传文件为空")
|
|
|
|
text = _decode_text_bytes(content)
|
|
normalized_region = _normalize_str(region_id)
|
|
normalized_location = _normalize_str(location_tag)
|
|
normalized_city = _normalize_str(city)
|
|
normalized_notes = _normalize_str(notes)
|
|
source_file_name = _normalize_str(file.filename)
|
|
parsed_event_time = datetime(event_year, 1, 1) if event_year is not None else None
|
|
|
|
rows_to_insert: list[dict[str, Any]] = []
|
|
inserted_count = 0
|
|
skipped_count = 0
|
|
warnings: list[str] = []
|
|
now = utcnow()
|
|
|
|
insert_stmt = insert(LightningCurrentEvent)
|
|
for line_no, raw_line in enumerate(text.splitlines(), start=1):
|
|
if inserted_count + skipped_count >= MAX_DISTRIBUTION_ROWS:
|
|
raise HTTPException(
|
|
status_code=status.HTTP_400_BAD_REQUEST,
|
|
detail=f"记录数超过上限（>{MAX_DISTRIBUTION_ROWS}），请拆分文件后重试",
|
|
)
|
|
|
|
parsed = _parse_distribution_line(raw_line)
|
|
if not parsed:
|
|
stripped = raw_line.strip()
|
|
if stripped and len(warnings) < 50:
|
|
warnings.append(f"第 {line_no} 行无法解析为“纬度 经度 电流”，已跳过")
|
|
skipped_count += 1
|
|
continue
|
|
|
|
latitude, longitude, current_ka = parsed
|
|
if not (-90.0 <= latitude <= 90.0 and -180.0 <= longitude <= 180.0):
|
|
skipped_count += 1
|
|
if len(warnings) < 50:
|
|
warnings.append(f"第 {line_no} 行坐标越界（lat={latitude}, lon={longitude}），已跳过")
|
|
continue
|
|
|
|
if math.isclose(current_ka, 0.0, abs_tol=1e-12):
|
|
polarity: LightningPolarity = "unknown"
|
|
else:
|
|
polarity = "positive" if current_ka > 0 else "negative"
|
|
|
|
event_id = _generate_distribution_event_id()
|
|
stroke_peaks = [
|
|
{
|
|
"seq_no": 1,
|
|
"sample_index": 1,
|
|
"time_us": 0.0,
|
|
"current_ka": current_ka,
|
|
}
|
|
]
|
|
rows_to_insert.append(
|
|
{
|
|
"id": uuid4().hex,
|
|
"event_id": event_id,
|
|
"source_file_name": source_file_name,
|
|
"event_time": parsed_event_time,
|
|
"sample_count": 1,
|
|
"sample_interval_us": None,
|
|
"sampling_frequency_hz": None,
|
|
"peak_current_ka": current_ka,
|
|
"peak_abs_current_ka": abs(current_ka),
|
|
"wavefront_time_t1_us": None,
|
|
"half_value_time_t2_us": None,
|
|
"steepness_ka_per_us": None,
|
|
"action_integral_j_ohm": None,
|
|
"wave_shape": None,
|
|
"polarity": polarity,
|
|
"stroke_count": 1,
|
|
"stroke_peaks_json": stroke_peaks,
|
|
"region_id": normalized_region,
|
|
"location_tag": normalized_location,
|
|
"city": normalized_city,
|
|
"longitude": longitude,
|
|
"latitude": latitude,
|
|
"altitude_m": None,
|
|
"sensor_model": None,
|
|
"install_position": None,
|
|
"weather_level": None,
|
|
"pressure_hpa": None,
|
|
"humidity_percent": None,
|
|
"is_synthetic": bool(is_synthetic),
|
|
"feature_json": {
|
|
"source": "distribution-import",
|
|
"line_no": line_no,
|
|
},
|
|
"notes": normalized_notes,
|
|
"create_user": actor_user_id,
|
|
"update_user": actor_user_id,
|
|
"create_date": now,
|
|
"update_date": now,
|
|
}
|
|
)
|
|
if len(rows_to_insert) >= INSERT_CHUNK_SIZE:
|
|
db.execute(insert_stmt, rows_to_insert)
|
|
inserted_count += len(rows_to_insert)
|
|
rows_to_insert.clear()
|
|
|
|
if rows_to_insert:
|
|
db.execute(insert_stmt, rows_to_insert)
|
|
inserted_count += len(rows_to_insert)
|
|
|
|
if inserted_count == 0:
|
|
raise HTTPException(status_code=status.HTTP_400_BAD_REQUEST, detail="未解析到可导入的有效雷电分布记录")
|
|
|
|
db.commit()
|
|
_publish_lightning_change(
|
|
"lightning.distribution.imported",
|
|
{"action": "distribution_imported", "imported_count": inserted_count},
|
|
)
|
|
return LightningDistributionImportResponse(
|
|
imported_count=inserted_count,
|
|
skipped_count=skipped_count,
|
|
warning_count=len(warnings),
|
|
warnings=warnings,
|
|
)
|
|
|
|
|
|
def get_lightning_distribution_stats(
|
|
db: Session,
|
|
*,
|
|
min_lat: float | None,
|
|
max_lat: float | None,
|
|
min_lon: float | None,
|
|
max_lon: float | None,
|
|
region_id: str | None,
|
|
city: str | None,
|
|
location_tag: str | None,
|
|
polarity: str | None,
|
|
is_synthetic: bool | None,
|
|
grid_size_km: float,
|
|
years: float | None,
|
|
grid_limit: int,
|
|
scatter_limit: int,
|
|
thresholds_ka: list[float] | None,
|
|
) -> LightningDistributionStatsResponse:
|
|
filters = _build_distribution_filters(
|
|
min_lat=min_lat,
|
|
max_lat=max_lat,
|
|
min_lon=min_lon,
|
|
max_lon=max_lon,
|
|
region_id=region_id,
|
|
city=city,
|
|
location_tag=location_tag,
|
|
polarity=polarity,
|
|
is_synthetic=is_synthetic,
|
|
)
|
|
|
|
summary_row = db.execute(
|
|
select(
|
|
func.count().label("total_records"),
|
|
func.min(LightningCurrentEvent.latitude).label("min_lat"),
|
|
func.max(LightningCurrentEvent.latitude).label("max_lat"),
|
|
func.min(LightningCurrentEvent.longitude).label("min_lon"),
|
|
func.max(LightningCurrentEvent.longitude).label("max_lon"),
|
|
func.max(LightningCurrentEvent.peak_abs_current_ka).label("max_abs"),
|
|
func.avg(LightningCurrentEvent.peak_abs_current_ka).label("avg_abs"),
|
|
func.sum(case((LightningCurrentEvent.polarity == "positive", 1), else_=0)).label("positive_count"),
|
|
func.sum(case((LightningCurrentEvent.polarity == "negative", 1), else_=0)).label("negative_count"),
|
|
func.sum(case((LightningCurrentEvent.polarity == "mixed", 1), else_=0)).label("mixed_count"),
|
|
func.sum(case((LightningCurrentEvent.polarity == "unknown", 1), else_=0)).label("unknown_count"),
|
|
func.sum(case((LightningCurrentEvent.is_synthetic.is_(True), 1), else_=0)).label("synthetic_count"),
|
|
).where(*filters)
|
|
).one()
|
|
|
|
total_records = int(summary_row.total_records or 0)
|
|
if total_records == 0:
|
|
return LightningDistributionStatsResponse(
|
|
summary=LightningDistributionSummary(
|
|
total_records=0,
|
|
area_km2=0.0,
|
|
data_years=years or 1.0,
|
|
grid_size_km=grid_size_km,
|
|
overall_ng_per_km2_year=0.0,
|
|
max_abs_current_ka=None,
|
|
avg_abs_current_ka=None,
|
|
),
|
|
polarity=LightningPolarityStats(),
|
|
sources=LightningSourceStats(),
|
|
grid_cells=[],
|
|
scatter_points=[],
|
|
p_curve=[],
|
|
)
|
|
|
|
min_lat_value = float(summary_row.min_lat)
|
|
max_lat_value = float(summary_row.max_lat)
|
|
min_lon_value = float(summary_row.min_lon)
|
|
max_lon_value = float(summary_row.max_lon)
|
|
center_lat = (min_lat_value + max_lat_value) / 2.0
|
|
km_per_lon = _safe_km_per_lon(center_lat)
|
|
|
|
width_km = max((max_lon_value - min_lon_value) * km_per_lon, grid_size_km)
|
|
height_km = max((max_lat_value - min_lat_value) * DEGREE_TO_KM, grid_size_km)
|
|
area_km2 = max(width_km * height_km, grid_size_km * grid_size_km)
|
|
data_years = years if years is not None else _infer_data_years(db, filters)
|
|
data_years = max(data_years, 1e-6)
|
|
overall_ng = total_records / (area_km2 * data_years)
|
|
|
|
positive_count = int(summary_row.positive_count or 0)
|
|
negative_count = int(summary_row.negative_count or 0)
|
|
mixed_count = int(summary_row.mixed_count or 0)
|
|
unknown_count = int(summary_row.unknown_count or 0)
|
|
synthetic_count = int(summary_row.synthetic_count or 0)
|
|
|
|
cell_x_expr = func.floor(((LightningCurrentEvent.longitude - min_lon_value) * km_per_lon) / grid_size_km)
|
|
cell_y_expr = func.floor(((LightningCurrentEvent.latitude - min_lat_value) * DEGREE_TO_KM) / grid_size_km)
|
|
grid_rows = db.execute(
|
|
select(
|
|
cell_x_expr.label("cell_x"),
|
|
cell_y_expr.label("cell_y"),
|
|
func.count().label("strike_count"),
|
|
func.max(LightningCurrentEvent.peak_abs_current_ka).label("i_max"),
|
|
func.avg(LightningCurrentEvent.peak_abs_current_ka).label("i_avg"),
|
|
func.sum(case((LightningCurrentEvent.polarity == "positive", 1), else_=0)).label("positive_count"),
|
|
)
|
|
.where(*filters)
|
|
.group_by(cell_x_expr, cell_y_expr)
|
|
.order_by(func.count().desc())
|
|
.limit(grid_limit)
|
|
).all()
|
|
|
|
grid_cells: list[LightningDistributionGridCell] = []
|
|
for row in grid_rows:
|
|
grid_x = int(float(row.cell_x))
|
|
grid_y = int(float(row.cell_y))
|
|
x0_km = grid_x * grid_size_km
|
|
x1_km = (grid_x + 1) * grid_size_km
|
|
y0_km = grid_y * grid_size_km
|
|
y1_km = (grid_y + 1) * grid_size_km
|
|
|
|
cell_min_lat = min_lat_value + y0_km / DEGREE_TO_KM
|
|
cell_max_lat = min_lat_value + y1_km / DEGREE_TO_KM
|
|
cell_center_lat = (cell_min_lat + cell_max_lat) / 2.0
|
|
cell_km_per_lon = _safe_km_per_lon(cell_center_lat)
|
|
cell_min_lon = min_lon_value + x0_km / cell_km_per_lon
|
|
cell_max_lon = min_lon_value + x1_km / cell_km_per_lon
|
|
cell_center_lon = (cell_min_lon + cell_max_lon) / 2.0
|
|
|
|
strike_count = int(row.strike_count or 0)
|
|
positive_in_cell = int(row.positive_count or 0)
|
|
grid_cells.append(
|
|
LightningDistributionGridCell(
|
|
grid_x=grid_x,
|
|
grid_y=grid_y,
|
|
min_lat=cell_min_lat,
|
|
max_lat=cell_max_lat,
|
|
min_lon=cell_min_lon,
|
|
max_lon=cell_max_lon,
|
|
center_lat=cell_center_lat,
|
|
center_lon=cell_center_lon,
|
|
strike_count=strike_count,
|
|
ng_per_km2_year=strike_count / (grid_size_km * grid_size_km * data_years),
|
|
i_max_ka=float(row.i_max) if row.i_max is not None else None,
|
|
i_avg_ka=float(row.i_avg) if row.i_avg is not None else None,
|
|
positive_ratio=(positive_in_cell / strike_count) if strike_count > 0 else 0.0,
|
|
)
|
|
)
|
|
|
|
scatter_rows = db.execute(
|
|
select(
|
|
LightningCurrentEvent.id,
|
|
LightningCurrentEvent.event_id,
|
|
LightningCurrentEvent.longitude,
|
|
LightningCurrentEvent.latitude,
|
|
LightningCurrentEvent.peak_current_ka,
|
|
LightningCurrentEvent.peak_abs_current_ka,
|
|
LightningCurrentEvent.polarity,
|
|
LightningCurrentEvent.region_id,
|
|
LightningCurrentEvent.city,
|
|
LightningCurrentEvent.location_tag,
|
|
LightningCurrentEvent.event_time,
|
|
)
|
|
.where(*filters)
|
|
.order_by(
|
|
LightningCurrentEvent.peak_abs_current_ka.desc(),
|
|
LightningCurrentEvent.update_date.desc(),
|
|
LightningCurrentEvent.id.desc(),
|
|
)
|
|
.limit(scatter_limit)
|
|
).all()
|
|
scatter_points = [
|
|
LightningDistributionScatterPoint(
|
|
id=str(row.id),
|
|
event_id=str(row.event_id),
|
|
longitude=float(row.longitude),
|
|
latitude=float(row.latitude),
|
|
current_ka=float(row.peak_current_ka) if row.peak_current_ka is not None else None,
|
|
abs_current_ka=float(row.peak_abs_current_ka) if row.peak_abs_current_ka is not None else None,
|
|
polarity=str(row.polarity), # type: ignore[arg-type]
|
|
region_id=row.region_id,
|
|
city=row.city,
|
|
location_tag=row.location_tag,
|
|
event_time=row.event_time,
|
|
)
|
|
for row in scatter_rows
|
|
]
|
|
|
|
p_curve = _build_p_curve_with_filters(
|
|
db,
|
|
filters=filters,
|
|
total_events=total_records,
|
|
max_peak=float(summary_row.max_abs) if summary_row.max_abs is not None else 0.0,
|
|
thresholds_ka=thresholds_ka,
|
|
)
|
|
|
|
return LightningDistributionStatsResponse(
|
|
summary=LightningDistributionSummary(
|
|
total_records=total_records,
|
|
area_km2=area_km2,
|
|
data_years=data_years,
|
|
grid_size_km=grid_size_km,
|
|
overall_ng_per_km2_year=overall_ng,
|
|
max_abs_current_ka=float(summary_row.max_abs) if summary_row.max_abs is not None else None,
|
|
avg_abs_current_ka=float(summary_row.avg_abs) if summary_row.avg_abs is not None else None,
|
|
),
|
|
polarity=LightningPolarityStats(
|
|
positive_count=positive_count,
|
|
negative_count=negative_count,
|
|
mixed_count=mixed_count,
|
|
unknown_count=unknown_count,
|
|
positive_ratio=(positive_count / total_records) if total_records > 0 else 0.0,
|
|
negative_ratio=(negative_count / total_records) if total_records > 0 else 0.0,
|
|
),
|
|
sources=LightningSourceStats(
|
|
measured_count=total_records - synthetic_count,
|
|
synthetic_count=synthetic_count,
|
|
),
|
|
grid_cells=grid_cells,
|
|
scatter_points=scatter_points,
|
|
p_curve=p_curve,
|
|
)
|
|
|
|
|
|
def get_tower_buffer_stats(
|
|
db: Session,
|
|
*,
|
|
tower_id: str | None,
|
|
longitude: float | None,
|
|
latitude: float | None,
|
|
radius_km: float,
|
|
design_current_ka: float,
|
|
years: float | None,
|
|
region_id: str | None,
|
|
is_synthetic: bool | None,
|
|
include_events_limit: int,
|
|
) -> LightningTowerBufferStatsResponse:
|
|
if radius_km <= 0:
|
|
raise HTTPException(status_code=status.HTTP_400_BAD_REQUEST, detail="radius_km 必须大于 0")
|
|
if design_current_ka <= 0:
|
|
raise HTTPException(status_code=status.HTTP_400_BAD_REQUEST, detail="design_current_ka 必须大于 0")
|
|
|
|
tower: LineTower | None = None
|
|
center_lon: float
|
|
center_lat: float
|
|
if tower_id:
|
|
tower = db.execute(select(LineTower).where(LineTower.id == tower_id)).scalar_one_or_none()
|
|
if not tower:
|
|
raise HTTPException(status_code=status.HTTP_404_NOT_FOUND, detail="杆塔不存在")
|
|
if tower.longitude is None or tower.latitude is None:
|
|
raise HTTPException(status_code=status.HTTP_400_BAD_REQUEST, detail="该杆塔缺少经纬度，无法执行缓冲区分析")
|
|
center_lon = float(tower.longitude)
|
|
center_lat = float(tower.latitude)
|
|
else:
|
|
if longitude is None or latitude is None:
|
|
raise HTTPException(status_code=status.HTTP_400_BAD_REQUEST, detail="请提供 tower_id 或经纬度坐标")
|
|
center_lon = float(longitude)
|
|
center_lat = float(latitude)
|
|
|
|
lat_delta = radius_km / DEGREE_TO_KM
|
|
lon_delta = radius_km / _safe_km_per_lon(center_lat)
|
|
filters = _build_distribution_filters(
|
|
min_lat=center_lat - lat_delta,
|
|
max_lat=center_lat + lat_delta,
|
|
min_lon=center_lon - lon_delta,
|
|
max_lon=center_lon + lon_delta,
|
|
region_id=region_id,
|
|
city=None,
|
|
location_tag=None,
|
|
polarity=None,
|
|
is_synthetic=is_synthetic,
|
|
)
|
|
|
|
candidate_rows = db.execute(
|
|
select(
|
|
LightningCurrentEvent.id,
|
|
LightningCurrentEvent.event_id,
|
|
LightningCurrentEvent.longitude,
|
|
LightningCurrentEvent.latitude,
|
|
LightningCurrentEvent.peak_current_ka,
|
|
LightningCurrentEvent.peak_abs_current_ka,
|
|
LightningCurrentEvent.polarity,
|
|
LightningCurrentEvent.event_time,
|
|
LightningCurrentEvent.location_tag,
|
|
LightningCurrentEvent.city,
|
|
).where(*filters)
|
|
).all()
|
|
|
|
events: list[LightningTowerBufferEventItem] = []
|
|
positive_count = 0
|
|
max_abs = 0.0
|
|
sum_abs = 0.0
|
|
exceed_count = 0
|
|
event_timestamps: list[datetime] = []
|
|
for row in candidate_rows:
|
|
if row.longitude is None or row.latitude is None:
|
|
continue
|
|
distance = _haversine_km(center_lat, center_lon, float(row.latitude), float(row.longitude))
|
|
if distance > radius_km:
|
|
continue
|
|
abs_current = float(row.peak_abs_current_ka) if row.peak_abs_current_ka is not None else 0.0
|
|
if str(row.polarity) == "positive":
|
|
positive_count += 1
|
|
if abs_current > max_abs:
|
|
max_abs = abs_current
|
|
sum_abs += abs_current
|
|
if abs_current >= design_current_ka:
|
|
exceed_count += 1
|
|
if row.event_time is not None:
|
|
event_timestamps.append(row.event_time)
|
|
|
|
events.append(
|
|
LightningTowerBufferEventItem(
|
|
id=str(row.id),
|
|
event_id=str(row.event_id),
|
|
longitude=float(row.longitude),
|
|
latitude=float(row.latitude),
|
|
current_ka=float(row.peak_current_ka) if row.peak_current_ka is not None else None,
|
|
abs_current_ka=abs_current,
|
|
polarity=str(row.polarity), # type: ignore[arg-type]
|
|
event_time=row.event_time,
|
|
location_tag=row.location_tag,
|
|
city=row.city,
|
|
distance_km=distance,
|
|
)
|
|
)
|
|
|
|
strike_count = len(events)
|
|
avg_abs = (sum_abs / strike_count) if strike_count > 0 else None
|
|
data_years = years if years is not None else _resolve_data_years_from_timestamps(event_timestamps)
|
|
data_years = max(data_years, 1e-6)
|
|
area_km2 = math.pi * (radius_km ** 2)
|
|
ng = strike_count / (area_km2 * data_years) if strike_count > 0 else 0.0
|
|
positive_ratio = (positive_count / strike_count) if strike_count > 0 else 0.0
|
|
terrain_metrics = _build_tower_terrain_metrics_from_tower(tower)
|
|
terrain_exposure = (
|
|
terrain_metrics.terrain_exposure_index
|
|
if terrain_metrics is not None and terrain_metrics.terrain_exposure_index is not None
|
|
else 0.0
|
|
)
|
|
ng_for_risk = ng * (1.0 + 0.25 * max(0.0, min(1.0, terrain_exposure)))
|
|
risk_level = _derive_tower_risk_level(
|
|
strike_count=strike_count,
|
|
exceed_design_count=exceed_count,
|
|
max_abs_current_ka=max_abs if strike_count > 0 else None,
|
|
design_current_ka=design_current_ka,
|
|
ng_per_km2_year=ng_for_risk,
|
|
)
|
|
recommended_action = _tower_risk_recommendation(risk_level)
|
|
if terrain_metrics is not None and terrain_exposure >= 0.7:
|
|
recommended_action = f"{recommended_action} 地形暴露指数较高，建议同步复核边坡防护与接地体布置。"
|
|
|
|
sorted_events = sorted(events, key=lambda item: ((item.abs_current_ka or 0.0), -item.distance_km), reverse=True)
|
|
return LightningTowerBufferStatsResponse(
|
|
tower_id=tower.id if tower else None,
|
|
tower_no=tower.tower_no if tower else None,
|
|
line_id=tower.line_id if tower else None,
|
|
center_longitude=center_lon,
|
|
center_latitude=center_lat,
|
|
radius_km=radius_km,
|
|
design_current_ka=design_current_ka,
|
|
strike_count=strike_count,
|
|
exceed_design_count=exceed_count,
|
|
max_abs_current_ka=(max_abs if strike_count > 0 else None),
|
|
avg_abs_current_ka=avg_abs,
|
|
ng_per_km2_year=ng,
|
|
positive_ratio=positive_ratio,
|
|
risk_level=risk_level,
|
|
recommended_action=recommended_action,
|
|
events=sorted_events[:include_events_limit],
|
|
terrain_metrics=terrain_metrics,
|
|
)
|
|
|
|
|
|
def compute_tower_terrain_metrics(
|
|
db: Session,
|
|
*,
|
|
payload: LightningTowerTerrainComputeRequest,
|
|
actor_user_id: str,
|
|
can_persist: bool,
|
|
) -> LightningTowerTerrainComputeResponse:
|
|
warnings: list[str] = []
|
|
tower: LineTower | None = None
|
|
|
|
if payload.tower_id:
|
|
tower = db.execute(select(LineTower).where(LineTower.id == payload.tower_id)).scalar_one_or_none()
|
|
if not tower:
|
|
raise HTTPException(status_code=status.HTTP_404_NOT_FOUND, detail="杆塔不存在")
|
|
|
|
center_lon = payload.longitude if payload.longitude is not None else (float(tower.longitude) if tower and tower.longitude is not None else None)
|
|
center_lat = payload.latitude if payload.latitude is not None else (float(tower.latitude) if tower and tower.latitude is not None else None)
|
|
if center_lon is None or center_lat is None:
|
|
raise HTTPException(status_code=status.HTTP_400_BAD_REQUEST, detail="缺少中心经纬度，请传入 tower_id 或经纬度")
|
|
|
|
z = payload.dem_grid_m
|
|
cell_size_m = float(payload.cell_size_m)
|
|
dem_resolution_m = float(payload.dem_resolution_m) if payload.dem_resolution_m is not None else cell_size_m
|
|
|
|
dz_dx, dz_dy = _compute_horn_gradient(z, cell_size_m=cell_size_m)
|
|
slope_rad = math.atan(math.sqrt(dz_dx * dz_dx + dz_dy * dz_dy))
|
|
slope_deg = math.degrees(slope_rad)
|
|
aspect_deg = _compute_aspect_deg(dz_dx, dz_dy)
|
|
relief_m_50 = max(max(row) for row in z) - min(min(row) for row in z)
|
|
|
|
neighbor_slopes = _compute_neighbor_slopes(z, cell_size_m=cell_size_m)
|
|
slope_mean_deg = (sum(neighbor_slopes) / len(neighbor_slopes)) if neighbor_slopes else None
|
|
slope_p95_deg = _percentile(neighbor_slopes, 0.95) if neighbor_slopes else None
|
|
slope_max_deg = max(neighbor_slopes) if neighbor_slopes else None
|
|
|
|
line_azimuth_deg = _infer_tower_line_azimuth_deg(db, tower) if tower else None
|
|
slope_along_line_deg: float | None = None
|
|
slope_cross_line_deg: float | None = None
|
|
if line_azimuth_deg is not None:
|
|
slope_along_line_deg = math.degrees(math.atan(_directional_gradient(dz_dx, dz_dy, line_azimuth_deg)))
|
|
slope_cross_line_deg = math.degrees(math.atan(_directional_gradient(dz_dx, dz_dy, (line_azimuth_deg + 90.0) % 360.0)))
|
|
elif tower is not None and tower.line_id:
|
|
warnings.append("未找到可用相邻杆塔坐标，无法推导线路方向纵坡/横坡")
|
|
|
|
windward_factor: float | None = None
|
|
if payload.wind_direction_deg is not None and aspect_deg is not None:
|
|
diff = _angle_difference_deg(aspect_deg, float(payload.wind_direction_deg))
|
|
windward_factor = max(0.0, math.cos(math.radians(diff)))
|
|
|
|
terrain_exposure_index = min(
|
|
1.0,
|
|
max(0.0, slope_deg / 45.0) * ((0.6 + 0.4 * windward_factor) if windward_factor is not None else 1.0),
|
|
)
|
|
|
|
quality_score, quality_level = _evaluate_terrain_quality(
|
|
dem_resolution_m=dem_resolution_m,
|
|
search_radius_m=float(payload.search_radius_m),
|
|
warnings=warnings,
|
|
)
|
|
|
|
metrics = LightningTowerTerrainMetrics(
|
|
slope_deg=round(slope_deg, 6),
|
|
aspect_deg=(round(aspect_deg, 6) if aspect_deg is not None else None),
|
|
slope_mean_deg=(round(slope_mean_deg, 6) if slope_mean_deg is not None else None),
|
|
slope_p95_deg=(round(slope_p95_deg, 6) if slope_p95_deg is not None else None),
|
|
slope_max_deg=(round(slope_max_deg, 6) if slope_max_deg is not None else None),
|
|
slope_along_line_deg=(round(slope_along_line_deg, 6) if slope_along_line_deg is not None else None),
|
|
slope_cross_line_deg=(round(slope_cross_line_deg, 6) if slope_cross_line_deg is not None else None),
|
|
relief_m_50=round(relief_m_50, 6),
|
|
dem_source=(payload.dem_source or "manual-grid"),
|
|
dem_resolution_m=round(dem_resolution_m, 6),
|
|
quality_score=round(quality_score, 2),
|
|
quality_level=quality_level,
|
|
terrain_exposure_index=round(terrain_exposure_index, 6),
|
|
windward_factor=(round(windward_factor, 6) if windward_factor is not None else None),
|
|
algorithm_version=TERRAIN_ALGORITHM_VERSION,
|
|
computed_at=utcnow(),
|
|
land_cover_type=payload.land_cover_type,
|
|
)
|
|
|
|
persisted = False
|
|
if payload.persist:
|
|
if tower is None:
|
|
raise HTTPException(status_code=status.HTTP_400_BAD_REQUEST, detail="persist=true 时必须传入 tower_id")
|
|
if not can_persist:
|
|
raise HTTPException(status_code=status.HTTP_403_FORBIDDEN, detail="缺少 tower.manage/lightning.manage 权限，无法持久化")
|
|
|
|
raw_extra = dict(tower.raw_extra_json or {})
|
|
terrain_payload = metrics.model_dump(mode="json", exclude_none=True)
|
|
terrain_payload["search_radius_m"] = float(payload.search_radius_m)
|
|
if line_azimuth_deg is not None:
|
|
terrain_payload["line_azimuth_deg"] = round(line_azimuth_deg, 6)
|
|
raw_extra["terrain_metrics"] = terrain_payload
|
|
|
|
tower.raw_extra_json = raw_extra
|
|
if slope_along_line_deg is not None:
|
|
tower.slope_1 = abs(float(slope_along_line_deg))
|
|
if slope_cross_line_deg is not None:
|
|
tower.slope_2 = abs(float(slope_cross_line_deg))
|
|
if payload.altitude_m is not None:
|
|
tower.altitude_m = payload.altitude_m
|
|
tower.update_user = actor_user_id
|
|
tower.update_date = utcnow()
|
|
db.commit()
|
|
persisted = True
|
|
|
|
return LightningTowerTerrainComputeResponse(
|
|
tower_id=tower.id if tower else None,
|
|
tower_no=tower.tower_no if tower else None,
|
|
line_id=tower.line_id if tower else None,
|
|
center_longitude=float(center_lon),
|
|
center_latitude=float(center_lat),
|
|
method="horn_3x3",
|
|
persisted=persisted,
|
|
terrain_metrics=metrics,
|
|
warnings=warnings,
|
|
)
|
|
|
|
|
|
def compare_measured_and_synthetic_distribution(
|
|
db: Session,
|
|
*,
|
|
min_lat: float | None,
|
|
max_lat: float | None,
|
|
min_lon: float | None,
|
|
max_lon: float | None,
|
|
region_id: str | None,
|
|
city: str | None,
|
|
location_tag: str | None,
|
|
grid_size_km: float,
|
|
years: float | None,
|
|
) -> LightningSyntheticCompareResponse:
|
|
base_filters = _build_distribution_filters(
|
|
min_lat=min_lat,
|
|
max_lat=max_lat,
|
|
min_lon=min_lon,
|
|
max_lon=max_lon,
|
|
region_id=region_id,
|
|
city=city,
|
|
location_tag=location_tag,
|
|
polarity=None,
|
|
is_synthetic=None,
|
|
)
|
|
|
|
extent_row = db.execute(
|
|
select(
|
|
func.min(LightningCurrentEvent.latitude).label("min_lat"),
|
|
func.max(LightningCurrentEvent.latitude).label("max_lat"),
|
|
func.min(LightningCurrentEvent.longitude).label("min_lon"),
|
|
func.max(LightningCurrentEvent.longitude).label("max_lon"),
|
|
).where(*base_filters)
|
|
).one()
|
|
|
|
if extent_row.min_lat is None or extent_row.max_lat is None or extent_row.min_lon is None or extent_row.max_lon is None:
|
|
empty_stats = LightningSyntheticDatasetStats(count=0)
|
|
return LightningSyntheticCompareResponse(
|
|
grid_size_km=grid_size_km,
|
|
data_years=years or 1.0,
|
|
measured=empty_stats,
|
|
synthetic=empty_stats,
|
|
grid_cosine_similarity=None,
|
|
note="当前筛选条件下没有可对比的数据",
|
|
)
|
|
|
|
min_lat_value = float(extent_row.min_lat)
|
|
max_lat_value = float(extent_row.max_lat)
|
|
min_lon_value = float(extent_row.min_lon)
|
|
max_lon_value = float(extent_row.max_lon)
|
|
center_lat = (min_lat_value + max_lat_value) / 2.0
|
|
km_per_lon = _safe_km_per_lon(center_lat)
|
|
width_km = max((max_lon_value - min_lon_value) * km_per_lon, grid_size_km)
|
|
height_km = max((max_lat_value - min_lat_value) * DEGREE_TO_KM, grid_size_km)
|
|
area_km2 = max(width_km * height_km, grid_size_km * grid_size_km)
|
|
data_years = years if years is not None else _infer_data_years(db, base_filters)
|
|
data_years = max(data_years, 1e-6)
|
|
|
|
measured_filters = [*base_filters, LightningCurrentEvent.is_synthetic.is_(False)]
|
|
synthetic_filters = [*base_filters, LightningCurrentEvent.is_synthetic.is_(True)]
|
|
measured = _aggregate_dataset_stats(db, measured_filters, area_km2=area_km2, data_years=data_years)
|
|
synthetic = _aggregate_dataset_stats(db, synthetic_filters, area_km2=area_km2, data_years=data_years)
|
|
|
|
similarity = _compute_grid_cosine_similarity(
|
|
db,
|
|
measured_filters=measured_filters,
|
|
synthetic_filters=synthetic_filters,
|
|
min_lat=min_lat_value,
|
|
min_lon=min_lon_value,
|
|
km_per_lon=km_per_lon,
|
|
grid_size_km=grid_size_km,
|
|
)
|
|
note = None
|
|
if measured.count == 0 or synthetic.count == 0:
|
|
note = "实测或合成数据为空，无法计算网格相似度"
|
|
|
|
return LightningSyntheticCompareResponse(
|
|
grid_size_km=grid_size_km,
|
|
data_years=data_years,
|
|
measured=measured,
|
|
synthetic=synthetic,
|
|
grid_cosine_similarity=similarity,
|
|
note=note,
|
|
)
|
|
|
|
|
|
def build_lightning_distribution_report(
|
|
db: Session,
|
|
*,
|
|
period: str,
|
|
anchor_time: datetime | None,
|
|
min_lat: float | None,
|
|
max_lat: float | None,
|
|
min_lon: float | None,
|
|
max_lon: float | None,
|
|
region_id: str | None,
|
|
city: str | None,
|
|
location_tag: str | None,
|
|
is_synthetic: bool | None,
|
|
) -> LightningDistributionReportResponse:
|
|
now = anchor_time or utcnow()
|
|
if period == "month":
|
|
days = 30
|
|
else:
|
|
period = "week"
|
|
days = 7
|
|
start_time = now - timedelta(days=days)
|
|
|
|
filters = _build_distribution_filters(
|
|
min_lat=min_lat,
|
|
max_lat=max_lat,
|
|
min_lon=min_lon,
|
|
max_lon=max_lon,
|
|
region_id=region_id,
|
|
city=city,
|
|
location_tag=location_tag,
|
|
polarity=None,
|
|
is_synthetic=is_synthetic,
|
|
)
|
|
event_ts_expr = func.coalesce(LightningCurrentEvent.event_time, LightningCurrentEvent.create_date)
|
|
filters.append(and_(event_ts_expr >= start_time, event_ts_expr <= now))
|
|
|
|
summary_row = db.execute(
|
|
select(
|
|
func.count().label("count"),
|
|
func.max(LightningCurrentEvent.peak_abs_current_ka).label("max_abs"),
|
|
func.avg(LightningCurrentEvent.peak_abs_current_ka).label("avg_abs"),
|
|
func.sum(case((LightningCurrentEvent.polarity == "positive", 1), else_=0)).label("positive_count"),
|
|
func.min(LightningCurrentEvent.latitude).label("min_lat"),
|
|
func.max(LightningCurrentEvent.latitude).label("max_lat"),
|
|
func.min(LightningCurrentEvent.longitude).label("min_lon"),
|
|
func.max(LightningCurrentEvent.longitude).label("max_lon"),
|
|
).where(*filters)
|
|
).one()
|
|
|
|
strike_count = int(summary_row.count or 0)
|
|
area_km2 = 1.0
|
|
if (
|
|
summary_row.min_lat is not None
|
|
and summary_row.max_lat is not None
|
|
and summary_row.min_lon is not None
|
|
and summary_row.max_lon is not None
|
|
):
|
|
center_lat = (float(summary_row.min_lat) + float(summary_row.max_lat)) / 2.0
|
|
km_per_lon = _safe_km_per_lon(center_lat)
|
|
width_km = max((float(summary_row.max_lon) - float(summary_row.min_lon)) * km_per_lon, 1.0)
|
|
height_km = max((float(summary_row.max_lat) - float(summary_row.min_lat)) * DEGREE_TO_KM, 1.0)
|
|
area_km2 = width_km * height_km
|
|
|
|
data_years = max(days / 365.25, 1e-6)
|
|
ng = strike_count / (area_km2 * data_years) if strike_count > 0 else 0.0
|
|
|
|
most_severe_row = db.execute(
|
|
select(
|
|
LightningCurrentEvent.id,
|
|
LightningCurrentEvent.event_id,
|
|
LightningCurrentEvent.longitude,
|
|
LightningCurrentEvent.latitude,
|
|
LightningCurrentEvent.peak_current_ka,
|
|
LightningCurrentEvent.peak_abs_current_ka,
|
|
LightningCurrentEvent.polarity,
|
|
LightningCurrentEvent.event_time,
|
|
LightningCurrentEvent.location_tag,
|
|
LightningCurrentEvent.city,
|
|
)
|
|
.where(*filters)
|
|
.order_by(
|
|
LightningCurrentEvent.peak_abs_current_ka.desc(),
|
|
event_ts_expr.desc(),
|
|
LightningCurrentEvent.id.desc(),
|
|
)
|
|
.limit(1)
|
|
).one_or_none()
|
|
|
|
severe_event: LightningDistributionEventBrief | None = None
|
|
if most_severe_row is not None:
|
|
severe_event = LightningDistributionEventBrief(
|
|
id=str(most_severe_row.id),
|
|
event_id=str(most_severe_row.event_id),
|
|
longitude=float(most_severe_row.longitude) if most_severe_row.longitude is not None else None,
|
|
latitude=float(most_severe_row.latitude) if most_severe_row.latitude is not None else None,
|
|
current_ka=float(most_severe_row.peak_current_ka) if most_severe_row.peak_current_ka is not None else None,
|
|
abs_current_ka=float(most_severe_row.peak_abs_current_ka) if most_severe_row.peak_abs_current_ka is not None else None,
|
|
polarity=str(most_severe_row.polarity), # type: ignore[arg-type]
|
|
event_time=most_severe_row.event_time,
|
|
location_tag=most_severe_row.location_tag,
|
|
city=most_severe_row.city,
|
|
)
|
|
|
|
positive_count = int(summary_row.positive_count or 0)
|
|
return LightningDistributionReportResponse(
|
|
period=period, # type: ignore[arg-type]
|
|
start_time=start_time,
|
|
end_time=now,
|
|
strike_count=strike_count,
|
|
max_abs_current_ka=float(summary_row.max_abs) if summary_row.max_abs is not None else None,
|
|
avg_abs_current_ka=float(summary_row.avg_abs) if summary_row.avg_abs is not None else None,
|
|
positive_ratio=(positive_count / strike_count) if strike_count > 0 else 0.0,
|
|
ng_per_km2_year=ng,
|
|
most_severe_event=severe_event,
|
|
)
|
|
|
|
|
|
def get_peak_exceedance_curve(
|
|
db: Session,
|
|
*,
|
|
region_id: str | None,
|
|
polarity: str | None,
|
|
wave_shape: str | None,
|
|
is_synthetic: bool | None,
|
|
thresholds_ka: list[float] | None,
|
|
default_points: int = 12,
|
|
) -> LightningCurrentExceedanceResponse:
|
|
filters: list[Any] = [LightningCurrentEvent.peak_abs_current_ka.is_not(None)]
|
|
|
|
normalized_region = _normalize_str(region_id)
|
|
if normalized_region:
|
|
filters.append(LightningCurrentEvent.region_id == normalized_region)
|
|
normalized_polarity = _normalize_str(polarity)
|
|
if normalized_polarity:
|
|
filters.append(LightningCurrentEvent.polarity == normalized_polarity.lower())
|
|
normalized_wave_shape = _normalize_str(wave_shape)
|
|
if normalized_wave_shape:
|
|
filters.append(func.lower(LightningCurrentEvent.wave_shape) == normalized_wave_shape.lower())
|
|
if is_synthetic is not None:
|
|
filters.append(LightningCurrentEvent.is_synthetic == is_synthetic)
|
|
|
|
values = db.execute(select(LightningCurrentEvent.peak_abs_current_ka).where(*filters)).scalars().all()
|
|
peaks = [float(item) for item in values if item is not None]
|
|
total = len(peaks)
|
|
if total == 0:
|
|
return LightningCurrentExceedanceResponse(total_events=0, thresholds=[])
|
|
|
|
manual_thresholds = sorted({float(value) for value in (thresholds_ka or []) if value > 0})
|
|
if manual_thresholds:
|
|
thresholds = manual_thresholds
|
|
else:
|
|
max_peak = max(peaks)
|
|
step = max(1.0, round(max_peak / max(default_points, 1), 3))
|
|
thresholds = [round(step * idx, 3) for idx in range(1, default_points + 1)]
|
|
|
|
points: list[LightningCurrentExceedancePoint] = []
|
|
for threshold in thresholds:
|
|
exceedance_count = sum(1 for value in peaks if value >= threshold)
|
|
probability = exceedance_count / total
|
|
points.append(
|
|
LightningCurrentExceedancePoint(
|
|
threshold_ka=threshold,
|
|
exceedance_probability=probability,
|
|
exceedance_count=exceedance_count,
|
|
)
|
|
)
|
|
|
|
return LightningCurrentExceedanceResponse(total_events=total, thresholds=points)
|
|
|
|
|
|
def _build_distribution_filters(
|
|
*,
|
|
min_lat: float | None,
|
|
max_lat: float | None,
|
|
min_lon: float | None,
|
|
max_lon: float | None,
|
|
region_id: str | None,
|
|
city: str | None,
|
|
location_tag: str | None,
|
|
polarity: str | None,
|
|
is_synthetic: bool | None,
|
|
) -> list[Any]:
|
|
filters: list[Any] = [
|
|
LightningCurrentEvent.latitude.is_not(None),
|
|
LightningCurrentEvent.longitude.is_not(None),
|
|
LightningCurrentEvent.peak_abs_current_ka.is_not(None),
|
|
]
|
|
if min_lat is not None:
|
|
filters.append(LightningCurrentEvent.latitude >= min_lat)
|
|
if max_lat is not None:
|
|
filters.append(LightningCurrentEvent.latitude <= max_lat)
|
|
if min_lon is not None:
|
|
filters.append(LightningCurrentEvent.longitude >= min_lon)
|
|
if max_lon is not None:
|
|
filters.append(LightningCurrentEvent.longitude <= max_lon)
|
|
|
|
normalized_region = _normalize_str(region_id)
|
|
if normalized_region:
|
|
filters.append(LightningCurrentEvent.region_id == normalized_region)
|
|
|
|
normalized_city = _normalize_str(city)
|
|
if normalized_city:
|
|
filters.append(func.lower(LightningCurrentEvent.city) == normalized_city.lower())
|
|
|
|
normalized_location = _normalize_str(location_tag)
|
|
if normalized_location:
|
|
filters.append(LightningCurrentEvent.location_tag.ilike(f"%{normalized_location}%"))
|
|
|
|
normalized_polarity = _normalize_str(polarity)
|
|
if normalized_polarity:
|
|
filters.append(LightningCurrentEvent.polarity == normalized_polarity.lower())
|
|
|
|
if is_synthetic is not None:
|
|
filters.append(LightningCurrentEvent.is_synthetic == is_synthetic)
|
|
|
|
return filters
|
|
|
|
|
|
def _build_p_curve_with_filters(
|
|
db: Session,
|
|
*,
|
|
filters: list[Any],
|
|
total_events: int,
|
|
max_peak: float,
|
|
thresholds_ka: list[float] | None,
|
|
) -> list[LightningCurrentExceedancePoint]:
|
|
if total_events <= 0:
|
|
return []
|
|
manual_thresholds = sorted({float(value) for value in (thresholds_ka or []) if value > 0})
|
|
if manual_thresholds:
|
|
thresholds = manual_thresholds
|
|
else:
|
|
default_points = 12
|
|
step = max(1.0, round(max_peak / max(default_points, 1), 3))
|
|
thresholds = [round(step * idx, 3) for idx in range(1, default_points + 1)]
|
|
|
|
points: list[LightningCurrentExceedancePoint] = []
|
|
for threshold in thresholds:
|
|
exceed_count = int(
|
|
db.scalar(
|
|
select(func.count()).select_from(LightningCurrentEvent).where(
|
|
*filters,
|
|
LightningCurrentEvent.peak_abs_current_ka >= threshold,
|
|
)
|
|
)
|
|
or 0
|
|
)
|
|
points.append(
|
|
LightningCurrentExceedancePoint(
|
|
threshold_ka=threshold,
|
|
exceedance_probability=(exceed_count / total_events) if total_events > 0 else 0.0,
|
|
exceedance_count=exceed_count,
|
|
)
|
|
)
|
|
return points
|
|
|
|
|
|
def _aggregate_dataset_stats(
|
|
db: Session,
|
|
filters: list[Any],
|
|
*,
|
|
area_km2: float,
|
|
data_years: float,
|
|
) -> LightningSyntheticDatasetStats:
|
|
row = db.execute(
|
|
select(
|
|
func.count().label("count"),
|
|
func.max(LightningCurrentEvent.peak_abs_current_ka).label("max_abs"),
|
|
func.avg(LightningCurrentEvent.peak_abs_current_ka).label("avg_abs"),
|
|
func.sum(case((LightningCurrentEvent.polarity == "positive", 1), else_=0)).label("positive_count"),
|
|
).where(*filters)
|
|
).one()
|
|
count = int(row.count or 0)
|
|
positive_count = int(row.positive_count or 0)
|
|
ng = count / (area_km2 * data_years) if count > 0 else 0.0
|
|
return LightningSyntheticDatasetStats(
|
|
count=count,
|
|
max_abs_current_ka=float(row.max_abs) if row.max_abs is not None else None,
|
|
avg_abs_current_ka=float(row.avg_abs) if row.avg_abs is not None else None,
|
|
positive_ratio=(positive_count / count) if count > 0 else 0.0,
|
|
ng_per_km2_year=ng,
|
|
)
|
|
|
|
|
|
def _compute_grid_cosine_similarity(
|
|
db: Session,
|
|
*,
|
|
measured_filters: list[Any],
|
|
synthetic_filters: list[Any],
|
|
min_lat: float,
|
|
min_lon: float,
|
|
km_per_lon: float,
|
|
grid_size_km: float,
|
|
) -> float | None:
|
|
cell_x_expr = func.floor(((LightningCurrentEvent.longitude - min_lon) * km_per_lon) / grid_size_km)
|
|
cell_y_expr = func.floor(((LightningCurrentEvent.latitude - min_lat) * DEGREE_TO_KM) / grid_size_km)
|
|
measured_rows = db.execute(
|
|
select(cell_x_expr.label("x"), cell_y_expr.label("y"), func.count().label("c"))
|
|
.where(*measured_filters)
|
|
.group_by(cell_x_expr, cell_y_expr)
|
|
).all()
|
|
synthetic_rows = db.execute(
|
|
select(cell_x_expr.label("x"), cell_y_expr.label("y"), func.count().label("c"))
|
|
.where(*synthetic_filters)
|
|
.group_by(cell_x_expr, cell_y_expr)
|
|
).all()
|
|
|
|
if not measured_rows or not synthetic_rows:
|
|
return None
|
|
|
|
measured_map = {(int(float(item.x)), int(float(item.y))): float(item.c) for item in measured_rows}
|
|
synthetic_map = {(int(float(item.x)), int(float(item.y))): float(item.c) for item in synthetic_rows}
|
|
|
|
all_keys = set(measured_map.keys()) | set(synthetic_map.keys())
|
|
dot = sum(measured_map.get(key, 0.0) * synthetic_map.get(key, 0.0) for key in all_keys)
|
|
norm_m = math.sqrt(sum(value * value for value in measured_map.values()))
|
|
norm_s = math.sqrt(sum(value * value for value in synthetic_map.values()))
|
|
if norm_m <= 0 or norm_s <= 0:
|
|
return None
|
|
return max(0.0, min(1.0, dot / (norm_m * norm_s)))
|
|
|
|
|
|
def _parse_distribution_line(raw_line: str) -> tuple[float, float, float] | None:
|
|
stripped = raw_line.strip()
|
|
if not stripped:
|
|
return None
|
|
tokens = [token for token in TOKEN_SPLIT_PATTERN.split(stripped) if token]
|
|
if len(tokens) < 3:
|
|
return None
|
|
values: list[float] = []
|
|
for token in tokens:
|
|
parsed = _parse_float(token)
|
|
if parsed is not None:
|
|
values.append(parsed)
|
|
if len(values) >= 3:
|
|
break
|
|
if len(values) < 3:
|
|
return None
|
|
return (values[0], values[1], values[2])
|
|
|
|
|
|
def _safe_km_per_lon(latitude: float) -> float:
|
|
factor = abs(math.cos(math.radians(latitude)))
|
|
return max(DEGREE_TO_KM * factor, 1e-6)
|
|
|
|
|
|
def _haversine_km(lat1: float, lon1: float, lat2: float, lon2: float) -> float:
|
|
radius_km = 6371.0
|
|
lat1_rad = math.radians(lat1)
|
|
lat2_rad = math.radians(lat2)
|
|
delta_lat = lat2_rad - lat1_rad
|
|
delta_lon = math.radians(lon2 - lon1)
|
|
a = (
|
|
math.sin(delta_lat / 2.0) ** 2
|
|
+ math.cos(lat1_rad) * math.cos(lat2_rad) * (math.sin(delta_lon / 2.0) ** 2)
|
|
)
|
|
c = 2.0 * math.atan2(math.sqrt(a), math.sqrt(max(1.0 - a, 0.0)))
|
|
return radius_km * c
|
|
|
|
|
|
def _resolve_data_years_from_timestamps(values: list[datetime]) -> float:
|
|
if len(values) < 2:
|
|
return 1.0
|
|
minimum = min(values)
|
|
maximum = max(values)
|
|
delta_days = (maximum - minimum).total_seconds() / 86400.0
|
|
return max(delta_days / 365.25, 1.0)
|
|
|
|
|
|
def _infer_data_years(db: Session, filters: list[Any]) -> float:
|
|
row = db.execute(
|
|
select(
|
|
func.min(LightningCurrentEvent.event_time).label("min_time"),
|
|
func.max(LightningCurrentEvent.event_time).label("max_time"),
|
|
).where(*filters, LightningCurrentEvent.event_time.is_not(None))
|
|
).one()
|
|
if row.min_time is None or row.max_time is None:
|
|
return 1.0
|
|
delta_days = (row.max_time - row.min_time).total_seconds() / 86400.0
|
|
return max(delta_days / 365.25, 1.0)
|
|
|
|
|
|
def _build_tower_terrain_metrics_from_tower(tower: LineTower | None) -> LightningTowerTerrainMetrics | None:
|
|
if tower is None:
|
|
return None
|
|
|
|
raw_extra = tower.raw_extra_json or {}
|
|
terrain_payload = raw_extra.get("terrain_metrics") if isinstance(raw_extra, dict) else None
|
|
if isinstance(terrain_payload, dict):
|
|
try:
|
|
return LightningTowerTerrainMetrics.model_validate(terrain_payload)
|
|
except Exception:
|
|
pass
|
|
|
|
slope_along = float(tower.slope_1) if tower.slope_1 is not None else None
|
|
slope_cross = float(tower.slope_2) if tower.slope_2 is not None else None
|
|
if slope_along is None and slope_cross is None:
|
|
return None
|
|
|
|
slope_candidates = [abs(value) for value in (slope_along, slope_cross) if value is not None]
|
|
slope_deg = max(slope_candidates) if slope_candidates else None
|
|
terrain_exposure = min(1.0, (slope_deg / 45.0)) if slope_deg is not None else None
|
|
return LightningTowerTerrainMetrics(
|
|
slope_deg=slope_deg,
|
|
slope_max_deg=slope_deg,
|
|
slope_along_line_deg=slope_along,
|
|
slope_cross_line_deg=slope_cross,
|
|
terrain_exposure_index=terrain_exposure,
|
|
algorithm_version="tower-legacy",
|
|
)
|
|
|
|
|
|
def _compute_horn_gradient(grid_m: list[list[float]], *, cell_size_m: float) -> tuple[float, float]:
|
|
if len(grid_m) != 3 or any(len(row) != 3 for row in grid_m):
|
|
raise HTTPException(status_code=status.HTTP_400_BAD_REQUEST, detail="DEM 网格必须是 3x3")
|
|
if cell_size_m <= 0:
|
|
raise HTTPException(status_code=status.HTTP_400_BAD_REQUEST, detail="cell_size_m 必须大于 0")
|
|
|
|
z1, z2, z3 = float(grid_m[0][0]), float(grid_m[0][1]), float(grid_m[0][2])
|
|
z4, z5, z6 = float(grid_m[1][0]), float(grid_m[1][1]), float(grid_m[1][2])
|
|
z7, z8, z9 = float(grid_m[2][0]), float(grid_m[2][1]), float(grid_m[2][2])
|
|
|
|
dz_dx = ((z3 + 2.0 * z6 + z9) - (z1 + 2.0 * z4 + z7)) / (8.0 * cell_size_m)
|
|
dz_dy = ((z1 + 2.0 * z2 + z3) - (z7 + 2.0 * z8 + z9)) / (8.0 * cell_size_m)
|
|
return dz_dx, dz_dy
|
|
|
|
|
|
def _compute_aspect_deg(dz_dx: float, dz_dy: float) -> float | None:
|
|
if math.isclose(dz_dx, 0.0, abs_tol=1e-12) and math.isclose(dz_dy, 0.0, abs_tol=1e-12):
|
|
return None
|
|
downslope_east = -dz_dx
|
|
downslope_north = -dz_dy
|
|
return (math.degrees(math.atan2(downslope_east, downslope_north)) + 360.0) % 360.0
|
|
|
|
|
|
def _compute_neighbor_slopes(grid_m: list[list[float]], *, cell_size_m: float) -> list[float]:
|
|
center = float(grid_m[1][1])
|
|
slopes: list[float] = []
|
|
for row in range(3):
|
|
for col in range(3):
|
|
if row == 1 and col == 1:
|
|
continue
|
|
delta_h = float(grid_m[row][col]) - center
|
|
distance = cell_size_m * math.sqrt(float((row - 1) ** 2 + (col - 1) ** 2))
|
|
if distance <= 0:
|
|
continue
|
|
slopes.append(math.degrees(math.atan(abs(delta_h) / distance)))
|
|
return slopes
|
|
|
|
|
|
def _percentile(values: list[float], quantile: float) -> float:
|
|
if not values:
|
|
return 0.0
|
|
if quantile <= 0:
|
|
return min(values)
|
|
if quantile >= 1:
|
|
return max(values)
|
|
|
|
sorted_values = sorted(values)
|
|
position = (len(sorted_values) - 1) * quantile
|
|
lower = int(math.floor(position))
|
|
upper = int(math.ceil(position))
|
|
if lower == upper:
|
|
return float(sorted_values[lower])
|
|
lower_value = float(sorted_values[lower])
|
|
upper_value = float(sorted_values[upper])
|
|
weight = position - lower
|
|
return lower_value * (1.0 - weight) + upper_value * weight
|
|
|
|
|
|
def _infer_tower_line_azimuth_deg(db: Session, tower: LineTower | None) -> float | None:
|
|
if (
|
|
tower is None
|
|
or tower.line_id is None
|
|
or tower.seq_no is None
|
|
or tower.longitude is None
|
|
or tower.latitude is None
|
|
):
|
|
return None
|
|
|
|
rows = db.execute(
|
|
select(LineTower.seq_no, LineTower.longitude, LineTower.latitude)
|
|
.where(
|
|
LineTower.line_id == tower.line_id,
|
|
LineTower.longitude.is_not(None),
|
|
LineTower.latitude.is_not(None),
|
|
)
|
|
.order_by(LineTower.seq_no.asc())
|
|
).all()
|
|
if not rows:
|
|
return None
|
|
|
|
previous: tuple[float, float] | None = None
|
|
next_point: tuple[float, float] | None = None
|
|
for row in rows:
|
|
seq_no = int(row.seq_no)
|
|
point = (float(row.longitude), float(row.latitude))
|
|
if seq_no < tower.seq_no:
|
|
previous = point
|
|
continue
|
|
if seq_no > tower.seq_no:
|
|
next_point = point
|
|
break
|
|
|
|
if previous is not None and next_point is not None:
|
|
start_lon, start_lat = previous
|
|
end_lon, end_lat = next_point
|
|
elif next_point is not None:
|
|
start_lon = float(tower.longitude)
|
|
start_lat = float(tower.latitude)
|
|
end_lon, end_lat = next_point
|
|
elif previous is not None:
|
|
start_lon, start_lat = previous
|
|
end_lon = float(tower.longitude)
|
|
end_lat = float(tower.latitude)
|
|
else:
|
|
return None
|
|
|
|
if math.isclose(start_lon, end_lon, abs_tol=1e-12) and math.isclose(start_lat, end_lat, abs_tol=1e-12):
|
|
return None
|
|
|
|
start_lat_rad = math.radians(start_lat)
|
|
end_lat_rad = math.radians(end_lat)
|
|
delta_lon_rad = math.radians(end_lon - start_lon)
|
|
x = math.sin(delta_lon_rad) * math.cos(end_lat_rad)
|
|
y = (
|
|
math.cos(start_lat_rad) * math.sin(end_lat_rad)
|
|
- math.sin(start_lat_rad) * math.cos(end_lat_rad) * math.cos(delta_lon_rad)
|
|
)
|
|
if math.isclose(x, 0.0, abs_tol=1e-12) and math.isclose(y, 0.0, abs_tol=1e-12):
|
|
return None
|
|
return (math.degrees(math.atan2(x, y)) + 360.0) % 360.0
|
|
|
|
|
|
def _directional_gradient(dz_dx: float, dz_dy: float, azimuth_deg: float) -> float:
|
|
azimuth_rad = math.radians(azimuth_deg % 360.0)
|
|
direction_east = math.sin(azimuth_rad)
|
|
direction_north = math.cos(azimuth_rad)
|
|
return dz_dx * direction_east + dz_dy * direction_north
|
|
|
|
|
|
def _angle_difference_deg(angle_1_deg: float, angle_2_deg: float) -> float:
|
|
return abs((angle_1_deg - angle_2_deg + 180.0) % 360.0 - 180.0)
|
|
|
|
|
|
def _evaluate_terrain_quality(
|
|
*,
|
|
dem_resolution_m: float,
|
|
search_radius_m: float,
|
|
warnings: list[str],
|
|
) -> tuple[float, str]:
|
|
score = 100.0
|
|
|
|
def add_warning(message: str) -> None:
|
|
if message not in warnings:
|
|
warnings.append(message)
|
|
|
|
if dem_resolution_m <= 5.0:
|
|
score -= 0.0
|
|
elif dem_resolution_m <= 10.0:
|
|
score -= 5.0
|
|
elif dem_resolution_m <= 12.5:
|
|
score -= 10.0
|
|
elif dem_resolution_m <= 30.0:
|
|
score -= 25.0
|
|
add_warning("DEM 分辨率大于 12.5m，微地形倾角结果可能偏平滑。")
|
|
elif dem_resolution_m <= 90.0:
|
|
score -= 45.0
|
|
add_warning("DEM 分辨率较低（30-90m），建议使用 10m 及以上数据复核。")
|
|
else:
|
|
score -= 60.0
|
|
add_warning("DEM 分辨率超过 90m，当前结果仅可用于粗略评估。")
|
|
|
|
if search_radius_m < 30.0:
|
|
score -= 15.0
|
|
add_warning("地形检索半径过小（<30m），建议提高到 50m 左右。")
|
|
elif search_radius_m < 50.0:
|
|
score -= 8.0
|
|
elif search_radius_m > 300.0:
|
|
score -= 10.0
|
|
add_warning("地形检索半径较大（>300m），局部坡度可能被稀释。")
|
|
elif search_radius_m > 150.0:
|
|
score -= 5.0
|
|
|
|
score = max(0.0, min(100.0, score))
|
|
if score >= 80.0:
|
|
quality_level = "HIGH"
|
|
elif score >= 60.0:
|
|
quality_level = "MEDIUM"
|
|
else:
|
|
quality_level = "LOW"
|
|
return score, quality_level
|
|
|
|
|
|
def _derive_tower_risk_level(
|
|
*,
|
|
strike_count: int,
|
|
exceed_design_count: int,
|
|
max_abs_current_ka: float | None,
|
|
design_current_ka: float,
|
|
ng_per_km2_year: float,
|
|
) -> str:
|
|
if strike_count <= 0 or max_abs_current_ka is None:
|
|
return "LOW"
|
|
if exceed_design_count >= 3 or max_abs_current_ka >= design_current_ka * 1.5 or ng_per_km2_year >= 8.0:
|
|
return "HIGH"
|
|
if exceed_design_count >= 1 or max_abs_current_ka >= design_current_ka or ng_per_km2_year >= 4.0:
|
|
return "MEDIUM"
|
|
return "LOW"
|
|
|
|
|
|
def _tower_risk_recommendation(risk_level: str) -> str:
|
|
if risk_level == "HIGH":
|
|
return "建议立即复核绝缘配置与避雷器参数，并安排现场巡检。"
|
|
if risk_level == "MEDIUM":
|
|
return "建议纳入近期巡检计划，核查接地与线路防雷参数。"
|
|
return "当前风险可控，建议按常规周期继续监测。"
|
|
|
|
|
|
def _generate_distribution_event_id() -> str:
|
|
return f"LD-{datetime.now().strftime('%Y%m%d%H%M%S')}-{uuid4().hex[:8]}"
|
|
|
|
|
|
def _parse_current_series(text: str) -> ParsedSeries:
|
|
rows: list[tuple[float | None, float]] = []
|
|
warnings: list[str] = []
|
|
explicit_time_rows = 0
|
|
|
|
for line_no, raw_line in enumerate(text.splitlines(), start=1):
|
|
stripped = raw_line.strip()
|
|
if not stripped:
|
|
continue
|
|
|
|
tokens = [token for token in TOKEN_SPLIT_PATTERN.split(stripped) if token]
|
|
if not tokens:
|
|
continue
|
|
|
|
parsed_pair = _parse_two_numbers(tokens)
|
|
if parsed_pair is not None:
|
|
rows.append(parsed_pair)
|
|
explicit_time_rows += 1
|
|
continue
|
|
|
|
current_value = _parse_last_number(tokens)
|
|
if current_value is None:
|
|
if len(warnings) < 20:
|
|
warnings.append(f"第 {line_no} 行未识别为数值，已跳过")
|
|
continue
|
|
rows.append((None, current_value))
|
|
|
|
if not rows:
|
|
raise HTTPException(status_code=status.HTTP_400_BAD_REQUEST, detail="文件中未解析到有效数值序列")
|
|
|
|
currents = [item[1] for item in rows]
|
|
time_axis: list[float] | None = None
|
|
if explicit_time_rows == len(rows):
|
|
time_axis = [float(item[0]) for item in rows if item[0] is not None]
|
|
elif explicit_time_rows > 0:
|
|
warnings.append("检测到部分行存在时间列但并不完整，已统一按采样间隔重建时间轴")
|
|
|
|
return ParsedSeries(currents_ka=currents, time_us=time_axis, warnings=warnings)
|
|
|
|
|
|
def _parse_two_numbers(tokens: list[str]) -> tuple[float, float] | None:
|
|
if len(tokens) < 2:
|
|
return None
|
|
first = _parse_float(tokens[0])
|
|
second = _parse_float(tokens[1])
|
|
if first is None or second is None:
|
|
return None
|
|
return (first, second)
|
|
|
|
|
|
def _parse_last_number(tokens: list[str]) -> float | None:
|
|
for token in reversed(tokens):
|
|
parsed = _parse_float(token)
|
|
if parsed is not None:
|
|
return parsed
|
|
return None
|
|
|
|
|
|
def _build_time_axis(
|
|
raw_time_us: list[float] | None,
|
|
*,
|
|
sample_interval_us: float | None,
|
|
sample_count: int,
|
|
) -> tuple[list[float], float, float | None, list[str]]:
|
|
warnings: list[str] = []
|
|
if raw_time_us and len(raw_time_us) == sample_count:
|
|
base_time = raw_time_us[0]
|
|
normalized = [float(item - base_time) for item in raw_time_us]
|
|
positive_deltas: list[float] = []
|
|
monotonic = True
|
|
for idx in range(1, len(normalized)):
|
|
delta = normalized[idx] - normalized[idx - 1]
|
|
if delta <= 0:
|
|
monotonic = False
|
|
break
|
|
positive_deltas.append(delta)
|
|
|
|
if monotonic and positive_deltas:
|
|
inferred_interval = _median(positive_deltas)
|
|
inferred_frequency = 1_000_000.0 / inferred_interval if inferred_interval > 0 else None
|
|
return normalized, inferred_interval, inferred_frequency, warnings
|
|
|
|
warnings.append("时间列非严格递增，已改用固定采样间隔重建时间轴")
|
|
|
|
interval = sample_interval_us if sample_interval_us is not None else 1.0
|
|
if interval <= 0:
|
|
warnings.append("采样间隔必须大于 0，已回退到 1.0us")
|
|
interval = 1.0
|
|
if sample_interval_us is None:
|
|
warnings.append("未提供采样间隔，默认按 1.0us 计算")
|
|
|
|
axis = [float(idx) * interval for idx in range(sample_count)]
|
|
fs = 1_000_000.0 / interval if interval > 0 else None
|
|
return axis, interval, fs, warnings
|
|
|
|
|
|
def _extract_wave_features(currents_ka: list[float], time_axis_us: list[float], sample_interval_us: float) -> WaveFeatures:
|
|
sample_count = len(currents_ka)
|
|
if sample_count == 0:
|
|
return WaveFeatures(
|
|
peak_current_ka=None,
|
|
peak_abs_current_ka=None,
|
|
wavefront_time_t1_us=None,
|
|
half_value_time_t2_us=None,
|
|
steepness_ka_per_us=None,
|
|
action_integral_j_ohm=None,
|
|
wave_shape=None,
|
|
polarity="unknown",
|
|
stroke_count=0,
|
|
stroke_peaks_json=[],
|
|
)
|
|
|
|
max_val = max(currents_ka)
|
|
min_val = min(currents_ka)
|
|
|
|
polarity: LightningPolarity
|
|
if math.isclose(max_val, 0.0, abs_tol=1e-12) and math.isclose(min_val, 0.0, abs_tol=1e-12):
|
|
polarity = "unknown"
|
|
elif abs(max_val) > abs(min_val):
|
|
polarity = "positive"
|
|
elif abs(min_val) > abs(max_val):
|
|
polarity = "negative"
|
|
else:
|
|
polarity = "mixed"
|
|
|
|
peak_idx = max(range(sample_count), key=lambda idx: abs(currents_ka[idx]))
|
|
peak_current = float(currents_ka[peak_idx])
|
|
peak_abs = abs(peak_current)
|
|
|
|
dominant_sign = 1.0 if peak_current >= 0 else -1.0
|
|
aligned = [dominant_sign * value for value in currents_ka]
|
|
aligned_peak = aligned[peak_idx]
|
|
|
|
steepness = _compute_max_steepness(aligned, time_axis_us, peak_idx)
|
|
t10 = _find_up_crossing_time(time_axis_us, aligned, 0.1 * aligned_peak, peak_idx)
|
|
t90 = _find_up_crossing_time(time_axis_us, aligned, 0.9 * aligned_peak, peak_idx)
|
|
wavefront_time_t1: float | None = None
|
|
t0_virtual: float | None = None
|
|
if t10 is not None and t90 is not None and t90 > t10:
|
|
wavefront_time_t1 = 1.25 * (t90 - t10)
|
|
t0_virtual = t10 - 0.125 * (t90 - t10)
|
|
|
|
t50 = _find_down_crossing_time(time_axis_us, aligned, 0.5 * aligned_peak, peak_idx)
|
|
half_value_time_t2: float | None = None
|
|
if t50 is not None:
|
|
if t0_virtual is not None:
|
|
half_value_time_t2 = max(t50 - t0_virtual, 0.0)
|
|
else:
|
|
half_value_time_t2 = max(t50 - time_axis_us[0], 0.0)
|
|
|
|
action_integral = _compute_action_integral(currents_ka, time_axis_us, fallback_interval=sample_interval_us)
|
|
wave_shape = _classify_wave_shape(wavefront_time_t1, half_value_time_t2)
|
|
stroke_peaks = _detect_multiple_strokes(
|
|
currents_ka=currents_ka,
|
|
time_axis_us=time_axis_us,
|
|
dominant_sign=dominant_sign,
|
|
peak_abs=peak_abs,
|
|
sample_interval_us=sample_interval_us,
|
|
)
|
|
|
|
return WaveFeatures(
|
|
peak_current_ka=peak_current,
|
|
peak_abs_current_ka=peak_abs,
|
|
wavefront_time_t1_us=wavefront_time_t1,
|
|
half_value_time_t2_us=half_value_time_t2,
|
|
steepness_ka_per_us=steepness,
|
|
action_integral_j_ohm=action_integral,
|
|
wave_shape=wave_shape,
|
|
polarity=polarity,
|
|
stroke_count=len(stroke_peaks),
|
|
stroke_peaks_json=stroke_peaks,
|
|
)
|
|
|
|
|
|
def _compute_max_steepness(aligned: list[float], time_axis_us: list[float], peak_idx: int) -> float | None:
|
|
max_slope: float | None = None
|
|
for idx in range(1, peak_idx + 1):
|
|
dt = time_axis_us[idx] - time_axis_us[idx - 1]
|
|
if dt <= 0:
|
|
continue
|
|
slope = (aligned[idx] - aligned[idx - 1]) / dt
|
|
if max_slope is None or slope > max_slope:
|
|
max_slope = slope
|
|
return max_slope
|
|
|
|
|
|
def _find_up_crossing_time(
|
|
time_axis_us: list[float],
|
|
signal: list[float],
|
|
threshold: float,
|
|
end_idx: int,
|
|
) -> float | None:
|
|
if end_idx <= 0:
|
|
return None
|
|
for idx in range(1, end_idx + 1):
|
|
left = signal[idx - 1]
|
|
right = signal[idx]
|
|
if left == threshold:
|
|
return time_axis_us[idx - 1]
|
|
if left < threshold <= right:
|
|
return _interpolate_time(
|
|
time_axis_us[idx - 1],
|
|
time_axis_us[idx],
|
|
left,
|
|
right,
|
|
threshold,
|
|
)
|
|
return None
|
|
|
|
|
|
def _find_down_crossing_time(
|
|
time_axis_us: list[float],
|
|
signal: list[float],
|
|
threshold: float,
|
|
start_idx: int,
|
|
) -> float | None:
|
|
if start_idx >= len(signal) - 1:
|
|
return None
|
|
for idx in range(start_idx + 1, len(signal)):
|
|
left = signal[idx - 1]
|
|
right = signal[idx]
|
|
if left == threshold:
|
|
return time_axis_us[idx - 1]
|
|
if left >= threshold > right:
|
|
return _interpolate_time(
|
|
time_axis_us[idx - 1],
|
|
time_axis_us[idx],
|
|
left,
|
|
right,
|
|
threshold,
|
|
)
|
|
return None
|
|
|
|
|
|
def _interpolate_time(x1: float, x2: float, y1: float, y2: float, target: float) -> float:
|
|
if math.isclose(y1, y2):
|
|
return x2
|
|
ratio = (target - y1) / (y2 - y1)
|
|
return x1 + ratio * (x2 - x1)
|
|
|
|
|
|
def _compute_action_integral(currents_ka: list[float], time_axis_us: list[float], *, fallback_interval: float) -> float:
|
|
if len(currents_ka) < 2:
|
|
return 0.0
|
|
|
|
total = 0.0
|
|
for idx in range(1, len(currents_ka)):
|
|
dt = time_axis_us[idx] - time_axis_us[idx - 1]
|
|
if dt <= 0:
|
|
dt = fallback_interval
|
|
average_i2 = ((currents_ka[idx - 1] ** 2) + (currents_ka[idx] ** 2)) / 2.0
|
|
total += average_i2 * dt
|
|
return total
|
|
|
|
|
|
def _classify_wave_shape(t1_us: float | None, t2_us: float | None) -> str | None:
|
|
if t1_us is None or t2_us is None or t1_us <= 0 or t2_us <= 0:
|
|
return None
|
|
|
|
best_name: str | None = None
|
|
best_score = float("inf")
|
|
for name, std_t1, std_t2 in STANDARD_WAVE_SHAPES:
|
|
score = ((t1_us - std_t1) / std_t1) ** 2 + ((t2_us - std_t2) / std_t2) ** 2
|
|
if score < best_score:
|
|
best_name = name
|
|
best_score = score
|
|
return best_name
|
|
|
|
|
|
def _detect_multiple_strokes(
|
|
*,
|
|
currents_ka: list[float],
|
|
time_axis_us: list[float],
|
|
dominant_sign: float,
|
|
peak_abs: float,
|
|
sample_interval_us: float,
|
|
) -> list[dict[str, Any]]:
|
|
if not currents_ka:
|
|
return []
|
|
|
|
aligned = [dominant_sign * item for item in currents_ka]
|
|
threshold = max(peak_abs * 0.3, 5.0)
|
|
min_gap_us = 80.0
|
|
min_gap_samples = max(1, int(round(min_gap_us / max(sample_interval_us, 1e-9))))
|
|
|
|
peak_indices: list[int] = []
|
|
for idx in range(1, len(aligned) - 1):
|
|
current = aligned[idx]
|
|
if current < threshold:
|
|
continue
|
|
if current >= aligned[idx - 1] and current > aligned[idx + 1]:
|
|
if peak_indices and idx - peak_indices[-1] < min_gap_samples:
|
|
if current > aligned[peak_indices[-1]]:
|
|
peak_indices[-1] = idx
|
|
continue
|
|
peak_indices.append(idx)
|
|
|
|
if not peak_indices:
|
|
dominant_idx = max(range(len(currents_ka)), key=lambda idx: abs(currents_ka[idx]))
|
|
peak_indices = [dominant_idx]
|
|
|
|
limited = peak_indices[:20]
|
|
return [
|
|
{
|
|
"seq_no": order + 1,
|
|
"sample_index": idx + 1,
|
|
"time_us": time_axis_us[idx],
|
|
"current_ka": currents_ka[idx],
|
|
}
|
|
for order, idx in enumerate(limited)
|
|
]
|
|
|
|
|
|
def _decode_text_bytes(content: bytes) -> str:
|
|
for encoding in TEXT_ENCODINGS:
|
|
try:
|
|
return content.decode(encoding)
|
|
except UnicodeDecodeError:
|
|
continue
|
|
raise HTTPException(status_code=status.HTTP_400_BAD_REQUEST, detail="文件编码不受支持")
|
|
|
|
|
|
def _parse_float(value: Any) -> float | None:
|
|
if value is None:
|
|
return None
|
|
normalized = str(value).strip()
|
|
if not normalized:
|
|
return None
|
|
try:
|
|
return float(normalized)
|
|
except ValueError:
|
|
return None
|
|
|
|
|
|
def _normalize_str(value: Any) -> str | None:
|
|
if value is None:
|
|
return None
|
|
normalized = str(value).strip()
|
|
if not normalized:
|
|
return None
|
|
return normalized
|
|
|
|
|
|
def _median(values: list[float]) -> float:
|
|
if not values:
|
|
return 0.0
|
|
sorted_values = sorted(values)
|
|
n = len(sorted_values)
|
|
mid = n // 2
|
|
if n % 2 == 1:
|
|
return float(sorted_values[mid])
|
|
return float((sorted_values[mid - 1] + sorted_values[mid]) / 2.0)
|
|
|
|
|
|
def _generate_event_id() -> str:
|
|
now = datetime.now().strftime("%Y%m%d%H%M%S")
|
|
return f"LC-{now}-{uuid4().hex[:6]}"
|
|
|
|
|
|
def _publish_lightning_change(event_name: str, payload: dict[str, Any]) -> None:
|
|
_fire_and_forget(
|
|
publish_topic(
|
|
LIGHTNING_TOPIC,
|
|
name=event_name,
|
|
payload=payload,
|
|
requires_refetch=["/api/v1/lightning-currents"],
|
|
dedupe_key=f"{event_name}:{payload.get('event_id', 'unknown')}",
|
|
)
|
|
)
|
|
|
|
|
|
def _fire_and_forget(coro: object) -> None:
|
|
try:
|
|
loop = asyncio.get_running_loop()
|
|
except RuntimeError:
|
|
return
|
|
loop.create_task(coro)
|