Przeglądaj źródła

优化土地数据计算接口,支持动态边界和插值控制,更新日志记录信息,确保参数验证和错误处理更加完善。重构数据处理流程,直接从数据库获取边界数据

drggboy 1 tydzień temu
rodzic
commit
c27c0e0410
2 zmienionych plików z 201 dodań i 57 usunięć
  1. 30 11
      app/api/water.py
  2. 171 46
      app/services/water_service.py

+ 30 - 11
app/api/water.py

@@ -16,13 +16,17 @@ from ..services.water_service import (
     get_base_dir
 )
 
-# 配置日志
+# 配置日志 - 避免重复日志输出
 logger = logging.getLogger(__name__)
-logger.setLevel(logging.INFO)
-handler = logging.StreamHandler()
-formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s')
-handler.setFormatter(formatter)
-logger.addHandler(handler)
+# 避免重复添加处理器导致重复日志
+if not logger.handlers:
+    logger.setLevel(logging.INFO)
+    handler = logging.StreamHandler()
+    formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s')
+    handler.setFormatter(formatter)
+    logger.addHandler(handler)
+    # 关闭日志传播,避免与父级日志处理器重复输出
+    logger.propagate = False
 
 router = APIRouter()
 
@@ -110,18 +114,28 @@ async def get_default_histogram(
 
 @router.post("/calculate",
              summary="重新计算土地数据",
-             description="根据输入的土地类型和系数重新计算土地数据")
+             description="根据输入的土地类型和系数重新计算土地数据,支持通过area和level参数控制地图边界,支持插值控制")
 async def recalculate_land_data(
         background_tasks: BackgroundTasks,
         land_type: str = Form(..., description="土地类型,如:'水田'、'旱地'或'水浇地'"),
         param1: float = Form(711, description="土地类型系数的第一个参数"),
         param2: float = Form(0.524, description="土地类型系数的第二个参数"),
         color_map_name: str = Form("绿-黄-红-紫", description="使用的色彩方案"),
-        output_size: int = Form(8, description="输出图片的尺寸")
+        output_size: int = Form(8, description="输出图片的尺寸"),
+        area: Optional[str] = Form(None, description="可选的地区名称,如:'乐昌市',用于动态控制地图边界"),
+        level: Optional[str] = Form(None, description="可选的行政层级,必须为: county | city | province"),
+        enable_interpolation: Optional[bool] = Form(False, description="是否启用空间插值,默认启用"),
+        interpolation_method: Optional[str] = Form("linear", description="插值方法: nearest | linear | cubic"),
+        resolution_factor: Optional[float] = Form(4.0, description="分辨率因子,默认4.0,越大分辨率越高")
 ) -> Dict[str, Any]:
-    """重新计算土地数据并返回结果路径"""
+    """重新计算土地数据并返回结果路径,支持动态边界控制和插值控制"""
     try:
         logger.info(f"重新计算土地数据: {land_type}")
+        if area and level:
+            logger.info(f"使用动态边界: {area} ({level})")
+        else:
+            logger.info("使用默认边界")
+        logger.info(f"插值设置: 启用={enable_interpolation}, 方法={interpolation_method}, 分辨率因子={resolution_factor}")
 
         # 获取默认目录
         base_dir = get_base_dir()
@@ -135,12 +149,17 @@ async def recalculate_land_data(
             land_type: (param1, param2)
         }
 
-        # 调用服务层函数处理土地数据
+        # 调用服务层函数处理土地数据(支持动态边界和插值控制)
         results = process_land_to_visualization(
             land_type=land_type,
             coefficient_params=coefficient_params,
             color_map_name=color_map_name,
-            output_size=output_size
+            output_size=output_size,
+            area=area,
+            level=level,
+            enable_interpolation=enable_interpolation,
+            interpolation_method=interpolation_method,
+            resolution_factor=resolution_factor
         )
 
         if not results or any(r is None for r in results):

+ 171 - 46
app/services/water_service.py

@@ -15,6 +15,11 @@ import sys
 from ..utils.mapping_utils import MappingUtils, csv_to_raster_workflow
 # 导入土地利用服务
 from .land_use_service import land_use_service
+# 导入边界服务
+from .admin_boundary_service import get_boundary_gdf_by_name
+from ..database import SessionLocal
+import geopandas as gpd
+import tempfile
 
 # 配置日志
 logger = logging.getLogger(__name__)
@@ -58,6 +63,32 @@ def get_base_dir():
         return os.path.dirname(os.path.abspath(__file__))
 
 
+def get_boundary_gdf_from_database(area: str, level: str) -> Optional[gpd.GeoDataFrame]:
+    """
+    直接从数据库获取边界数据作为GeoDataFrame
+    
+    @description: 与其他模块保持一致的边界获取方法
+    @param area: 地区名称
+    @param level: 行政层级(county | city | province)
+    @returns: 边界GeoDataFrame或None
+    """
+    try:
+        with SessionLocal() as db:
+            norm_area = area.strip()
+            boundary_gdf = get_boundary_gdf_by_name(db, norm_area, level=level)
+            if boundary_gdf is not None:
+                logger.info(f"从数据库获取边界数据: {norm_area} ({level})")
+            return boundary_gdf
+                
+    except Exception as e:
+        logger.warning(f"从数据库获取边界数据失败: {str(e)}")
+        
+    return None
+
+
+
+
+
 # 土地数据处理函数
 def process_land_data(land_type, coefficient_params=None):
     """处理土地类型数据并生成清洗后的简化数据"""
@@ -161,9 +192,20 @@ def process_land_data(land_type, coefficient_params=None):
     return cleaned_csv, Num
 
 
-# 可视化函数(使用MappingUtils)
+# 可视化函数(完全使用统一的MappingUtils接口
 def mapping_raster(shp_path, tif_path, color_map_name, title_name, output_path, output_size=8):
-    """生成栅格地图可视化"""
+    """
+    生成栅格地图可视化
+    
+    @description: 使用统一的MappingUtils接口生成栅格地图
+    @param shp_path: 边界shapefile路径
+    @param tif_path: 栅格数据文件路径  
+    @param color_map_name: 色彩方案名称(中文)
+    @param title_name: 地图标题
+    @param output_path: 输出图片路径
+    @param output_size: 图片尺寸
+    @returns: 输出图片路径或None(失败时)
+    """
     try:
         logger.info(f"生成栅格地图: {title_name}")
 
@@ -176,11 +218,11 @@ def mapping_raster(shp_path, tif_path, color_map_name, title_name, output_path,
         # 转换颜色方案名称
         colormap_key = COLORMAP_MAPPING.get(color_map_name, "green_yellow_red_purple")
 
-        # 调用MappingUtils的高级绘图功能
-        mapper.create_raster_map(
+        # 直接调用统一的绘图接口,无需额外包装
+        generated_path = mapper.create_raster_map(
             shp_path=shp_path,
             tif_path=tif_path,
-            output_path=os.path.splitext(output_path)[0],  # 去掉扩展名
+            output_path=os.path.splitext(output_path)[0],  # 去掉扩展名,MappingUtils会自动添加.jpg
             colormap=colormap_key,
             title=title_name,
             output_size=output_size,
@@ -190,13 +232,13 @@ def mapping_raster(shp_path, tif_path, color_map_name, title_name, output_path,
             interpolation_method='linear'
         )
 
-        # MappingUtils会自动添加.jpg扩展名,重命名文件以匹配原始输出路径
-        generated_path = os.path.splitext(output_path)[0] + ".jpg"
+        # 如果生成的路径与期望路径不同,则重命名
         if generated_path != output_path and os.path.exists(generated_path):
             shutil.move(generated_path, output_path)
             logger.info(f"重命名图片: {generated_path} -> {output_path}")
-
-        return output_path
+            return output_path
+        
+        return generated_path
 
     except Exception as e:
         logger.error(f"栅格地图生成失败: {str(e)}")
@@ -206,15 +248,26 @@ def mapping_raster(shp_path, tif_path, color_map_name, title_name, output_path,
 def plot_tif_histogram(file_path, output_path, figsize=(8, 8),
                        xlabel='Cd(ug/L)', ylabel='frequency density',
                        title='Irrigation water input flux'):
-    """生成TIFF文件的直方图"""
+    """
+    生成TIFF文件的直方图
+    
+    @description: 使用统一的MappingUtils接口生成直方图
+    @param file_path: TIFF文件路径
+    @param output_path: 输出图片路径
+    @param figsize: 图片尺寸元组
+    @param xlabel: X轴标签
+    @param ylabel: Y轴标签  
+    @param title: 图表标题
+    @returns: 输出图片路径或None(失败时)
+    """
     try:
         logger.info(f"生成直方图: {file_path}")
 
         # 创建MappingUtils实例
         mapper = MappingUtils()
 
-        # 调用MappingUtils的高级直方图功能
-        mapper.create_histogram(
+        # 直接调用统一的绘图接口,无需额外包装
+        generated_path = mapper.create_histogram(
             file_path=file_path,
             save_path=output_path,
             figsize=figsize,
@@ -225,25 +278,45 @@ def plot_tif_histogram(file_path, output_path, figsize=(8, 8),
             dpi=300
         )
 
-        return output_path
+        return generated_path
 
     except Exception as e:
         logger.error(f"直方图生成失败: {str(e)}")
         return None
 
 
-# 完整的处理流程
+# 完整的处理流程(已更新为完全使用统一绘图接口,支持动态边界和插值控制)
 def process_land_to_visualization(land_type, coefficient_params=None,
                                   color_map_name="绿-黄-红-紫",
-                                  output_size=8):
+                                  output_size=8,
+                                  area: Optional[str] = None,
+                                  level: Optional[str] = None,
+                                  enable_interpolation: Optional[bool] = True,
+                                  interpolation_method: Optional[str] = "linear",
+                                  resolution_factor: Optional[float] = 4.0):
     """
-    完整的土地数据处理可视化流程:
-    1. 生成清洗后CSV
-    2. 使用csv_to_raster_workflow转换为GeoTIFF
-    3. 生成栅格地图
-    4. 生成直方图
-
-    返回所有生成文件的路径
+    完整的土地数据处理可视化流程(使用统一的MappingUtils接口,支持动态边界和插值控制)
+    
+    @description: 水样数据处理与可视化的完整工作流程,已优化为完全使用统一的绘图接口,
+                  支持通过area和level参数动态控制地图边界,支持插值控制参数
+    
+    工作流程:
+    1. 生成清洗后CSV - 从数据库获取土地利用数据并计算Cd含量
+    2. 获取动态边界 - 根据area和level参数从数据库获取边界数据
+    3. 使用csv_to_raster_workflow转换为GeoTIFF - 调用统一的栅格转换工具,支持插值控制
+    4. 生成栅格地图 - 直接调用MappingUtils.create_raster_map(),支持动态边界
+    5. 生成直方图 - 直接调用MappingUtils.create_histogram()
+
+    @param land_type: 土地类型(水田、旱地、水浇地)
+    @param coefficient_params: 可选的系数参数字典
+    @param color_map_name: 色彩方案名称(中文)
+    @param output_size: 输出图片尺寸
+    @param area: 可选的地区名称,如"乐昌市",用于动态控制地图边界
+    @param level: 可选的行政层级,county | city | province
+    @param enable_interpolation: 是否启用空间插值,默认True
+    @param interpolation_method: 插值方法,nearest | linear | cubic,默认linear
+    @param resolution_factor: 分辨率因子,默认4.0,越大分辨率越高
+    @returns: 包含所有生成文件路径的元组
     """
     base_dir = get_base_dir()
     logger.info(f"开始处理: {land_type}")
@@ -254,51 +327,103 @@ def process_land_to_visualization(land_type, coefficient_params=None,
         logger.error(f"处理土地数据失败: {land_type}")
         return None, None, None, None, None, None
 
-    # 2. 使用csv_to_raster_workflow转换为GeoTIFF
+    # 2. 获取动态边界数据(提前获取用于栅格转换)
+    boundary_gdf = None
+    if area and level:
+        logger.info(f"尝试获取动态边界: area='{area}', level='{level}'")
+        boundary_gdf = get_boundary_gdf_from_database(area, level)
+        if boundary_gdf is not None:
+            logger.info(f"成功获取动态边界用于栅格转换: {area} ({level}) - 包含 {len(boundary_gdf)} 个几何体")
+            logger.info(f"边界数据范围: {boundary_gdf.bounds.to_dict('records')[0] if len(boundary_gdf) > 0 else 'N/A'}")
+        else:
+            logger.warning(f"未能获取动态边界,栅格转换将使用完整范围: {area} ({level})")
+            logger.warning("可能的原因: 1) 数据库中没有该地区数据 2) 地区名称拼写错误 3) level参数不正确")
+    else:
+        logger.info("未指定area和level参数,尝试获取默认边界(乐昌市)")
+        boundary_gdf = get_boundary_gdf_from_database("乐昌市", "county")
+        if boundary_gdf is not None:
+            logger.info(f"成功获取默认边界(乐昌市)用于栅格转换 - 包含 {len(boundary_gdf)} 个几何体")
+        else:
+            logger.info("未能获取默认边界,栅格转换将使用数据点范围")
+            
+    # 3. 使用csv_to_raster_workflow转换为GeoTIFF(支持动态边界)
     raster_dir = os.path.join(base_dir, "..", "static", "water", "Raster")
     template_tif = os.path.join(raster_dir, "meanTemp.tif")
     output_dir = os.path.dirname(cleaned_csv)  # 使用CSV所在目录作为输出目录
 
-    # 调用csv_to_raster_workflow
+    # 调用csv_to_raster_workflow,支持动态边界和插值控制
     workflow_result = csv_to_raster_workflow(
         csv_file=cleaned_csv,
         template_tif=template_tif,
         output_dir=output_dir,
-        boundary_shp=None,
-        resolution_factor=4.0,
-        interpolation_method='linear',
+        boundary_shp=None,  # 不使用Shapefile
+        boundary_gdf=boundary_gdf,  # 直接使用GeoDataFrame
+        resolution_factor=resolution_factor,  # 可配置的分辨率因子
+        interpolation_method=interpolation_method,  # 可配置的插值方法
         field_name='Prediction',
         lon_col=0,  # CSV中经度列索引
         lat_col=1,  # CSV中纬度列索引
         value_col=2,  # CSV中数值列索引
-        enable_interpolation=True
+        enable_interpolation=enable_interpolation  # 可配置的插值开关
     )
 
     # 获取输出的栅格文件路径
     output_tif = workflow_result['raster']
     logger.info(f"生成栅格文件: {output_tif}")
 
-    # 3. 生成栅格地图
+    # 4. 生成栅格地图(直接使用统一的MappingUtils接口,支持动态边界)
     map_output = os.path.join(raster_dir, f"{land_type}_Cd含量地图.jpg")
-    county_shp = os.path.join(raster_dir, "Lechang.shp")  # 县界SHP
-
-    mapping_raster(
-        shp_path=county_shp,
-        tif_path=output_tif,
-        color_map_name=color_map_name,
-        title_name=f"Irrigation water input flux",
-        output_path=map_output,
-        output_size=output_size
-
-    )
+    
+    try:
+        # 创建MappingUtils实例并直接调用
+        mapper = MappingUtils()
+        
+        # 转换颜色方案名称
+        colormap_key = COLORMAP_MAPPING.get(color_map_name, "green_yellow_red_purple")
+        
+        # 统一使用boundary_gdf,不再依赖Shapefile
+        map_result = mapper.create_raster_map(
+            shp_path=None,  # 不使用Shapefile
+            tif_path=output_tif,
+            output_path=os.path.splitext(map_output)[0],  # 去掉扩展名
+            colormap=colormap_key,
+            title="Irrigation water input flux",
+            output_size=output_size,
+            dpi=300,
+            enable_interpolation=enable_interpolation,  # 可配置的插值开关
+            interpolation_method=interpolation_method,  # 可配置的插值方法
+            boundary_gdf=boundary_gdf  # 使用从数据库获取的边界数据,如果为None则使用数据点范围
+        )
+        
+        # 如果生成的路径与期望不同,则重命名
+        if map_result != map_output and os.path.exists(map_result):
+            shutil.move(map_result, map_output)
+            logger.info(f"栅格地图已生成: {map_output}")
+            
+    except Exception as e:
+        logger.error(f"栅格地图生成失败: {str(e)}")
 
-    # 4. 生成直方图
+    # 5. 生成直方图(直接使用统一的MappingUtils接口)
     hist_output = os.path.join(raster_dir, f"{land_type}_Cd含量直方图.jpg")
-    plot_tif_histogram(
-        file_path=output_tif,
-        output_path=hist_output,
-        title=f"Irrigation water input flux"
-    )
+    
+    try:
+        # 直接调用统一绘图接口生成直方图
+        hist_result = mapper.create_histogram(
+            file_path=output_tif,
+            save_path=hist_output,
+            figsize=(8, 8),
+            xlabel='Cd(ug/L)',
+            ylabel='frequency density',
+            title='Irrigation water input flux',
+            bins=100,
+            dpi=300
+        )
+        
+        if hist_result:
+            logger.info(f"直方图已生成: {hist_output}")
+            
+    except Exception as e:
+        logger.error(f"直方图生成失败: {str(e)}")
 
     return cleaned_csv, workflow_result['shapefile'], output_tif, map_output, hist_output, used_coeff