雷达组网拼图3.0数据掌握和python解析处理

废话不多说，先展示雷达图

以反射率为例：

核对数据格式

Z_RADA_C_BABJ_20240705043615_P_DOR_ACHN_CREF_20240705_043000.bin

数据分析认识

 1. 组网产品分类：
	组网产品包括组网混合扫描反射率（HSR），组网组合反射率（CR）、组网最大反射率高度（CRH）、组网回波顶高产品（ET）、组网垂直累积液态水含量（VIL）、组网一小时降水（OHP）、组网等高面反射率（CAP）。数据频次为10分钟
 2. 文件名结构：
	Z_RADA_C_BABJ_{YYYYMMDDhhmmss}_P_DOR_{AREACODE}_{DTYPE}_{YYYYMMDD_hhmmss}.bin，{}内表示为命名规则中的变量，定义依次如下：
	YYYYMMDDhhmmss：产品生成时间，年月日时分秒
	AREACODE：区域号，如全国区域：ACHN
	DTYPE：雷达产品类型简写，如QREF
	YYYYMMDD_hhmmss：资料观测时间，年月日_时分秒
 3. 文件格式：
	组网拼图产品文件格式定义为：文件头 + 数据块。文件头中定义了文件卷标（文件固定标识、文件格式版本代码、文件字节数）、
	产品描述（拼图产品编号、坐标类型、产品代码、产品描述、产品数据起始位置、产品数据字节数）、数据时间（数据时钟、观测时间年、月、日等）、
	数据区信息（数据区的南、西、北、东边界等），具体定义如下表所示。文件头共占用256字节，实际定义变量占176字节，
	后面的80字节为预留空间。数据块为nX*nY（文件头中定义变量）个short类型数组，根据文件头中的Compress标志位判断是否压缩

总体数据包括文件头和数据块，数据为bin格式二进制，根据文件格式将二进制转化为实际的数据类型，其余不过多分析，到这里，我们基本了解了组网产品的基本情况，接下来我们直接读取和出雷达图。

#  雷达拼图3.0数据绘制雷达图*************

# 设置日志
logger = logging.getLogger('RADA_L3_MST_CREF_QC')
logging.basicConfig(level=logging.INFO)


# 假设雷达颜色映射表 radarColoMap 已定义，并且是Python列表格式

def group(array, sub_group_length):
    for i in range(0, len(array), sub_group_length):
        yield array[i:i + sub_group_length]


def read_string(buffer, length, encoding='gb2312'):
    return buffer[:length].decode(encoding).replace('\x00', '').replace(' ', '')


def unbzip2_decompress(data):
    try:
        return bz2.decompress(data)
    except Exception as e:
        logger.error(f"Error decompressing bz2 data: {e}")
        return None


def get_path(path, path1):
    # 找到path1在path中的结束位置
    end_index = path.find(path1) + len(path1)
    # 截取从path1之后的部分
    sub_path = path[end_index:]
    return sub_path


class BinaryReader:
    def __init__(self, buffer, is_little_endian=True):
        self.buffer = buffer
        self.index = 0
        self.is_little_endian = is_little_endian

    def seek(self, offset, origin=0):
        if origin == 0:  # Set absolute position
            self.index = offset
        elif origin == 1:  # Set relative to current position
            self.index += offset
        elif origin == 2:  # Set relative to end of file
            self.index = len(self.buffer) + offset

    def read_bytes(self, length):
        data = self.buffer[self.index:self.index + length]
        self.index += length
        return data

    def read_string(self, length):
        return read_string(self.read_bytes(length), length)

    def read_int(self):
        fmt = '<i' if self.is_little_endian else '>i'
        value, = struct.unpack_from(fmt, self.buffer, self.index)
        self.index += struct.calcsize(fmt)
        return value

    def read_int16(self):
        fmt = '<h' if self.is_little_endian else '>h'
        value, = struct.unpack_from(fmt, self.buffer, self.index)
        self.index += struct.calcsize(fmt)
        return value

    def readHeader(br):
        # 读取文件头信息
        label = br.read_string(4)  # 文件固定标识
        version = br.read_string(4)  # 文件格式版本代码
        file_bytes = br.read_int()  # 文件字节数
        mosaic_id = br.read_int16()  # 拼图产品编号
        coordinate = br.read_int16()
        varname = br.read_string(8)
        description = br.read_string(64)
        BlockPos = br.read_int()
        BlockLen = br.read_int()
        TimeZone = br.read_int()
        yr = br.read_int16()
        mon = br.read_int16()
        day = br.read_int16()
        hr = br.read_int16()
        min = br.read_int16()
        sec = br.read_int16()
        ObsSeconds = br.read_int()
        ObsDates = br.read_int16()
        GenDates = br.read_int16()
        GenSeconds = br.read_int()
        latMin = br.read_int() / 1000
        lonMin = br.read_int() / 1000
        latMax = br.read_int() / 1000
        lonMax = br.read_int() / 1000
        cx = br.read_int() / 1000
        cy = br.read_int() / 1000
        # 读取数据区前的信息
        nX = br.read_int()  # 列数
        nY = br.read_int()  # 行数
        dx = br.read_int() / 10000
        dy = br.read_int() / 10000
        height = br.read_int16()
        compress = br.read_int16()
        num_of_radars = br.read_int()
        UnZipBytes = br.read_int()
        scale = br.read_int16()
        unUsed = br.read_int16()
        RgnID = br.read_string(8)
        units = br.read_string(8)
        reserved = br.read_string(8)
        return yr, mon, day, nX, nY, latMin, latMax, lonMin, lonMax, scale

    # 读取数据字节
    def readData(buffer):
        # 读取数据区
        # br.seek(256)  # 跳过未压缩的256字节头部
        compressed_data = unbzip2_decompress(buffer[256:])
        return compressed_data

    # 解压bz2数据
    def unbzip2_decompress(data):
        try:
            return bz2.decompress(data)
        except Exception as e:
            logger.error(f"Error decompressing bz2 data: {e}")
            return None


class RadarWarnRequest(BaseModel):
    filePath: Any
    radarValue: Any
    lonMin: Any
    lonMax: Any
    latMin: Any
    latMax: Any
    provinceId: Any


class RADA_L3_MST_CREF_QC:
    @staticmethod
    def parse(file_path, radarValue, lonMin, lonMax, latMin, latMax, provinceId):
        with open(file_path, 'rb') as f:
            buffer = f.read()
        br = BinaryReader(buffer, is_little_endian=True)
        yr, mon, day, nX, nY, lat_Min, lat_Max, lon_Min, lon_Max, scale = readHeader(br)
        compressed_data = readData(buffer)
        if compressed_data:
            data = np.frombuffer(compressed_data, dtype=np.int16)
            array_2d = data.reshape((nY, nX))
            array_2d = array_2d / scale  # 数据放大倍数
            array_2d = np.where((array_2d >= 0) & (array_2d <= 75), array_2d, np.nan)

            # array_2d = ma.masked_where((array_2d > 75) | (array_2d < -5), array_2d)  # 筛除无效值
            # array_2d = ma.filled(array_2d)
            # array_2d = np.round(array_2d, decimals=0)  # 保留两位
            # 使用flipud函数进行翻转
            array_2d = np.flipud(array_2d)
            lon = np.linspace(lon_Min, lon_Max, nX)
            lat = np.linspace(lat_Min, lat_Max, nY)
            # # 裁切完的数据
            cropped_array_2d, cropped_lon, cropped_lat = crop_masked_array_by_bounds(array_2d, lon, lat, lonMin, lonMax,
                                                                              latMin, latMax)
            # 绘制河南的
            proj = ccrs.PlateCarree()  # 创建投影，选择cartopy的platecarree投影
            fig = plt.figure(figsize=(10, 8))  # 创建页面，可以自己选择大小
            # ax = fig.subplots(subplot_kw={'projection': proj})  #子图
            ax = fig.add_subplot(1, 1, 1, projection=proj)
            levs = [0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 55.0, 60.0, 65.0, 70.0]
            cols = [(0, 0, 239 / 255), (65 / 255, 157 / 255, 241 / 255), (100 / 255, 231 / 255, 235 / 255),
                    (109 / 255, 250 / 255, 61 / 255), (0, 216 / 255, 0),
                    (1 / 255, 144 / 255, 0), (255 / 255, 255 / 255, 0), (231 / 255, 192 / 255, 0),
                    (255 / 255, 144 / 255, 0), (255 / 255, 0, 0), (214 / 255, 0, 0), (192 / 255, 0, 0),
                    (255 / 255, 0, 240 / 255),
                    (150 / 255, 0, 180 / 255), (173 / 255, 144 / 255, 240 / 255)]
            # 创建自定义颜色映射
            custom_cmap = ListedColormap(cols)
            norm = BoundaryNorm(levs, custom_cmap.N, clip=False)
            # 绘制等高线填充图
            plt.contourf(lon, lat, array_2d, cmap=custom_cmap, norm=norm, levels=levs)
            # 设置透明背景并保存图片
            plt.savefig('../v3/china.png', transparent=True, bbox_inches='tight',
                        pad_inches=0)
        else:
            logger.error("Failed to decompress data")
            return