- 积分
- 493
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2013-4-2
- 最后登录
- 1970-1-1
|
20金钱
matlab画FY2G降水估计6小时,图1降水落区和量级都没问题放大就是图2的样子,这图要不了。
求解,附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
云顶温度同样的情况,云南省附件云顶温度对比国家卫星中心,我画的偏东,偏南,放大都是小格子。
m文件基本与降水估计一样,
data1(data1>50)=nan;%仅取小于50的部分
问过国家卫星中心说0-10是表征量,250表空,把250踢了画出来,都是小格子,
具体k值不会换,卫星中心说是卫星内部编码,让写信去到卫星中心问,写去了也一直内回音。
求救!
|
|