Matlab读取天气雷达基数据实现反射率因子三维可视化

rabialee

%读取雷达基数据，并按坐标输出基本反射率，画反射率因子图和三维锥形回波图
%读数据,一个基数据共4032条径向，每条径向2432字节，共9805824字节。
[FileName,PathName] = uigetfile('.bin','选择你的文件');                    %自行选择数据
fid=fopen(strcat(PathName,FileName),'r');
A = fread(fid,'uchar');
B=zeros(2432,4032);                                                        %将雷达基数据按径向存放，共有4032条，每一行代表一条径向
for a=1:4032                                                               %每一行代表一条径向
    b=1:2432;    %每个径向有2432个字节
    B(b,a)=A(b+(a-1)*2432);
end  
%XLSWRITE('C:\Users\Rabia Lee\Desktop\radar.xls',B);
%C=B(:,1:100);
%B2=B(1:2432,1200:1300)
switch B(73,1)                                                             %读取体扫模式，73，74位为体扫数据
    case 11
        phi=[0.50,0.50,1.45,1.45,2.40,3.35,4.30,5.25,6.2,7.5,8.7,10,12,14,16.7,19.5]  
    case 21
        phi=[0.50,0.50,1.45,1.45,2.40,3.35,4.30,6.00,9.00,14.6,19.5]       %体扫21模式的具体仰角，共有11个仰角，单位为度
    case 31
        phi=[0.50,0.50,1.50,1.50,2.50,2.50,3.50,4.50]
    case 32
        phi=[0.50,0.50,2.50,3.50,4.50]
end

f1=1;                                                                      %计数单位
g1=zeros(4032,460);                                                        %提取4032条径向的仰角
h1=zeros(4032,460);                                                        %提取4032条径向的方位角
i1=zeros(4032,460);                                                        %4032条径向每个距离点的距离
j1=zeros(4032,460);                                                        %提取4032条径向的反射率
%每一条径向460个点的仰角、方位角、距离的确定
for a1=1:4032
    b1=B(45,a1)+256*B(46,a1);                                              %每一条径向的仰角序数，1到11
    c1=(B(37,a1)+256*B(38,a1))/8*180/4096*2*pi/360;                        %每一条径向的方位角，弧度单位
    d1=B(55,a1)+256*B(56,a1);                                              %每一条径向的距离库数，1到460
    if d1==0
       continue
    end
    for e1=1:460                                                           %提取一条径向的每一个距离点的位置坐标
        g1(a1,e1)=phi(b1)*2*pi/360;                                        %提取仰角，由仰角序数计算仰角，弧度单位
        h1(a1,e1)=c1;                                                      %提取方位角，弧度为单位
        i1(a1,e1)=0.5+e1-1;                                                %计算距离，每一个距离点相隔0.5公里，共有460个距离点
    %提取反射率数据
    if e1>d1                                                               %距离库数以外的反射率为0
    j1(a1,e1)=0;                                                           %有4032条径向，每条径向有460个距离点
    else
        if B(128+e1,a1)==0                                                 %无数据
           j1(a1,e1)=0;
        elseif B(128+e1,a1)==1                                             %距离模糊
           j1(a1,e1)=0;
        else
            j1(a1,e1)=(B(128+e1,a1)-2)/2-32;                               %正常情况，反射率因子以dbz为单位，换算公式参见雷达基数据说明
        end
        f1=f1+1;
    end
    end
end
%提取每1个仰角序数的反射率因子，共有11个仰角序数
n=2                                                                        %第2个仰角序数,代表第1个仰角
for a2=1:4032
    if B(45,a2)>0                                                          %如果每条径向内的仰角数大于0，循环终止                                       
        break
    end
end
for a3=a2:4032                                                             %如果每条径向内的仰角数大于2，循环终止
    if B(45,a3)>2
        break
    end
end
%由前面的仰角序数确认坐标提取范围，对于第1个仰角，仰角序数的范围为1-2
g1=g1(a2:(a3-1),:);%g1=g1(:,1:(B(55,a2)+256*B(56,a2)));仰角坐标
h1=h1(a2:(a3-1),:);%h1=h1(:,1:(B(55,a2)+256*B(56,a2)));方位角坐标
i1=i1(a2:(a3-1),:);%i1=i1(:,1:(B(55,a2)+256*B(56,a2)));距离坐标
j1=j1(a2:(a3-1),:);%j1=j1(:,1:(B(55,a2)+256*B(56,a2)));反射率

[X,Y,Z]=sph2cart((pi/2-h1),g1,i1);                                         %将球面坐标转换成笛卡尔坐标
%设置雷达色标，与PUP产品色标对应，13个色标等级
mycolor2=[255,255,255;123,104,238;
    0,0,255;152,251,152;
    50,205,50;0,100,0;
    255,255,0;192,192,0;
    184,134,11;255,192,192;
    255,0,0;192,0,0;128,0,0]/255.;
colormap(mycolor2);
pcolor(X,Y,j1);
shading interp;
caxis([0,65]);                                                             %色标的范围0-65dbz，分为13级，每级间隔5dbz
colorbar('EastOutside');
    axis image;
    %将反射率因子按坐标输出为文本文件,坐标为相对于雷达的坐标
   [m1,n1]=size(X);
    Xcor=reshape(X,m1*n1,1);
    Ycor=reshape(Y,m1*n1,1);
    j1cor=reshape(j1,m1*n1,1);
    r_1(:,1)=Xcor;
    r_1(:,2)=Ycor;
    r_1(:,3)=j1cor;
fid=fopen('C:\Users\Rabia Lee\Desktop\Z_RADR_I_Z9737_20200820093300_lev1.txt','wt');                   %将重新排列后的数据另存为txt文件
[m2,n2]=size(r_1);
for i=1:1:m2
    for j=1:1:n2
       if j==n2
         fprintf(fid,'%6.3f\n',r_1(i,j));                                     %数值格式设置为双精度，3位小数                          
      else
        fprintf(fid,'%6.3f\t',r_1(i,j));                                      %数值格式设置为双精度，3位小数
       end
    end
end
fclose(fid);
（以上这一部分是读取基数据每一个仰角的反射率因子）

p1=load('F:\雷达回波三维立体数据\Z_RADR_I_Z9737_20200821000500_lev1.txt');
p2=load('F:\雷达回波三维立体数据\Z_RADR_I_Z9737_20200821000500_lev2.txt');
p3=load('F:\雷达回波三维立体数据\Z_RADR_I_Z9737_20200821000500_lev3.txt');
p4=load('F:\雷达回波三维立体数据\Z_RADR_I_Z9737_20200821000500_lev4.txt');
p5=load('F:\雷达回波三维立体数据\Z_RADR_I_Z9737_20200821000500_lev5.txt');
p6=load('F:\雷达回波三维立体数据\Z_RADR_I_Z9737_20200821000500_lev6.txt');
p7=load('F:\雷达回波三维立体数据\Z_RADR_I_Z9737_20200821000500_lev7.txt');
p8=load('F:\雷达回波三维立体数据\Z_RADR_I_Z9737_20200821000500_lev8.txt');
p9=load('F:\雷达回波三维立体数据\Z_RADR_I_Z9737_20200821000500_lev9.txt');
P=[p1;p2;p3;p4;p5;p6;p7;p8;p9];
Point=P(:,1:3);
v=P(:,4);;
F = scatteredInterpolant(Point,v);
x=-100:1:100;
y=-100:1:100;
z=0:0.5:10;
[xq,yq,zq] = meshgrid(x,y,z);
vq = F(xq,yq,zq);
vq=smooth3(vq);
%绘制三维包络面图
myColor=[1 1 1];                                      %地形图色标
mycolor2=[255,255,0;192,192,0;                        %雷达回波色标，根据需要修改
    184,134,11;255,192,192;
    255,0,0;192,0,0;255,0,255]/255.;
mycolor3=[255,255,255;123,104,238;
    0,0,255;152,251,152;
    50,205,50;0,100,0;
    255,255,0;192,192,0;
    184,134,11;255,192,192;
    255,0,0;192,0,0;128,0,0]/255.;
figure(1);
hiso=patch(isosurface(xq,yq,zq,vq,25)); hold on;           %绘制雷达回波三维包络面，包络面值为40
isonormals(xq,yq,zq,vq,hiso);
h=slice(xq,yq,zq,vq,-50,50,2);
colormap(mycolor3);
colorbar('EastOutside');
caxis([0,65]);
shading flat;
hiso.FaceColor='green';                                 %包络面表面颜色设置
hiso.EdgeColor='none';
set(gca,'BoxStyle','full','Box','on')   
view(3);
camlight;
lighting gouraud;
alpha(0.5);
%%
figure(2);
h=slice(xq,yq,zq,vq,-50,50,2);
view(3);
daspect([2 2 40]);                                                         %坐标轴缩放
axis([-100 100 -100 100 0 10]);                                         %坐标轴范围
colormap(mycolor3);
colorbar('EastOutside');
caxis([0,65]);
camlight;
lighting gouraud;
alpha(0.9);
set([h(1),h(2)],'ambientstrength',0.6);
（这一部分是利用读取的9个仰角的反射率因子进行空间插值，然后实现三维可视化，附件图片已经叠加地形，代码中没有叠加地形需要另外自己加入）

simon_han

图也太酷了吧，学习学习

rabialee

F:\雷达回波三维立体数据\三维图\叠加地形的雷达回波切片图.tif

小其其格

我用matlab做过单多普勒和双偏振多普勒雷达的二维显示，三维没做过。
点赞楼主，我学习一下！

小其其格

P=[p1;p2;p3;p4;p5;p6;p7;p8;p9];
Point=P(:,1:3);
v=P(:,4);;

这里P是3维的，v=P(:,4)不会报错么？

rabialee

小其其格 发表于 2020-11-24 21:40
P=;
Point=P(:,1:3);
v=P(:,4);;

第一行P其实是由9页的列数为4的矩阵组成，每1页的1到3列为3维坐标，第4列为反射率因子，不好意思没有把数据图放出来。

yrovl

要不要这么酷炫啊，点赞大神{:eb502:}{:eb502:}

Luoli2195

Cool！！！

小其其格

rabialee 发表于 2020-11-25 17:08
第一行P其实是由9页的列数为4的矩阵组成，每1页的1到3列为3维坐标，第4列为反射率因子，不好意思没有把数 ...

r_1(:,1)=Xcor;
    r_1(:,2)=Ycor;
    r_1(:,3)=j1cor;

那这里少了Zcor的赋值和写入

rabialee

小其其格 发表于 2020-11-26 09:15
r_1(:,1)=Xcor;
    r_1(:,2)=Ycor;
    r_1(:,3)=j1cor;

j1cor就是z方向的坐标