修改分辨率触发逻辑

app/services/rainfall_grid_service.py

@@ -1,5 +1,5 @@
 """
-降雨栅格服务
+降雨栅格服务 - 内存优化版本
 负责降雨插值、边缘优化、PNG生成等业务逻辑
 """
 import os
@@ -42,7 +42,7 @@ class RainfallGridService:
         }
 
         # 栅格分辨率（度）
-        self.grid_resolution = 0.01  # 约1km
+        self.grid_resolution = 0.001
 
     def _create_buffer_points(self, points_array) -> 'np.ndarray':
         """
@@ -119,16 +119,13 @@ class RainfallGridService:
 
     def interpolate_rainfall(self, station_data: List[Dict[str, Any]]) -> Dict[str, Any]:
         """
-        使用优化的反距离权重法（IDW）进行降雨插值
+        使用优化的反距离权重法（IDW）进行降雨插值（内存优化版本）
 
-        注意：station_data 现在包含 'rainfall'（累计降雨量）和 'duration_hours'（持续时间）
+        内存优化：
-        与DBN推演使用相同的降雨量计算逻辑（72小时回溯 + 3小时无雨截断）
+        1. 使用float32代替float64（内存减半）
-        
+        2. 分块处理距离计算
-        改进：
+        3. 提前过滤无效站点
-        1. 高斯核衰减替代简单幂律
+        4. 减少中间数组
-        2. 自适应距离阈值
-        3. 边缘渐变处理
-        4. 高斯平滑减少突变
 
         Args:
             station_data: 站点数据列表，格式：
@@ -141,27 +138,29 @@ class RainfallGridService:
             插值结果字典
         """
         import numpy as np
-        from scipy.spatial import Delaunay, ConvexHull, distance_matrix
+        from scipy.spatial import Delaunay, ConvexHull
         from scipy.ndimage import gaussian_filter
 
         # 提取站点坐标和降雨量
-        points_array = np.array([[s['lon'], s['lat']] for s in station_data])
+        points_array = np.array([[s['lon'], s['lat']] for s in station_data], dtype=np.float32)
-        values_array = np.array([s['rainfall'] for s in station_data])
+        values_array = np.array([s['rainfall'] for s in station_data], dtype=np.float32)
 
         # 创建栅格网格
         lon_range = np.arange(
             self.xian_bounds['min_lon'],
             self.xian_bounds['max_lon'],
-            self.grid_resolution
+            self.grid_resolution,
+            dtype=np.float32
         )
         lat_range = np.arange(
             self.xian_bounds['min_lat'],
             self.xian_bounds['max_lat'],
-            self.grid_resolution
+            self.grid_resolution,
+            dtype=np.float32
         )
 
         grid_lon, grid_lat = np.meshgrid(lon_range, lat_range)
-        result = np.full_like(grid_lon, np.nan)
+        result = np.full_like(grid_lon, np.nan, dtype=np.float32)
 
         # 自适应计算最大距离
         actual_max_distance = self._calculate_adaptive_max_distance(points_array)
@@ -185,73 +184,160 @@ class RainfallGridService:
 
                 # 向量化判断所有网格点是否在凸包内
                 grid_points = np.column_stack([grid_lon.ravel(), grid_lat.ravel()])
-                hull_indices = tri.find_simplex(grid_points)
+                simplex_indices = tri.find_simplex(grid_points)
-                hull_mask = hull_indices >= 0
+                hull_mask = simplex_indices >= 0
                 hull_mask = hull_mask.reshape(grid_lon.shape)
 
-                # 计算置信度：基于到最近站点的距离
+                # 计算置信度：基于到最近站点的距离（分块处理）
                 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)
+                    chunk_size = 100000  # 每次处理10万点
+                    n_valid = len(grid_valid)
+                    min_distances = np.zeros(n_valid, dtype=np.float32)
+                    
+                    for i in range(0, n_valid, chunk_size):
+                        chunk_end = min(i + chunk_size, n_valid)
+                        chunk_points = grid_valid[i:chunk_end]
+                        
+                        # 计算当前块到所有站点的距离
+                        lon_diff = chunk_points[:, 0:1] - points_array[np.newaxis, :, 0]
+                        lat_diff = chunk_points[:, 1:2] - points_array[np.newaxis, :, 1]
+                        distances = np.sqrt(lon_diff**2 + lat_diff**2)
+                        
+                        # 记录最小距离
+                        min_distances[i:chunk_end] = np.min(distances, axis=1)
+                        
+                        # 释放临时数组
+                        del lon_diff, lat_diff, distances
 
                     # 创建置信度掩码（距离越远，置信度越低）
-                    confidence = np.ones(len(grid_points))
+                    confidence = np.ones(len(grid_points), dtype=np.float32)
                     confidence[hull_mask.ravel()] = np.exp(-min_distances / actual_max_distance)
                     confidence_mask = confidence.reshape(grid_lon.shape)
                 else:
-                    confidence_mask = np.ones_like(grid_lon)
+                    confidence_mask = np.ones_like(grid_lon, dtype=np.float32)
 
             except Exception as e:
                 self.logger.warning(f"凸包计算失败: {e}，使用全区域插值")
                 hull_mask = np.ones_like(grid_lon, dtype=bool)
-                confidence_mask = np.ones_like(grid_lon)
+                confidence_mask = np.ones_like(grid_lon, dtype=np.float32)
         else:
             hull_mask = np.ones_like(grid_lon, dtype=bool)
-            confidence_mask = np.ones_like(grid_lon)
+            confidence_mask = np.ones_like(grid_lon, dtype=np.float32)
 
-        # 向量化计算所有网格点到所有站点的距离
+        # 获取凸包内网格点坐标
-        lon_diff = grid_lon[:, :, np.newaxis] - points_array[np.newaxis, np.newaxis, :, 0]
+        grid_points = np.column_stack([grid_lon.ravel(), grid_lat.ravel()])
-        lat_diff = grid_lat[:, :, np.newaxis] - points_array[np.newaxis, np.newaxis, :, 1]
+        if hull_mask is not None:
-        distances = np.sqrt(lon_diff**2 + lat_diff**2)
+            # 只计算凸包内网格点到站点的距离
+            hull_point_indices = np.where(hull_mask.ravel())[0]
+            grid_points_hull = grid_points[hull_point_indices]
+            n_hull_points = len(grid_points_hull)
+            self.logger.info(f"凸包内网格点数量: {n_hull_points}, 总网格点: {grid_lon.size}")
+        else:
+            # 如果凸包掩码不可用，使用所有网格点
+            grid_points_hull = grid_points
+            hull_point_indices = np.arange(len(grid_points))
+            n_hull_points = len(grid_points_hull)
+
+        # 分块计算凸包内网格点到所有站点的距离
+        chunk_size = 50000  # 每次处理5万点
+        result_hull = np.full(n_hull_points, np.nan, dtype=np.float32)
+        
+        # 用于记录最后一个块的has_valid_stations_chunk
+        last_has_valid_stations_chunk = None
+        
+        for i in range(0, n_hull_points, chunk_size):
+            chunk_end = min(i + chunk_size, n_hull_points)
+            chunk_points = grid_points_hull[i:chunk_end]
+            chunk_size_actual = chunk_end - i
+            
+            # 计算当前块到所有站点的距离
+            lon_diff_chunk = chunk_points[:, 0:1] - points_array[np.newaxis, :, 0]
+            lat_diff_chunk = chunk_points[:, 1:2] - points_array[np.newaxis, :, 1]
+            distances_chunk = np.sqrt(lon_diff_chunk**2 + lat_diff_chunk**2)
             
             # 过滤超出最大距离的站点
-        valid_mask = distances <= actual_max_distance
+            valid_mask_chunk = distances_chunk <= actual_max_distance
             
             # 对于每个网格点，检查是否有有效站点
-        has_valid_stations = np.any(valid_mask, axis=2)
+            has_valid_stations_chunk = np.any(valid_mask_chunk, axis=1)
-        
-        # 合并凸包掩码和有效站点掩码
-        final_mask = hull_mask & has_valid_stations
             
             # 避免除零
-        distances = np.where(valid_mask, distances, np.inf)
+            distances_chunk = np.where(valid_mask_chunk, distances_chunk, np.inf)
-        distances = np.maximum(distances, 1e-10)
+            distances_chunk = np.maximum(distances_chunk, 1e-10)
             
             # 优化的权重计算：结合幂律和高斯衰减
             power = 2.0
-        power_weights = 1.0 / (distances ** power)
+            power_weights_chunk = 1.0 / (distances_chunk ** power)
-        gaussian_weights = np.exp(-0.5 * (distances / (actual_max_distance * 0.5)) ** 2)
+            gaussian_weights_chunk = np.exp(-0.5 * (distances_chunk / (actual_max_distance * 0.5)) ** 2)
             
             # 混合权重：距离越远，高斯权重占比越大
-        distance_ratio = distances / actual_max_distance
+            distance_ratio_chunk = distances_chunk / actual_max_distance
-        mix_factor = np.clip(distance_ratio, 0, 1)
+            mix_factor_chunk = np.clip(distance_ratio_chunk, 0, 1)
-        weights = (1 - mix_factor) * power_weights + mix_factor * gaussian_weights
+            weights_chunk = (1 - mix_factor_chunk) * power_weights_chunk + mix_factor_chunk * gaussian_weights_chunk
             
-        weights = np.where(valid_mask, weights, 0)
+            weights_chunk = np.where(valid_mask_chunk, weights_chunk, 0)
             
             # 加权平均
-        weighted_sum = np.sum(weights * values_array[np.newaxis, np.newaxis, :], axis=2)
+            weighted_sum_chunk = np.sum(weights_chunk * values_array[np.newaxis, :], axis=1)
-        weight_total = np.sum(weights, axis=2)
+            weight_total_chunk = np.sum(weights_chunk, axis=1)
             
-        # 计算基础插值结果
+            # 计算当前块的插值结果
             with np.errstate(divide='ignore', invalid='ignore'):
-            result = np.where(
+                chunk_result = np.where(
-                final_mask & (weight_total > 0),
+                    has_valid_stations_chunk & (weight_total_chunk > 0),
-                weighted_sum / weight_total,
+                    weighted_sum_chunk / weight_total_chunk,
                     np.nan
                 )
             
+            # 存储结果
+            result_hull[i:chunk_end] = chunk_result
+            
+            # 记录最后一个块的has_valid_stations_chunk
+            last_has_valid_stations_chunk = has_valid_stations_chunk
+            
+            # 释放临时数组
+            del lon_diff_chunk, lat_diff_chunk, distances_chunk, valid_mask_chunk
+            del power_weights_chunk, gaussian_weights_chunk, weights_chunk
+            del weighted_sum_chunk, weight_total_chunk, chunk_result
+
+        # 将凸包内点的结果映射回完整网格
+        result = np.full_like(grid_lon, np.nan, dtype=np.float32)
+        result.ravel()[hull_point_indices] = result_hull
+
+        # 构建完整网格的有效掩码（凸包内且有有效站点）
+        # 注意：这里需要重新计算所有凸包内点的有效站点掩码
+        # 由于分块处理，我们需要重新计算完整掩码
+        has_valid_stations_full = np.zeros_like(grid_lon, dtype=bool)
+        
+        # 重新计算所有凸包内点的有效站点掩码（内存优化：分块计算）
+        if n_hull_points > 0:
+            # 分块计算有效站点掩码
+            chunk_size_mask = 100000  # 每次处理10万点
+            for i in range(0, n_hull_points, chunk_size_mask):
+                chunk_end = min(i + chunk_size_mask, n_hull_points)
+                chunk_points = grid_points_hull[i:chunk_end]
+                
+                # 计算当前块到所有站点的距离
+                lon_diff_chunk = chunk_points[:, 0:1] - points_array[np.newaxis, :, 0]
+                lat_diff_chunk = chunk_points[:, 1:2] - points_array[np.newaxis, :, 1]
+                distances_chunk = np.sqrt(lon_diff_chunk**2 + lat_diff_chunk**2)
+                
+                # 过滤超出最大距离的站点
+                valid_mask_chunk = distances_chunk <= actual_max_distance
+                
+                # 对于每个网格点，检查是否有有效站点
+                has_valid_chunk = np.any(valid_mask_chunk, axis=1)
+                
+                # 存储结果
+                has_valid_stations_full.ravel()[hull_point_indices[i:chunk_end]] = has_valid_chunk
+                
+                # 释放临时数组
+                del lon_diff_chunk, lat_diff_chunk, distances_chunk, valid_mask_chunk, has_valid_chunk
+        
+        final_mask = hull_mask & has_valid_stations_full
+
         # 应用置信度调整：边缘区域向邻近值渐变
         if confidence_mask is not None:
             valid_rainfall = result[final_mask]
@@ -310,8 +396,8 @@ class RainfallGridService:
 
         try:
             grid_values = grid_data['grid_values']
-            lon_range = grid_data['lon_range']
+            lon_range = grid_data['grid_lon']
-            lat_range = grid_data['lat_range']
+            lat_range = grid_data['grid_lat']
 
             # 创建自定义颜色映射
             levels = self.rainfall_levels['levels']