.net解析雷达拼图3.0组合反射率产品,并按经纬度裁剪绘制组合反射率图

.net解析雷达拼图产品并绘制反射率图

1.组合反射率产品文件名称:

Z_RADA_C_BABJ_2025XXXXXXXXXX_P_DOR_ACHN_CREF_2025XXXX_XXXXXX.bin2

2.文件格式

拼图的雷达文件分为两个部分:文件头 + 数据块,文件头包括一下数据

文件固定标识 (4字节)
文件格式版本代码 (4字节)
文件字节数 (4字节)
拼图产品编号 (2字节)
坐标系类型 (2字节)
变量名 (8字节)
描述信息 (64字节)
块位置 (4字节)
块长度 (4字节)
时区 (4字节)
年份 (2字节)
月份 (2字节)
日期 (2字节)
 小时 (2字节)
 分钟 (2字节)
秒钟 (2字节)
观测秒数 (4字节)
观测日期 (2字节)
生成日期 (2字节)
生成秒数 (4字节)
最小纬度 (4字节,除以1000)12.2
最小经度 (4字节,除以1000)73
最大纬度 (4字节,除以1000)54.2
最大经度 (4字节,除以1000)135
中心点X (4字节,除以1000)104
中心点Y (4字节,除以1000)33.2
列数 (4字节)6200
行数 (4字节)4200
X方向分辨率 (4字节,除以10000)0.01
Y方向分辨率 (4字节,除以10000)0.01
高度 (2字节)
压缩标志 (2字节)1
雷达数量 (4字节)239
解压后字节数 (4字节)
数据缩放因子 (2字节)
未使用字段 (2字节)
区域ID (8字节)
单位 (8字节)
保留字段 (8字节)

数据块中在后续解压和解析中使用的数据有时间(年月日时分秒),经、纬度方向网格数,最大最小经纬度,数据缩放因子,创建对象将这些数据存储起来

3.完整代码

using System;
using System.Collections.Generic;
using System.Drawing;
using System.Drawing.Imaging;
using System.IO;
using System.Linq;
using System.Runtime.Intrinsics.X86;
using System.Text;
using ICSharpCode.SharpZipLib.BZip2;
using Newtonsoft.Json;
using Newtonsoft.Json.Linq;
public class RadarDataParser
{
    // 定义反射率等级 (dBZ)
    private static readonly double[] Levels = { 0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70 };

    // 定义与每个反射率等级对应的颜色 (RGBA)
    private static readonly Color[] Colors =
    {
        Color.FromArgb(0, 0, 239),     // 0-5 dBZ
        Color.FromArgb(65, 157, 241),  // 5-10 dBZ
        Color.FromArgb(100, 231, 235), // 10-15 dBZ
        Color.FromArgb(109, 250, 61),  // 15-20 dBZ
        Color.FromArgb(0, 216, 0),     // 20-25 dBZ
        Color.FromArgb(1, 144, 0),     // 25-30 dBZ
        Color.FromArgb(255, 255, 0),   // 30-35 dBZ
        Color.FromArgb(231, 192, 0),   // 35-40 dBZ
        Color.FromArgb(255, 144, 0),   // 40-45 dBZ
        Color.FromArgb(255, 0, 0),     // 45-50 dBZ
        Color.FromArgb(214, 0, 0),     // 50-55 dBZ
        Color.FromArgb(192, 0, 0),     // 55-60 dBZ
        Color.FromArgb(255, 0, 240),   // 60-65 dBZ
        Color.FromArgb(150, 0, 180),   // 65-70 dBZ
        Color.FromArgb(173, 144, 240)  // >70 dBZ
    };

    /// <summary>
    /// 解析雷达数据文件并生成雷达图像
    /// </summary>
    /// <param name="filePath">雷达数据文件路径</param>
    public static void Parse(string filePath, double lonMin, double lonMax, double latMin, double latMax, string savePath)
    {
        byte[] fileBytes = File.ReadAllBytes(filePath);

        using (var stream = new MemoryStream(fileBytes))
        using (var reader = new BinaryReader(stream))
        {
            var header = ReadHeader(reader);
            stream.Seek(256, SeekOrigin.Begin);

            byte[] compressedData = reader.ReadBytes(fileBytes.Length - 256);

            byte[] decompressedData = DecompressBZip2(compressedData);

            if (decompressedData == null)
            {
                Console.WriteLine("数据解压失败");
                return;
            }

            short[] data = new short[decompressedData.Length / 2];
            Buffer.BlockCopy(decompressedData, 0, data, 0, decompressedData.Length);

            double[,] array2D = new double[header.nY, header.nX];

            for (int i = 0; i < header.nY; i++)
            {
                for (int j = 0; j < header.nX; j++)
                {
                    short value = data[i * header.nX + j];

                    double scaledValue = value / (double)header.scale;

                    array2D[i, j] = (scaledValue >= 0 && scaledValue <= 75) ? scaledValue : double.NaN;
                }
            }

            // 生成经纬度网格
            double[] lon = new double[header.nX];
            double[] lat = new double[header.nY];

            // 经度网格 (从左到右)
            for (int i = 0; i < header.nX; i++)
            {
                lon[i] = header.lonMin + i * (header.lonMax - header.lonMin) / (header.nX - 1);
            }

            // 纬度网格 (从下到上)
            for (int i = 0; i < header.nY; i++)
            {
                lat[i] = header.latMin + i * (header.latMax - header.latMin) / (header.nY - 1);
            }

            // 裁剪裁剪地区数据
            var croppedData = CropDataForHenan(array2D, lon, lat, lonMin, lonMax, latMin, latMax);
            double[,] croppedArray = croppedData.Data;
            int cropWidth = croppedData.Width;
            int cropHeight = croppedData.Height;
            int startX = croppedData.StartX;
            int startY = croppedData.StartY;

            DateTime dt = new DateTime(header.yr, header.mon, header.day, header.hr, header.min, 0).AddHours(8);//雷达日期
            String ttpp = $"{Path.DirectorySeparatorChar}{dt.ToString("yyyy")}{Path.DirectorySeparatorChar}{dt.ToString("MM")}{Path.DirectorySeparatorChar}";
            //  创建位图并绘制裁剪雷达数据
            using (Bitmap bitmap = new Bitmap(cropWidth, cropHeight, PixelFormat.Format32bppArgb))
            {
                // 遍历每个像素点
                for (int y = 0; y < cropHeight; y++)
                {
                    for (int x = 0; x < cropWidth; x++)
                    {
                        // 获取数据值
                        double value = croppedArray[y, x];

                        // 根据值获取对应颜色
                        Color color = GetColorForValue(value);

                        // 设置像素颜色
                        bitmap.SetPixel(x, y, color);
                    }
                }

                //保存为透明背景的PNG图像
                string descPath = string.Format(@"{0}\{1}", savePath, ttpp);
                if (!Directory.Exists(descPath))
                {
                    Directory.CreateDirectory(descPath);
                }
                string outputPath = Path.Combine(descPath, dt.ToString("yyyyMMddHHmm") + ".png");
                
                if (File.Exists(outputPath)) { File.Delete(outputPath); }
                bitmap.Save(outputPath, ImageFormat.Png);
                Console.WriteLine($"裁剪雷达图像已成功保存为: {outputPath}");
               
            }
            try
            {
                File.Delete(filePath);
            }
            catch (Exception ex)
            {

            }
        }

    }
    /// <summary>
    /// 裁剪裁剪地区数据
    /// </summary>
    private static CroppedData CropDataForHenan(double[,] fullData, double[] lon, double[] lat, double lonMin, double lonMax, double latMin, double latMax)
    {
        int xStart = -1, xEnd = -1;
        for (int i = 0; i < lon.Length; i++)
        {
            if (lon[i] >= lonMin && xStart == -1) xStart = i;
            if (lon[i] <= lonMax) xEnd = i;
            else break;
        }

        int yStart = -1, yEnd = -1;
        for (int i = 0; i < lat.Length; i++)
        {
            if (lat[i] >= latMin && yStart == -1) yStart = i;
            if (lat[i] <= latMax) yEnd = i;
            else break;
        }

        if (xStart == -1 || xEnd == -1 || yStart == -1 || yEnd == -1)
        {
            return new CroppedData
            {
                Data = fullData,
                Width = fullData.GetLength(1),
                Height = fullData.GetLength(0),
                StartX = 0,
                StartY = 0
            };
        }

        int width = xEnd - xStart + 1;
        int height = yEnd - yStart + 1;

        double[,] cropped = new double[height, width];

        for (int y = yStart; y <= yEnd; y++)
        {
            for (int x = xStart; x <= xEnd; x++)
            {
                cropped[y - yStart, x - xStart] = fullData[y, x];
            }
        }

        return new CroppedData
        {
            Data = cropped,
            Width = width,
            Height = height,
            StartX = xStart,
            StartY = yStart
        };
    }
    // <summary>
    /// 裁剪数据结果容器
    /// </summary>
    private class CroppedData
    {
        public double[,] Data { get; set; }
        public int Width { get; set; }
        public int Height { get; set; }
        public int StartX { get; set; }
        public int StartY { get; set; }
    }
    /// <summary>
    /// 根据反射率值获取对应颜色
    /// </summary>
    /// <param name="value">反射率值 (dBZ)</param>
    /// <returns>对应的颜色,无效值为透明</returns>
    private static Color GetColorForValue(double value)
    {
        if (double.IsNaN(value))
            return Color.Transparent;

        for (int i = 0; i < Levels.Length - 1; i++)
        {
            if (value >= Levels[i] && value < Levels[i + 1])
            {
                return Colors[i];
            }
        }

        if (value >= Levels[Levels.Length - 1])
        {
            return Colors[Colors.Length - 1];
        }

        return Color.Transparent;
    }

    /// <summary>
    /// 解压BZip2格式数据
    /// </summary>
    /// <param name="data">压缩的字节数组</param>
    /// <returns>解压后的字节数组,失败返回null</returns>
    private static byte[] DecompressBZip2(byte[] data)
    {
        try
        {
            using (var input = new MemoryStream(data))
            using (var output = new MemoryStream())
            {
                BZip2.Decompress(input, output, false);
                return output.ToArray();
            }
        }
        catch (Exception ex)
        {
            Console.WriteLine($"BZip2解压错误: {ex.Message}");
            return null;
        }
    }

    /// <summary>
    /// 读取雷达数据文件头
    /// </summary>
    /// <param name="reader">二进制读取器</param>
    /// <returns>包含头信息的RadarHeader对象</returns>
    private static RadarHeader ReadHeader(BinaryReader reader)
    {
        var header = new RadarHeader();

        reader.BaseStream.Seek(0, SeekOrigin.Begin);

        string label = ReadString(reader, 4, "gb2312");
        string version = ReadString(reader, 4, "gb2312");
        int file_bytes = reader.ReadInt32();
        short mosaic_id = reader.ReadInt16();
        short coordinate = reader.ReadInt16();
        string varname = ReadString(reader, 8, "gb2312");
        string description = ReadString(reader, 64, "gb2312");
        int BlockPos = reader.ReadInt32();
        int BlockLen = reader.ReadInt32();
        int TimeZone = reader.ReadInt32();
        header.yr = reader.ReadInt16();
        header.mon = reader.ReadInt16();
        header.day = reader.ReadInt16();
        header.hr = reader.ReadInt16();
        header.min = reader.ReadInt16();
        short sec = reader.ReadInt16();
        int ObsSeconds = reader.ReadInt32();
        short ObsDates = reader.ReadInt16();
        short GenDates = reader.ReadInt16();
        int GenSeconds = reader.ReadInt32();
        header.latMin = reader.ReadInt32() / 1000.0;
        header.lonMin = reader.ReadInt32() / 1000.0; 
        header.latMax = reader.ReadInt32() / 1000.0; 
        header.lonMax = reader.ReadInt32() / 1000.0;
        double cx = reader.ReadInt32() / 1000.0;
        double cy = reader.ReadInt32() / 1000.0;
        header.nX = reader.ReadInt32();
        header.nY = reader.ReadInt32();
        double dx = reader.ReadInt32() / 10000.0;
        double dy = reader.ReadInt32() / 10000.0; 
        short height = reader.ReadInt16();
        short compress = reader.ReadInt16();
        int num_of_radars = reader.ReadInt32();
        int UnZipBytes = reader.ReadInt32();
        header.scale = reader.ReadInt16();
        short unUsed = reader.ReadInt16(); 
        string RgnID = ReadString(reader, 8, "gb2312");
        string units = ReadString(reader, 8, "gb2312");
        string reserved = ReadString(reader, 8, "gb2312");

        // 验证读取位置是否正确 (总头长度应为256字节,有部分保留字节)
        if (reader.BaseStream.Position != 256)
        {
            // 强制跳转到256字节位置
            reader.BaseStream.Seek(256, SeekOrigin.Begin);
        }

        return header;
    }
    /// <summary>
    /// 读取固定长度的GB2312编码字符串 
    /// </summary>
    /// <param name="reader">二进制读取器</param>
    /// <param name="length">字符串长度</param>
    /// <param name="encoding">编码格式</param>
    /// <returns>处理后的字符串</returns>
    private static string ReadString(BinaryReader reader, int length, string encoding = "gb2312")
    {
        byte[] bytes = reader.ReadBytes(length);
        Encoding.RegisterProvider(CodePagesEncodingProvider.Instance);
        Encoding enc = Encoding.GetEncoding(encoding);
        string str = enc.GetString(bytes);

        return str.Replace("\0", "").Replace(" ", "");
    }
    /// <summary>
    /// 雷达文件头信息结构 (包含所有必要字段)
    /// </summary>
    private class RadarHeader
    {
        public short yr;         // 年份
        public short mon;        // 月份
        public short day;        // 日期
        public short hr;         //小时
        public short min;       //分钟
        public int nX;           // 经度方向网格数
        public int nY;           // 纬度方向网格数
        public double latMin;   // 最小纬度
        public double latMax;   // 最大纬度
        public double lonMin;   // 最小经度
        public double lonMax;   // 最大经度
        public short scale;     // 数据缩放因子

        public override string ToString()
        {
            return $"日期: {yr}-{mon}-{day}, 网格: {nX}x{nY}, 纬度范围: [{latMin},{latMax}], 经度范围: [{lonMin},{lonMax}], 缩放因子: {scale}";
        }
    }
}

class Program
{
    static JObject valueconfig;
    static void Main(string[] args)
    {
        try
        {
            System.AppContext.SetSwitch("System.Drawing.EnableUnixSupport", true);
            string pathh = AppDomain.CurrentDomain.SetupInformation.ApplicationBase;

            using (System.IO.StreamReader file = new System.IO.StreamReader(AppDomain.CurrentDomain.BaseDirectory + "config.json", Encoding.UTF8))
            {
                valueconfig = JsonConvert.DeserializeObject<JObject>(file.ReadToEnd());
                file.Close();
            }
            Console.WriteLine("雷达数据解析");

            string[] files = Directory.GetFiles(valueconfig["pathfile"].ToString(), "*.*");
            foreach (string file in files)
            {
                if (!File.Exists(file))
                {
                    Console.WriteLine($"文件不存在: {file}");
                    continue;
                }
                double henanLonMin = 109.46, henanLonMax = 124.22, henanLatMin = 27.05, henanLatMax = 38.80;
                try
                {
                    henanLonMin = double.Parse(valueconfig["minlon"].ToString());  // 裁剪最小经度
                    henanLonMax = double.Parse(valueconfig["maxlon"].ToString()); // 裁剪最大经度
                    henanLatMin = double.Parse(valueconfig["minlat"].ToString());  // 裁剪最小纬度
                    henanLatMax = double.Parse(valueconfig["maxlat"].ToString());   // 裁剪最大纬度
                }
                catch (Exception ex) { }

                Console.WriteLine($"处理裁剪雷达数据 (范围: 经度[{henanLonMin}-{henanLonMax}], 纬度[{henanLatMin}-{henanLatMax}])");

                // 处理雷达数据并生成裁剪地区图像
                RadarDataParser.Parse(file, henanLonMin, henanLonMax, henanLatMin, henanLatMax, valueconfig["savePath"].ToString());
            }
        }catch(Exception ex) { Console.WriteLine(ex.Message); }
        

    }
}

配置文件格式

{
  "source": "天擎雷达拼图解析",
  "pathfile": "E:\\雷达文件",
  "savePath": "E:\\雷达解析",
  "minlon": 109.46, //最小73
  "maxlon": 124.22, //最大135
  "minlat": 27.05, //最小12.2
  "maxlat": 38.80 //最大54.2
}

代码中用到了包Newtonsoft.Json处理json文件(主要是配置文件),SharpZipLib解压文件,Syatem.Drawing.Common用来画图

4.整体思路:

 读取文件头信息

读取压缩数据部分

解压BZip2数据

创建二维数组存储处理后的数据

生成经纬度网格

裁剪裁剪地区数据(可以不用裁剪,直接用全国数据)

根据色标卡生成对应的色板图

删除文件(每6分钟一组数据,后期会比较占存储空间)

5.最后生成的雷达图

要对天气雷达拼图V3.0的CREF参数进行绘,首先需要导入必要的库设置一些基本参数。这里我们假设你已经有一个名为`Cref_data`的数据结构(可能是数组或pandas DataFrame),其中包含了CREF参数。以下是一段示例Python代码: ```python import numpy as np import cartopy.crs as ccrs import matplotlib.pyplot as plt from cinrad.visualize import PPI # 假设你已经有了CREF数据 # Cref_data = ... (实际数据填充) # 设置圆柱投影 projection = ccrs.PlateCarree() # 创建一个新的地 fig = plt.figure(figsize=(10, 10)) ax = fig.add_subplot(1, 1, 1, projection=projection) # 绘制CREF数据 clevs = np.arange(0, 61, 5) # 定义等值线范围 cs = ax.contourf(Cref_data[&#39;lon&#39;], Cref_data[&#39;lat&#39;], Cref_data[&#39;CREF&#39;], transform=projection, cmap=&#39;viridis&#39;, levels=clevs) # 添加标题和坐标轴标签 ax.set_title(&#39;CREF Reflectivity in V3.0 Radar Network&#39;) ax.coastlines(resolution=&#39;10m&#39;) # 显示海岸线 ax.gridlines(draw_labels=True) # 添加网格线 # 添加颜色bar plt.colorbar(cs, orientation="vertical", shrink=0.8, pad=0.05, aspect=40, extend="both") # 圆形坐标转换 circle = cartopy.feature.GeoFeature(cfeature.ARC_GEO, facecolor=&#39;none&#39;, edgecolor=&#39;black&#39;, linewidth=0.5) ax.add_feature(circle, extent=(-180, 180, -90, 90)) plt.show() ``` 这段代码首先设置了圆柱投影,然后创建了一个地添加了一个子。接着,它使用`contourf()`函数以CREF数据作为输入,在地绘制等值线。最后,添加了标题、海岸线和网格线,显示了一个圆形的地球轮廓。 请注意,你需要根据实际数据结构替换`Cref_data`变量,可能调整`clevs`、`colormap`等参数以适应你的需求。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值