爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 701|回复: 2

[程序设计] matlab画FY2G云顶温度和降水估计出现的问题

[复制链接]

新浪微博达人勋

发表于 2024-3-7 15:00:46 | 显示全部楼层 |阅读模式
20金钱
matlab画FY2G降水估计6小时,图1降水落区和量级都没问题放大就是图2的样子,这图要不了。 2022-05-20 08.png 2022-05-21 08.png
求解,附m文件:
clear all;clc;close all;
lj='E:\11\pre-6\22\';%数据存放路径
file=dir([lj '*.AWX']);%读数据
No_file=length(file);%awx文件数量
for i=1:No_file
     a=file(i).name;
    fid = fopen([lj a],'rb');
    name=fread(fid,12,'char');%文件名
    int_order=fread(fid,1,'uint16');%整型数的字节顺序
    head1_length=fread(fid,1,'uint16');%第一级文件头长度
    head2_length=fread(fid,1,'uint16');%第二级文件头长度
    tianchongdata_length=fread(fid,1,'uint16');%填充段数据长度
    record_length=fread(fid,1,'uint16');%记录长度
    head_record=fread(fid,1,'uint16');%文件头占用记录数
    data_record=fread(fid,1,'uint16');%产品数据占用记录数
    product_kind=fread(fid,1,'uint16');%产品类别
    zip_method=fread(fid,1,'uint16');%压缩方式
    form_string=fread(fid,8,'char');%格式说明字串
    data_quality=fread(fid,1,'uint16');%产品数据质量标记
    satellite_name=fread(fid,8,'char');%卫星名
    grid_factor=fread(fid,1,'uint16');%格点场要素
    grid_byte=fread(fid,1,'uint16');%格点数据字节
    grid_base=fread(fid,1,'uint16');%格点数据基准值
    grid_scale=fread(fid,1,'uint16');%格点数据比例因子
    time_scale=fread(fid,1,'uint16');%时间范围代码
    start_year=fread(fid,1,'uint16');%开始年
    start_month=fread(fid,1,'uint16');%开始月
    start_day=fread(fid,1,'uint16');%开始日
    start_hour=fread(fid,1,'uint16');%开始时
    start_minute=fread(fid,1,'uint16');%开始分
    end_year=fread(fid,1,'uint16');%结束年
    end_month=fread(fid,1,'uint16');%结束月
    end_day=fread(fid,1,'uint16');%结束日
    end_hour=fread(fid,1,'uint16');%结束时
    end_minute=fread(fid,1,'uint16');%结束分
    leftup_lat=fread(fid,1,'uint16');%网格左上角纬度
    leftup_lon=fread(fid,1,'uint16');%网格左上角经度
    rightdown_lat=fread(fid,1,'int16');%网格右下角纬度
    rightdown_lon=fread(fid,1,'uint16');%网格右下角经度
    grid_unit=fread(fid,1,'uint16');%格距单位
    heng_unit=fread(fid,1,'uint16');%横向格距
    shu_unit=fread(fid,1,'uint16');%纵向格距
    heng_num=fread(fid,1,'uint16');%横向格点数
    shu_num=fread(fid,1,'uint16');%纵向格点数
    is_land=fread(fid,1,'uint16');%有无陆地判释值
    land_juti=fread(fid,1,'uint16');%陆地具体判释值
    is_cloud=fread(fid,1,'uint16');%有无云判释值
    cloud_juti=fread(fid,1,'uint16');%云具体判释值
    is_water=fread(fid,1,'uint16');%有无水体判释值
    water_juti=fread(fid,1,'uint16');%水体具体判释值
    is_ice=fread(fid,1,'uint16');%有无冰体判释值
    ice_juti=fread(fid,1,'uint16');%冰体具体判释值
    is_quility=fread(fid,1,'uint16');%是否有质量控制值
    quility_up=fread(fid,1,'uint16');%质量控制值上限
    quility_down=fread(fid,1,'uint16');%质量控制值下限
    backup=fread(fid,1,'uint16');%备用
    tianchongduan=fread(fid,128,'uint8');%填充段
    tianchong_kuozhan=fread(fid,228,'uint8');%填充段扩展
    data0 = fread(fid,[heng_num shu_num]);% 读取pre-6小时数据
    data0 = data0';% 二进制读取读取之后需要进行一个转置操作,才能是正常经纬度放置方式
    data1=[data0(:,726:end),data0(:,1:725)];
    fclose(fid); % 关闭文件
    data1(data1>140)=nan;%仅取降水不超过140mm的部分

    %% 下面是画图部分
    latmax=leftup_lat/100;
    lonmin=leftup_lon/100;
    latmin=rightdown_lat/100;
    lonmax=rightdown_lon/100;
    f=figure('units','normalized','position',[0.2 0.1 0.6 0.8],'visible','on'); % 新建图形窗口,设为不可见
    %set(f,'position',[200 100 1000 800]);
    m_proj('miller','lat',[10 35],'lon',[85 120]);
    [LON,LAT]=meshgrid(lonmin:0.1:lonmax,latmax:-0.1:latmin);
    m_contourf(LON,LAT,data1,'LineColor','none');
    colormap(flipud(jet(800)));
    shading interp;
    h=colorbar('eastoutside');
    M=m_shaperead('E:\11\国家基础地理信息系统数据\国界与省界\bou2_4l');
    for k=1:length(M.ncst)
        AA1=find(M.ncst{k}(:,1)>0);BB1=find(M.ncst{k}(:,2)<18888);
        AB1=union(AA1,BB1);
        m_line(M.ncst{k}(AB1,1),M.ncst{k}(AB1,2),'Color','k');
        hold on
    end
    %添加其它国家地图
    M=m_shaperead('E:\11\中国为中心的世界地图shp\中国为中心的世界地图shp\世界地图shp\世界地图shp格式\世界各国和主要城市shp图\世界地图\country');
    for k=1:length(M.ncst)
        m_line(M.ncst{k}(:,1),M.ncst{k}(:,2),'Color','k');
        hold on
    end

    bt=char(a(49:61));%降水的时间选取
    tt=datetime(str2num(bt(1:4)),str2num(bt(5:6)),str2num(bt(7:8)),str2num(bt(10:11)),str2num(bt(12:13)),0);%时间选取
   t_BJ = tt+8/24;%换北京时
   t_BJ.Format='yyyy-MM-dd HH:mm:ss';
   at=char(t_BJ);
   at1=char(at(1:13));
   
    savepath='E:\11\pre-6\1\';
    title([char(t_BJ) ' FY2G PRE-6h'],'FontName','宋体','FontSize',16);
    text(0.57,1.1,'单位:℃','FontName','宋体','FontSize',10);
    xlabel('');
    ylabel('');
    m_grid;
    %设置图例颜色
    z_max=nanmax(max(data1));
    z_min=nanmin(min(data1));
   color_data=[];
    cmin=z_min;
    cmax=z_max;
    caxis([cmin cmax]);
    color_map=cmin:(cmax-cmin)/64:cmax;
    size_color=size(color_map);
      sb=[166 242 143;61 186 61;97 184 255;0 0 255;...
        250 0 250;128 0 64]./255;
    for i=1:(size_color(1,2)-1)
        if color_map(1,i)<=3.9
            color_data(i,:)=sb(1,:);
        elseif color_map(1,i)>3.9&color_map(1,i)<=12.9
            color_data(i,:)=sb(2,:);
        elseif color_map(1,i)>12.9&color_map(1,i)<=24.9
            color_data(i,:)=sb(3,:);
        elseif color_map(1,i)>24.9&color_map(1,i)<=59.9
            color_data(i,:)=sb(4,:);
        elseif color_map(1,i)>59.9&color_map(1,i)<=119.9
            color_data(i,:)=sb(5,:);
            else
            color_data(i,:)=sb(6,:);
     end
    end
    colormap(color_data);
    h=colorbar('eastoutside');
      saveas(gca,[savepath at1,'.png']);
   end

云顶温度同样的情况,云南省附件云顶温度对比国家卫星中心,我画的偏东,偏南,放大都是小格子。
国家卫星中心原图.jpg 2024-02-28 08.png
m文件基本与降水估计一样,
data1(data1>50)=nan;%仅取小于50的部分

问过国家卫星中心说0-10是表征量,250表空,把250踢了画出来,都是小格子,
具体k值不会换,卫星中心说是卫星内部编码,让写信去到卫星中心问,写去了也一直内回音。
求救!





密码修改失败请联系微信:mofangbao

新浪微博达人勋

 成长值: 32430
发表于 2024-3-7 23:46:05 | 显示全部楼层
本帖最后由 二爷名声在外 于 2024-3-7 23:52 编辑

感觉好像是用我读TBB的程序改的,不知道也没有看格式说明。方便的话加我vx,发我文件和对应的格式说明看看。
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2024-4-14 16:40:09 | 显示全部楼层
不错的帖子,收藏一下,可以好好学一学
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表