wzy-warehouse
606 lines
22 KiB
Python
606 lines
22 KiB
Python
"""
|
||
降雨数据Service - 业务逻辑层
|
||
"""
|
||
import numpy as np
|
||
from scipy.spatial import Delaunay, ConvexHull
|
||
from typing import List, Dict, Any, Tuple, Optional
|
||
from datetime import datetime
|
||
|
||
from app.repositories.rainfall_repository import RainfallRepository
|
||
from app.utils.logger import setup_logging
|
||
|
||
logger = setup_logging()
|
||
|
||
|
||
class InterpolationService:
|
||
"""插值服务类"""
|
||
|
||
@staticmethod
|
||
def _create_buffer_points(
|
||
points_array: np.ndarray
|
||
) -> np.ndarray:
|
||
"""
|
||
创建缓冲点：在原始站点外围生成虚拟点以扩展插值区域
|
||
|
||
Args:
|
||
points_array: 原始站点坐标数组
|
||
|
||
Returns:
|
||
缓冲点坐标数组
|
||
"""
|
||
# 计算站点分布的中心
|
||
center = np.mean(points_array, axis=0)
|
||
|
||
# 计算站点到中心的最大距离
|
||
distances_from_center = np.sqrt(np.sum((points_array - center) ** 2, axis=1))
|
||
np.max(distances_from_center)
|
||
|
||
# 在站点外围生成缓冲点（沿着各个方向扩展）
|
||
buffer_points = []
|
||
num_angles = 360 # 每隔1度生成一个缓冲点
|
||
|
||
for angle_deg in range(0, 360, 360 // num_angles):
|
||
angle_rad = np.radians(angle_deg)
|
||
# 在凸包边界外扩展
|
||
for scale in [1.05, 1.1, 1.15]:
|
||
# 找到该方向上最远的站点
|
||
direction = np.array([np.cos(angle_rad), np.sin(angle_rad)])
|
||
projections = points_array @ direction
|
||
max_idx = np.argmax(projections)
|
||
|
||
# 在该方向上扩展
|
||
base_point = points_array[max_idx]
|
||
buffer_point = center + (base_point - center) * scale
|
||
buffer_points.append(buffer_point)
|
||
|
||
return np.array(buffer_points)
|
||
|
||
@staticmethod
|
||
def gaussian_smoothing(
|
||
grid_data: np.ndarray,
|
||
sigma: float = 1.5
|
||
) -> np.ndarray:
|
||
"""
|
||
高斯平滑滤波，减少边缘突变
|
||
|
||
Args:
|
||
grid_data: 栅格数据
|
||
sigma: 高斯核标准差
|
||
|
||
Returns:
|
||
平滑后的栅格数据
|
||
"""
|
||
from scipy.ndimage import gaussian_filter
|
||
|
||
# 只对有效数据进行平滑
|
||
valid_mask = ~np.isnan(grid_data)
|
||
if not np.any(valid_mask):
|
||
return grid_data
|
||
|
||
# 填充NaN值以便平滑
|
||
filled_data = grid_data.copy()
|
||
mean_val = np.nanmean(grid_data)
|
||
filled_data[~valid_mask] = mean_val
|
||
|
||
# 应用高斯滤波
|
||
smoothed = gaussian_filter(filled_data, sigma=sigma)
|
||
|
||
# 恢复原始NaN区域
|
||
result = np.where(valid_mask, smoothed, np.nan)
|
||
|
||
return result
|
||
|
||
@staticmethod
|
||
def calculate_adaptive_max_distance(
|
||
points_array: np.ndarray,
|
||
base_distance: float = 0.3,
|
||
min_distance: float = 0.15,
|
||
max_distance: float = 0.5
|
||
) -> float:
|
||
"""
|
||
根据站点密度自适应计算最大影响距离
|
||
|
||
Args:
|
||
points_array: 站点坐标数组
|
||
base_distance: 基础距离
|
||
min_distance: 最小距离
|
||
max_distance: 最大距离
|
||
|
||
Returns:
|
||
自适应的最大影响距离
|
||
"""
|
||
if len(points_array) < 3:
|
||
return base_distance
|
||
|
||
# 计算站点间的平均距离
|
||
from scipy.spatial import distance_matrix
|
||
dist_matrix = distance_matrix(points_array, points_array)
|
||
|
||
# 排除对角线（自身距离为0）
|
||
np.fill_diagonal(dist_matrix, np.inf)
|
||
avg_distance = np.mean(np.min(dist_matrix, axis=1))
|
||
|
||
# 根据平均距离调整max_distance
|
||
adaptive_distance = avg_distance * 3 # 约3倍平均站点间距
|
||
|
||
# 限制在合理范围内
|
||
return np.clip(adaptive_distance, min_distance, max_distance)
|
||
|
||
@staticmethod
|
||
def inverse_distance_weighting(
|
||
points: List[Tuple[float, float]],
|
||
values: List[float],
|
||
grid_lon: np.ndarray,
|
||
grid_lat: np.ndarray,
|
||
power: float = 2.0,
|
||
max_distance: float = None,
|
||
use_adaptive_distance: bool = True,
|
||
apply_smoothing: bool = True,
|
||
smoothing_sigma: float = 1.0
|
||
) -> np.ndarray:
|
||
"""
|
||
反距离权重插值 (IDW) - 优化版本
|
||
改进：
|
||
1. 高斯核衰减替代简单幂律
|
||
2. 自适应距离阈值
|
||
3. 边缘渐变处理
|
||
4. 高斯平滑减少突变
|
||
|
||
Args:
|
||
points: 已知点坐标 [(lon, lat), ...]
|
||
values: 已知点的值 [rainfall, ...]
|
||
grid_lon: 网格经度数组
|
||
grid_lat: 网格纬度数组
|
||
power: 距离幂次（基础值）
|
||
max_distance: 最大影响距离（度），None则自适应计算
|
||
use_adaptive_distance: 是否使用自适应距离
|
||
apply_smoothing: 是否应用平滑
|
||
smoothing_sigma: 平滑强度
|
||
|
||
Returns:
|
||
插值后的栅格数据，无效区域为 NaN
|
||
"""
|
||
points_array = np.array(points)
|
||
values_array = np.array(values)
|
||
|
||
# 创建网格
|
||
lon_grid, lat_grid = np.meshgrid(grid_lon, grid_lat)
|
||
result = np.full_like(lon_grid, np.nan)
|
||
|
||
# 自适应计算最大距离
|
||
if use_adaptive_distance or max_distance is None:
|
||
actual_max_distance = InterpolationService.calculate_adaptive_max_distance(
|
||
points_array
|
||
)
|
||
if max_distance is not None:
|
||
actual_max_distance = min(actual_max_distance, max_distance)
|
||
else:
|
||
actual_max_distance = max_distance
|
||
|
||
logger.info(f"使用最大影响距离: {actual_max_distance:.3f} 度")
|
||
|
||
# 计算站点的凸包（带边缘缓冲）
|
||
hull_mask = None
|
||
confidence_mask = None # 置信度掩码
|
||
if len(points_array) >= 3:
|
||
try:
|
||
# 创建缓冲站点：在原始站点外围添加虚拟点
|
||
buffer_points = InterpolationService._create_buffer_points(
|
||
points_array
|
||
)
|
||
|
||
# 合并原始站点和缓冲站点
|
||
all_points = np.vstack([points_array, buffer_points])
|
||
|
||
# 计算凸包
|
||
hull = ConvexHull(all_points)
|
||
hull_points = all_points[hull.vertices]
|
||
tri = Delaunay(hull_points)
|
||
|
||
# 向量化判断所有网格点是否在凸包内
|
||
grid_points = np.column_stack([lon_grid.ravel(), lat_grid.ravel()])
|
||
hull_indices = tri.find_simplex(grid_points)
|
||
hull_mask = hull_indices >= 0
|
||
hull_mask = hull_mask.reshape(lon_grid.shape)
|
||
|
||
# 计算置信度：基于到最近站点的距离
|
||
# 在凸包内但远离站点的区域降低置信度
|
||
from scipy.spatial import distance_matrix
|
||
grid_valid = grid_points[hull_mask.ravel()]
|
||
if len(grid_valid) > 0:
|
||
dist_to_stations = distance_matrix(grid_valid, points_array)
|
||
min_distances = np.min(dist_to_stations, axis=1)
|
||
|
||
# 创建置信度掩码（距离越远，置信度越低）
|
||
confidence = np.ones(len(grid_points))
|
||
confidence[hull_mask.ravel()] = np.exp(-min_distances / actual_max_distance)
|
||
confidence_mask = confidence.reshape(lon_grid.shape)
|
||
else:
|
||
confidence_mask = np.ones_like(lon_grid)
|
||
|
||
except Exception as e:
|
||
logger.warning(f"凸包计算失败: {e}，使用全区域插值")
|
||
hull_mask = np.ones_like(lon_grid, dtype=bool)
|
||
confidence_mask = np.ones_like(lon_grid)
|
||
else:
|
||
hull_mask = np.ones_like(lon_grid, dtype=bool)
|
||
confidence_mask = np.ones_like(lon_grid)
|
||
|
||
# 向量化计算所有网格点到所有站点的距离
|
||
lon_diff = lon_grid[:, :, np.newaxis] - points_array[np.newaxis, np.newaxis, :, 0]
|
||
lat_diff = lat_grid[:, :, np.newaxis] - points_array[np.newaxis, np.newaxis, :, 1]
|
||
distances = np.sqrt(lon_diff**2 + lat_diff**2)
|
||
|
||
# 过滤超出最大距离的站点
|
||
valid_mask = distances <= actual_max_distance
|
||
|
||
# 对于每个网格点，检查是否有有效站点
|
||
has_valid_stations = np.any(valid_mask, axis=2)
|
||
|
||
# 合并凸包掩码和有效站点掩码
|
||
final_mask = hull_mask & has_valid_stations
|
||
|
||
# 避免除零
|
||
distances = np.where(valid_mask, distances, np.inf)
|
||
distances = np.maximum(distances, 1e-10)
|
||
|
||
# 优化的权重计算：结合幂律和高斯衰减
|
||
# 近处使用幂律，远处使用高斯衰减使过渡更平滑
|
||
power_weights = 1.0 / (distances ** power)
|
||
gaussian_weights = np.exp(-0.5 * (distances / (actual_max_distance * 0.5)) ** 2)
|
||
|
||
# 混合权重：距离越远，高斯权重占比越大
|
||
distance_ratio = distances / actual_max_distance
|
||
mix_factor = np.clip(distance_ratio, 0, 1)
|
||
weights = (1 - mix_factor) * power_weights + mix_factor * gaussian_weights
|
||
|
||
weights = np.where(valid_mask, weights, 0)
|
||
|
||
# 加权平均
|
||
weighted_sum = np.sum(weights * values_array[np.newaxis, np.newaxis, :], axis=2)
|
||
weight_total = np.sum(weights, axis=2)
|
||
|
||
# 计算基础插值结果（使用 errstate 忽略预期的除零警告，np.where 已安全过滤）
|
||
with np.errstate(divide='ignore', invalid='ignore'):
|
||
result = np.where(
|
||
final_mask & (weight_total > 0),
|
||
weighted_sum / weight_total,
|
||
np.nan
|
||
)
|
||
|
||
# 应用置信度调整：边缘区域向邻近值渐变
|
||
if confidence_mask is not None:
|
||
# 计算全局平均降雨量作为边缘区域的基准
|
||
valid_rainfall = result[final_mask]
|
||
if len(valid_rainfall) > 0:
|
||
mean_rainfall = np.mean(valid_rainfall)
|
||
# 边缘区域向平均值渐变
|
||
result = np.where(
|
||
final_mask,
|
||
result,
|
||
np.nan
|
||
)
|
||
# 根据置信度调整结果，低置信度区域向均值靠拢
|
||
adjusted_result = result * confidence_mask + mean_rainfall * (1 - confidence_mask)
|
||
result = np.where(final_mask, adjusted_result, np.nan)
|
||
|
||
# 应用高斯平滑减少边缘突变
|
||
if apply_smoothing:
|
||
result = InterpolationService.gaussian_smoothing(result, sigma=smoothing_sigma)
|
||
|
||
return result
|
||
|
||
@staticmethod
|
||
def get_rainfall_color(rainfall: float, duration: int = 12) -> str:
|
||
"""
|
||
根据降雨量获取颜色（按照国标）
|
||
|
||
Args:
|
||
rainfall: 降雨量(mm)
|
||
duration: 持续时间（12或24小时）
|
||
|
||
Returns:
|
||
颜色字符串 "rgba(r,g,b,a)"
|
||
"""
|
||
# 国标降雨等级颜色映射
|
||
if rainfall < 0.1:
|
||
return "rgba(200,200,200,0)" # 透明 - 微量降雨（零星小雨）
|
||
elif rainfall < 10 if duration == 12 else 9.9:
|
||
return "rgba(0,0,255,0.4)" # 浅蓝 - 小雨
|
||
elif rainfall < 15 if duration == 12 else 25:
|
||
return "rgba(0,255,255,0.5)" # 青色 - 中雨
|
||
elif rainfall < 30 if duration == 12 else 50:
|
||
return "rgba(0,255,0,0.6)" # 绿色 - 大雨
|
||
elif rainfall < 70 if duration == 12 else 100:
|
||
return "rgba(255,255,0,0.7)" # 黄色 - 暴雨
|
||
elif rainfall < 140 if duration == 12 else 250:
|
||
return "rgba(255,165,0,0.8)" # 橙色 - 大暴雨
|
||
else:
|
||
return "rgba(255,0,0,0.9)" # 红色 - 特大暴雨
|
||
|
||
|
||
class GeoJSONService:
|
||
"""GeoJSON生成服务"""
|
||
|
||
@staticmethod
|
||
def create_feature_collection(
|
||
grid_metadata: Dict[str, Any],
|
||
rainfall_array: np.ndarray,
|
||
grid_lon: np.ndarray,
|
||
grid_lat: np.ndarray,
|
||
duration: int = 12
|
||
) -> Dict[str, Any]:
|
||
"""
|
||
创建GeoJSON FeatureCollection用于Cesium渲染
|
||
|
||
Args:
|
||
grid_metadata: 栅格元数据
|
||
rainfall_array: 降雨量数组
|
||
grid_lon: 经度网格
|
||
grid_lat: 纬度网格
|
||
duration: 持续时间（12或24小时）
|
||
|
||
Returns:
|
||
GeoJSON格式的FeatureCollection
|
||
"""
|
||
features = []
|
||
|
||
# 将栅格数据转换为矩形要素
|
||
for i in range(len(grid_lat) - 1):
|
||
for j in range(len(grid_lon) - 1):
|
||
rainfall_value = float(rainfall_array[i, j])
|
||
|
||
# 跳过无数据的区域
|
||
if np.isnan(rainfall_value) or rainfall_value < 0:
|
||
continue
|
||
|
||
# 创建矩形多边形
|
||
lon_min = float(grid_lon[j])
|
||
lon_max = float(grid_lon[j + 1])
|
||
lat_min = float(grid_lat[i])
|
||
lat_max = float(grid_lat[i + 1])
|
||
|
||
feature = {
|
||
"type": "Feature",
|
||
"geometry": {
|
||
"type": "Polygon",
|
||
"coordinates": [[
|
||
[lon_min, lat_min],
|
||
[lon_max, lat_min],
|
||
[lon_max, lat_max],
|
||
[lon_min, lat_max],
|
||
[lon_min, lat_min]
|
||
]]
|
||
},
|
||
"properties": {
|
||
"rainfall": round(rainfall_value, 2),
|
||
"level": RainfallService._get_rainfall_level(rainfall_value, duration),
|
||
"color": InterpolationService.get_rainfall_color(rainfall_value, duration)
|
||
}
|
||
}
|
||
features.append(feature)
|
||
|
||
return {
|
||
"type": "FeatureCollection",
|
||
"features": features,
|
||
"metadata": {
|
||
"resolution": grid_metadata['resolution'],
|
||
"grid_size": [len(grid_lon), len(grid_lat)],
|
||
"bounds": {
|
||
"min_lon": float(grid_lon.min()),
|
||
"max_lon": float(grid_lon.max()),
|
||
"min_lat": float(grid_lat.min()),
|
||
"max_lat": float(grid_lat.max())
|
||
}
|
||
}
|
||
}
|
||
|
||
|
||
class RainfallService:
|
||
"""降雨数据业务服务类"""
|
||
|
||
def __init__(self):
|
||
self.repository = RainfallRepository()
|
||
self.interpolation_service = InterpolationService()
|
||
self.geojson_service = GeoJSONService()
|
||
|
||
def get_stations_data(
|
||
self,
|
||
query_time: datetime,
|
||
duration: int = 12
|
||
) -> List[Dict[str, Any]]:
|
||
"""
|
||
获取站点降雨数据
|
||
|
||
Args:
|
||
query_time: 查询时间（自动查询前12小时或24小时数据）
|
||
duration: 持续时间（12或24小时）
|
||
|
||
Returns:
|
||
站点数据列表
|
||
"""
|
||
return self.repository.query_stations_rainfall(query_time, duration)
|
||
|
||
def generate_rainfall_grid(
|
||
self,
|
||
query_time: datetime,
|
||
resolution: float = 0.01,
|
||
duration: int = 12
|
||
) -> Dict[str, Any]:
|
||
"""
|
||
生成降雨栅格数据
|
||
|
||
Args:
|
||
query_time: 查询时间（自动查询前12小时或24小时数据）
|
||
resolution: 栅格分辨率
|
||
duration: 持续时间（12或24小时）
|
||
|
||
Returns:
|
||
GeoJSON格式的栅格数据
|
||
"""
|
||
logger.info(f"查询降雨数据: {query_time}, 持续时间: {duration}小时")
|
||
|
||
# 查询站点数据（自动查询前12小时或24小时数据）
|
||
stations_data = self.get_stations_data(query_time, duration)
|
||
|
||
if not stations_data:
|
||
return None
|
||
|
||
# 提取站点坐标和降雨量（过滤空值）
|
||
valid_stations = [row for row in stations_data if row['rainfall'] is not None]
|
||
|
||
if not valid_stations:
|
||
logger.warning("所有站点的降雨量数据均为空")
|
||
return None
|
||
|
||
points = [(row['lon'], row['lat']) for row in valid_stations]
|
||
values = [float(row['rainfall']) for row in valid_stations]
|
||
|
||
# 确定栅格范围
|
||
lon_min, lon_max = 107, 110
|
||
lat_min, lat_max = 33, 35
|
||
|
||
# 创建栅格网格
|
||
num_lon = int((lon_max - lon_min) / resolution) + 1
|
||
num_lat = int((lat_max - lat_min) / resolution) + 1
|
||
|
||
grid_lon = np.linspace(lon_min, lon_max, num_lon)
|
||
grid_lat = np.linspace(lat_min, lat_max, num_lat)
|
||
|
||
logger.info(f"生成栅格: {num_lon}x{num_lat}, 分辨率: {resolution}")
|
||
|
||
# 执行IDW插值（优化版本：自适应距离、混合权重、平滑处理）
|
||
rainfall_grid = self.interpolation_service.inverse_distance_weighting(
|
||
points=points,
|
||
values=values,
|
||
grid_lon=grid_lon,
|
||
grid_lat=grid_lat,
|
||
power=2.0,
|
||
max_distance=0.35, # 最大影响距离0.35度（约35公里）
|
||
use_adaptive_distance=True, # 启用自适应距离
|
||
apply_smoothing=True, # 启用平滑处理
|
||
smoothing_sigma=1.2 # 平滑强度
|
||
)
|
||
|
||
# 创建栅格元数据
|
||
grid_metadata = {
|
||
"query_time": query_time.isoformat(),
|
||
"resolution": resolution,
|
||
"station_count": len(stations_data),
|
||
"grid_size": [num_lon, num_lat]
|
||
}
|
||
|
||
# 转换为GeoJSON格式
|
||
geojson_data = self.geojson_service.create_feature_collection(
|
||
grid_metadata, rainfall_grid, grid_lon, grid_lat, duration
|
||
)
|
||
|
||
logger.info("降雨栅格数据生成成功")
|
||
|
||
return geojson_data
|
||
|
||
def get_rainfall_at_point(
|
||
self,
|
||
longitude: float,
|
||
latitude: float,
|
||
query_time: datetime,
|
||
duration: int = 12
|
||
) -> Optional[Dict[str, Any]]:
|
||
"""
|
||
查询指定点位的降雨量（使用IDW插值）
|
||
|
||
Args:
|
||
longitude: 经度
|
||
latitude: 纬度
|
||
query_time: 查询时间（自动查询前12小时或24小时数据）
|
||
duration: 持续时间（12或24小时）
|
||
|
||
Returns:
|
||
点位降雨量信息
|
||
"""
|
||
# 获取站点数据（自动查询前12小时或24小时数据）
|
||
stations_data = self.get_stations_data(query_time, duration)
|
||
|
||
if not stations_data:
|
||
return None
|
||
|
||
# 提取站点坐标和降雨量
|
||
points = [(row['lon'], row['lat']) for row in stations_data]
|
||
values = [float(row['rainfall']) for row in stations_data]
|
||
|
||
# 使用IDW插值计算该点的降雨量
|
||
target_point = np.array([[longitude, latitude]])
|
||
points_array = np.array(points)
|
||
|
||
# 计算距离
|
||
distances = np.sqrt(
|
||
(points_array[:, 0] - longitude) ** 2 +
|
||
(points_array[:, 1] - latitude) ** 2
|
||
)
|
||
|
||
# 避免除零
|
||
min_dist = 1e-10
|
||
distances = np.maximum(distances, min_dist)
|
||
|
||
# IDW公式
|
||
power = 2.0
|
||
weights = 1.0 / (distances ** power)
|
||
rainfall_value = np.sum(weights * values) / np.sum(weights)
|
||
|
||
# 返回结果
|
||
return {
|
||
"longitude": longitude,
|
||
"latitude": latitude,
|
||
"rainfall": round(float(rainfall_value), 2),
|
||
"level": self._get_rainfall_level(rainfall_value, duration),
|
||
"color": InterpolationService.get_rainfall_color(rainfall_value, duration),
|
||
"station_count": len(stations_data),
|
||
"query_time": query_time.isoformat(),
|
||
"duration": duration
|
||
}
|
||
|
||
@staticmethod
|
||
def _get_rainfall_level(rainfall: float, duration: int = 12) -> str:
|
||
"""
|
||
获取降雨等级（按照国标）
|
||
|
||
Args:
|
||
rainfall: 降雨量(mm)
|
||
duration: 持续时间（12或24小时）
|
||
|
||
Returns:
|
||
降雨等级字符串
|
||
"""
|
||
if duration == 12:
|
||
# 12小时降雨等级标准
|
||
if rainfall < 0.1:
|
||
return "微量降雨"
|
||
elif rainfall < 5.0:
|
||
return "小雨"
|
||
elif rainfall < 15.0:
|
||
return "中雨"
|
||
elif rainfall < 30.0:
|
||
return "大雨"
|
||
elif rainfall < 70.0:
|
||
return "暴雨"
|
||
elif rainfall < 140.0:
|
||
return "大暴雨"
|
||
else:
|
||
return "特大暴雨"
|
||
else: # 24小时
|
||
# 24小时降雨等级标准
|
||
if rainfall < 0.1:
|
||
return "微量降雨"
|
||
elif rainfall < 10.0:
|
||
return "小雨"
|
||
elif rainfall < 25.0:
|
||
return "中雨"
|
||
elif rainfall < 50.0:
|
||
return "大雨"
|
||
elif rainfall < 100.0:
|
||
return "暴雨"
|
||
elif rainfall < 250.0:
|
||
return "大暴雨"
|
||
else:
|
||
return "特大暴雨"
|