雷达系列 | 拆解NMC出品的雷达库metradar 架构

┌─────────────────────────────────────────┐
│       Application Layer (Project)       │  业务流程封装
│   (拼图制作 / 临近预报 / 风场反演)        │
├─────────────────────────────────────────┤
│     Presentation Layer (Graph)           │  可视化接口
│   (雷达绘图 / 综合分析图 / 等值线图)       │
├─────────────────────────────────────────┤
│      Business Logic Layer (Core)        │  核心算法
│   (质量控制 / 拼图合并 / 散度涡度计算)     │
├─────────────────────────────────────────┤
│       Data Access Layer (IO)             │  数据I/O
│   (FMT解码 / PUP读取 / SWAN拼图)          │
├─────────────────────────────────────────┤
│      Utility Layer (Util)               │  通用工具
│   (地理变换 / 异常定义 / 色标解析)         │
├─────────────────────────────────────────┤
│     Configuration Layer (Config)        │  配置管理
│   (资源路径 / 全局参数)                   │
└─────────────────────────────────────────┘

class CNRADLevel2Reader:
    """中国雷达基数据读取器 - 适配器模式"""
    
    def read(self, file_path):
        # 1. 读取原始二进制数据
        raw_data = self._read_binary(file_path)
        
        # 2. 解析 FMT 格式（中国特有的二进制格式）
        parsed = self._parse_fmt(raw_data)
        
        # 3. 转换为 PyArt 标准对象
        return pyart.Radar(
            time=parsed['time'],
            fields=parsed['fields'],
            metadata=parsed['metadata'],
            scan_type=parsed['scan_type']
        )

def decode_fmt(file_path, version=None):
    """FMT 解码工厂函数"""
    if version is None:
        version = detect_fmt_version(file_path)
    
    if version == 'v2.0':
        decoder = FMTv2Decoder()
    elif version == 'v2.1':
        decoder = FMTv2_1Decoder()
    else:
        raise UnsupportedFMTVersion(version)
    
    return decoder.decode(file_path)

def calculate_divergence(u, v, dx, dy):
    """
    计算散度
    
    Args:
        u: 纬向风场分量
        v: 经向风场分量
        dx, dy: 网格间距
    
    Returns:
        divergence: 散度场
    """
    du_dx = np.gradient(u, dx, axis=1)
    dv_dy = np.gradient(v, dy, axis=0)
    return du_dx + dv_dy

def calculate_vorticity(u, v, dx, dy):
    """
    计算涡度
    """
    dv_dx = np.gradient(v, dx, axis=1)
    du_dy = np.gradient(u, dy, axis=0)
    return dv_dx - du_dy

class MosaicMergeStrategy:
    """拼图合并策略接口"""
    def merge(self, radars: List[Radar]) -> Radar:
        raise NotImplementedError

class MaxValueStrategy(MosaicMergeStrategy):
    """最大值策略：选择各点的最大反射率"""
    def merge(self, radars):
        # 选择最大反射率
        pass

class NearestStrategy(MosaicMergeStrategy):
    """最近邻策略：选择距离最近的雷达值"""
    def merge(self, radars):
        # 距离最近的雷达值
        pass

class WeightedAverageStrategy(MosaicMergeStrategy):
    """加权平均策略：按距离加权平均"""
    def merge(self, radars):
        # 按距离加权
        pass

# 使用示例
strategy = MaxValueStrategy()
mosaic = strategy.merge(radar_list)

class RadarMapBuilder:
    """雷达地图构建器"""
    
    def __init__(self):
        self.radar = None
        self.overlays = []
        self.colormap = None
        
    def set_radar(self, radar):
        """设置雷达数据"""
        self.radar = radar
        return self
    
    def add_aws_overlay(self, aws_data):
        """添加自动站数据叠加"""
        self.overlays.append(aws_data)
        return self
    
    def set_colormap(self, cmap):
        """设置色标"""
        self.colormap = cmap
        return self
    
    def build(self):
        """构建完整图形"""
        fig = self._create_base_map()
        self._add_radar_data(fig)
        self._add_overlays(fig)
        self._add_colorbar(fig)
        return fig

# 使用示例：链式调用，简洁优雅
fig = (RadarMapBuilder()
       .set_radar(radar)
       .add_aws_overlay(aws)
       .set_colormap('NWSRef')
       .build())

class BaseRadarPlotter:
    """雷达绘图基类"""
    
    def plot(self, radar, **kwargs):
        """模板方法：定义绘图流程骨架"""
        self._prepare_data(radar)
        self._setup_figure()
        self._draw_map_background()
        self._draw_radar_data()
        self._add_overlays()
        self._add_colorbar()
        return self.fig
        
    def _prepare_data(self, radar):
        """子类可覆盖的数据预处理"""
        pass

class RadarAWSPlotter(BaseRadarPlotter):
    """雷达+自动站绘图器"""
    
    def _add_overlays(self):
        """添加自动站观测数据"""
        self._draw_aws_symbols()

class MosaicPipeline:
    """拼图制作管道"""
    
    def __init__(self, config):
        self.config = config
        self.stages = []
        
    def add_stage(self, stage):
        """添加处理阶段"""
        self.stages.append(stage)
        return self
        
    def run(self):
        """执行完整管道"""
        data = None
        for stage in self.stages:
            data = stage.process(data)
        return data

# 使用示例：构建完整的拼图制作流程
pipeline = (MosaicPipeline(config)
            .add_stage(DataDownloadStage())     # 下载数据
            .add_stage(QCStage())                # 质量控制
            .add_stage(MergeStage())             # 拼图合并
            .add_stage(PlotStage())              # 生成图像
            .add_stage(OutputStage()))           # 输出产品

pipeline.run()

┌──────────────┐
│ 原始数据源    │
│ (FMT/PUP/SWAN)│
└──────┬───────┘
       │ IO 层适配器
       │ 解析并统一
       ▼
┌──────────────┐
│ PyArt 统一   │
│ Radar 对象    │
└──────┬───────┘
       │ Core 层算法
       │ QC/反演/分析
       ▼
┌──────────────┐
│ 处理结果     │
│ (质量控制的   │
│  雷达数据)   │
└──────┬───────┘
       │ Graph 层可视化
       │ 绘图/叠加
       ▼
┌──────────────┐
│ 图形产品     │
│ (PNG/动画)  │
└──────┬───────┘
       │ Project 层工作流
       │ 批处理/输出
       ▼
┌──────────────┐
│ 业务输出     │
└─────────────┘

┌─────────────────────────────────────┐
│        metradar 应用层              │
│  (Project + 业务流程封装)            │
├─────────────────────────────────────┤
│     metradar 核心层                  │
│  (Core + Graph + IO 本地化功能)      │
├─────────────────────────────────────┤
│       PyArt 统一抽象层              │
│  (Radar 对象统一接口)                │
├─────────────────────────────────────┤
│     第三方工具包层                   │
│  pyart | wradlib | pydda | pysteps  │
└─────────────────────────────────────┘

import wradlib.dp as dp
cluttermap = dp.process_clutter(data, method='histogram', 
                                 scale=1.0, max_range=200, 
                                 beamwidth=1.0, ... )

from metradar.core import remove_clutter
cluttermap = remove_clutter(data, method='histogram')

def chinese_radar_qc(radar):
    """中国雷达专用质量控制流程"""
    
    # 1. 调用 PyArt 标准质量控制
    radar = pyart.correct.correct_radar(radar)
    
    # 2. 中国特有的处理
    radar = _remove_sea_clutter(radar)      # 海洋杂波去除
    radar = _correct_beam_blockage(radar)   # 地形波束遮挡订正
    radar = _detect_precip_type(radar)      # 降水类型识别
    radar = _apply_china_qc_rules(radar)   # 应用中国 QC 规则
    
    return radar

def hybrid_wind_retrieval(radar, aws_data):
    """混合风场反演：结合多种方法的优点"""
    
    # 1. wradlib 单雷达风场反演
    wind1 = wradlib.dp.vel2z(radar)
    
    # 2. pydda 双多普勒分析（如果有多个雷达）
    if len(radar_list) >= 2:
        wind2 = pydda.Grid.get_dd_wind(radar_list)
        wind1 = _merge_winds(wind1, wind2)
    
    # 3. 自动站数据融合
    wind_final = _assimilate_aws(wind1, aws_data)
    
    return wind_final

metradar/
├── config.py              # 配置层
├── core/                  # 核心算法层
│   ├── get_cross_section.py    # 剖面图算法
│   ├── mosaic_merge.py         # 拼图合并算法
│   └── oa_dig_func.py          # 客观分析函数
├── io/                    # 数据 I/O 层
│   ├── cnrad_level2.py         # 中国雷达基数据
│   ├── decode_fmt_pyart.py     # FMT 格式解码
│   ├── decode_pup_rose.py      # PUP/ROSE 解码
│   ├── read_swan.py            # SWAN 拼图读取
│   └── rose_structer.py        # ROSE 数据结构
├── graph/                 # 可视化层
│   ├── draw_comp_mosaic.py     # 综合拼图绘制
│   ├── draw_latlon_func.py     # 经纬度投影
│   ├── draw_radar_aws.py       # 雷达+自动站
│   └── draw_radar_comp_func.py # 综合绘图函数
├── util/                  # 工具层
│   ├── comm_func.py            # 通用函数
│   ├── exceptions.py           # 自定义异常
│   ├── geo_transforms_pyart.py # 地理变换
│   ├── parse_pal.py            # 色标解析
│   ├── radar_common.py         # 雷达通用工具
│   └── trans_new_mosaic_nc.py  # 拼图格式转换
└── project/               # 应用层
    ├── make_mosaic/       # 拼图制作项目
    ├── make_vpr_aws/      # VPR+自动站项目
    ├── nowcasting/        # 临近预报项目
    ├── qpe/               # 降水估测项目
    └── wind_retrieval/    # 风场反演项目