此篇介绍如何爬取中国大陆的探空数据，并转成.FSL格式，或三捷AERMOD需要的.rao文件

最新推荐文章于 2024-09-25 11:45:01 发布

weixin_38999566

最新推荐文章于 2024-09-25 11:45:01 发布

阅读量1.2k

点赞数 30

文章标签： matlab

本文链接：https://blog.csdn.net/weixin_38999566/article/details/140084480

版权

数据网址：

University of Wyoming Atmospheric Science Radiosonde Archive

手动点击下载，只支持一次下载一个站一天中一次的数据。

2.中国大陆探空站信息

国家气象信息中心-中国气象数据网

可以查到所有站的站号、站名、经纬度、海拔高度（不是拔海）

3. 需转成的.FSL格式说明

每个数据块均包括控制信息和数据信息

（1）控制信息前4行

注1：LINTYP（在控制信息中LINTYP只能依次等于254、1、2、3）

254=新的观测数据开始，该行写入时间信息

1 =站点标识栏，该行写入站点标识信息

2 =观测检测栏 sounding checks line

3 = 站点识别或者其他标识station identifier and other indicators line

4 = 强制层mandatory level

5 = 重要层significant level

6 = wind level (PPBB) (GTS or merged data)

7 = tropopause level (GTS or merged data)

8 = maximum wind level (GTS or merged data)

9 = 地表面层，仅第一层写入此数，且必须写入。

注2：月份为JAN、FEB、MAR、APR、MAY、JUN、JUL、AUG、SEP、OCT、NOV、DEC。

注3：风速标识，可输入米/秒，也可 kt（海里/小时）。本代码按kt输入

（2）数据信息

注意单位转换 .FSL只用到下载数据列中涂黄色的。

下载数据块的排列和单位。

PRES HGHT TEMP DWPT RELH MIXR DRCT SPED THTA THTE THTV

hPa m C C % g/kg deg m/s K K K

.FSL中的数据列和单位。（注意表中的气压单位不是0.1百帕，而是百帕）

注1：各要素均为整数，不足位数，高位补空格

注2：温度和露点温度均有正负标识

注3：写入时风速单位固定为米/秒

4. 分享一个可直接实现，自动爬取探高数据，并转成.FSL格式，或三捷AERMOD需要的.rao文件的MATLAB M文件代码。用MATLAB打开M文件后，修改必须的站号、经纬度、海拔高度、起始时间后，直接点M文件的运行即可。生成的文件在当前目录文件夹下。

function upperair_download()
clear
clc
filename = 'upperair.txt'; % 生成的FSL文件，运行完后将后缀修改成 .rao
% 1. 可不修改，也可根据自身情况修改 &&&&&&@@@@@@****** 生成的FSL文件，运行完后将后缀修改成 .rao
fid = fopen(filename, 'a');
date_xulie= datetime(2023,01,01):caldays(1):datetime(2023,01,02); % 前为起始日期后为结束日期
% 2. 根据自身情况修改 &&&&&&@@@@@@****** 日期不一样，修改此处
id='56187';
% 3. 根据自身情况修改 &&&&&&@@@@@@****** 探空站编号
Latitude='30.75N';
Longitude='103.87E';
% 4. 根据自身情况修改 &&&&&&@@@@@@****** 探空站经纬度 ,小数位数保持一致
haiba='548';
% 5. 根据自身情况修改 &&&&&&@@@@@@****** 探空站海拔高度，整数，不带小数
[n,num]=size(date_xulie);
for i=1:num



cengshu=0;

tempstr=datestr(date_xulie(i));
tempstr=tempstr(1:11);
switch tempstr(4:6)
case 'Jan'
month='01';
case 'Feb'
month='02';
case 'Mar'
month='03';
case 'Apr'
month='04';
case 'May'
month='05';
case 'Jun'
month='06';
case 'Jul'
month='07';
case 'Aug'
month='08';
case 'Sep'
month='09';
case 'Oct'
month='10';
case 'Nov'
month='11';
case 'Dec'
month='12';
end

datetime1=strcat(tempstr(8:11),'-',month,'-',tempstr(1:2));
url1=strcat( 'http://weather.uwyo.edu/cgi-bin/bufrraob.py?src=bufr&datetime=',datetime1);
url1_last=strcat( url1,'%2000:00:00&id=',id,'&type=TEXT:LIST');
url2_last=strcat( url1,'%2012:00:00&id=',id,'&type=TEXT:LIST');
options = weboptions('CharacterEncoding','auto','Timeout',600,'RequestMethod','post');
date_num=num2str(str2num(tempstr(1:2)));
[long,long1]=size(date_num);
if long1==2
date_num=strcat(32,32,32,32,32,date_num);
else
date_num=strcat(32,32,32,32,32,32,date_num);
end
month_capital=upper(tempstr(4:6));


try
content1 = webread(url1_last,options);

positions1 = strfind(content1, '<PRE>');
positions2 = strfind(content1, '</PRE>');
content11=content1(positions1+5+314:positions2-1);



[long,long2]=size(content11);
if rem(long2,78)==0
hang=floor(long2/78);
else
hang=floor(long2/78)+1;
end

for j=1:hang-1

PRES_temp=content11((j-1)*78+1:(j-1)*78+7);
HGHT_temp=content11((j-1)*78+8:(j-1)*78+14);
TEMP_temp=content11((j-1)*78+15:(j-1)*78+21);
DWPT_temp=content11((j-1)*78+22:(j-1)*78+28);
DRCT_temp=content11((j-1)*78+43:(j-1)*78+49);
SPED_temp=content11((j-1)*78+50:(j-1)*78+56);

if PRES_temp== ' '
PRES_temp='99999';
num2=5;
else
PRES_temp=round(str2num(PRES_temp));
if PRES_temp<1000
PRES_temp=strcat(32,32,32,32,num2str(PRES_temp));
else
PRES_temp=strcat(32,32,32,num2str(PRES_temp));

end

num2=7;
end
switch num2
case 5
fprintf(fid_temp, ' ');
fprintf(fid_temp, PRES_temp);
case 7
fprintf(fid_temp, PRES_temp);
end

if HGHT_temp== ' '

HGHT_temp='32727';
num2=5;
else
num2=7;
end
switch num2
case 5
fprintf(fid_temp, ' ');
fprintf(fid_temp, HGHT_temp);
case 7
fprintf(fid_temp, HGHT_temp);
end

[num1,num2]=size(str2num(TEMP_temp));
if num2
TEMP_temp=str2num(TEMP_temp)*10;
TEMP_temp = round(TEMP_temp, 1);

TEMP_temp=num2str(TEMP_temp);
[num1,num2]=size(TEMP_temp);
else
TEMP_temp='99999';
num2=5;
end

switch num2
case 1
fprintf(fid_temp, ' ');
fprintf(fid_temp, TEMP_temp);
case 3
fprintf(fid_temp, ' ');
fprintf(fid_temp, TEMP_temp);
case 2
fprintf(fid_temp, ' ');
fprintf(fid_temp, TEMP_temp);
case 3
fprintf(fid_temp, ' ');
fprintf(fid_temp, TEMP_temp);
case 4
fprintf(fid_temp, ' ');
fprintf(fid_temp, TEMP_temp);
case 5
fprintf(fid_temp, ' ');
fprintf(fid_temp, TEMP_temp);
case 6
fprintf(fid_temp, ' ');
fprintf(fid_temp, TEMP_temp);
case 7
fprintf(fid_temp, TEMP_temp);
end

[num1,num2]=size(str2num(DWPT_temp));
if num2
DWPT_temp=str2num(DWPT_temp)*10;
DWPT_temp = round(DWPT_temp, 1);
DWPT_temp=num2str(DWPT_temp);
[num1,num2]=size(DWPT_temp);
else
DWPT_temp='32727';
num2=5;
end

switch num2
case 1
fprintf(fid_temp, ' ');
fprintf(fid_temp, DWPT_temp);
case 2
fprintf(fid_temp, ' ');
fprintf(fid_temp, DWPT_temp);
case 3
fprintf(fid_temp, ' ');
fprintf(fid_temp, DWPT_temp);
case 4
fprintf(fid_temp, ' ');
fprintf(fid_temp, DWPT_temp);
case 5
fprintf(fid_temp, ' ');
fprintf(fid_temp, DWPT_temp);
case 6
fprintf(fid_temp, ' ');
fprintf(fid_temp, DWPT_temp);
case 7
fprintf(fid_temp, DWPT_temp);
end

if DRCT_temp== ' '

DRCT_temp='99999';
num2=5;
else
num2=7;
end
switch num2
case 5
fprintf(fid_temp, ' ');
fprintf(fid_temp, DRCT_temp);
case 7
fprintf(fid_temp, DRCT_temp);
end

[num1,num2]=size(str2num(SPED_temp));
if num2
SPED_temp=str2num(SPED_temp)*1.9440;
SPED_temp = round(SPED_temp);
SPED_temp=num2str(SPED_temp);
[num1,num2]=size(SPED_temp);
else
SPED_temp='99999';
num2=5;
end

switch num2
case 1
fprintf(fid_temp, ' ');
fprintf(fid_temp, SPED_temp);
case 2
fprintf(fid_temp, ' ');
fprintf(fid_temp, SPED_temp);
case 3
fprintf(fid_temp, ' ');
fprintf(fid_temp, SPED_temp);
case 4
fprintf(fid_temp, ' ');
fprintf(fid_temp, SPED_temp);
case 5
fprintf(fid_temp, ' ');
fprintf(fid_temp, SPED_temp);
case 6
fprintf(fid_temp, ' ');
fprintf(fid_temp, SPED_temp);
case 7
fprintf(fid_temp, SPED_temp);
end
fprintf(fid_temp, '\n');
end
end
fclose(fid_temp);
cengshu=cengshu+4;

if cengshu>0 & cengshu<10
buge=6;
end
if cengshu>=10 & cengshu<100
buge=5;
end
if cengshu>=100 & cengshu<1000
buge=4;
end



head11=strcat(' 254',' 0',date_num,32,32,32,32,32,32,month_capital,32,32,32,32,tempstr(8:11));
if haiba<10000
head12=strcat(' 1',' 3190',32,32,id,32,32,Latitude,Longitude,32,32,32,haiba,' 2303');%% 3190 此处随便4位？？
elseif haiba<1000
head12=strcat(' 1',' 3190',32,32,id,32,32,Latitude,Longitude,32,32,32,32,haiba,' 2303');%% 3190 此处随便4位？？
elseif haiba<100
head12=strcat(' 1',' 3190',32,32,id,32,32,Latitude,Longitude,32,32,32,32,32,haiba,' 2303');%% 3190 此处随便4位？？
elseif haiba<10
head12=strcat(' 1',' 3190',32,32,id,32,32,Latitude,Longitude,32,32,32,32,32,32,haiba,' 2303');%% 3190 此处随便4位？？
end


if buge==4
   head13=strcat(' 2 100 100 194',32,32,32,32,num2str(cengshu),' 32767 3');%% '44'为层数*****
elseif buge==5
head13=strcat(' 2 100 100 194',32,32,32,32,32,num2str(cengshu),' 32767 3');%% '44'为层数*****
elseif buge==6
head13=strcat(' 2 100 100 194',32,32,32,32,32,32,num2str(cengshu),' 32767 3');%% '44'为层数*****
end

head14=' 3 32767 kt';
head24=' 3 32767 kt';


fprintf(fid, head11);
fprintf(fid, '\n');
fprintf(fid, head12);
fprintf(fid, '\n');
fprintf(fid, head13);
fprintf(fid, '\n');
fprintf(fid, head14);
fprintf(fid, '\n');


fid_temp = fopen('upperair_temp.txt', 'r');

while ~feof(fid_temp)
tline=fgetl(fid_temp);
fprintf(fid, '%s', tline);
fprintf(fid, '\n');
end
fclose(fid_temp);


% file1 = fgetl(fid_temp);

% fprintf(fid, '%s', file1);
% fprintf(fid, '\n');
% ii = 1:length(file1)
% for ii = 1:length(file1)
% tline =) fgetl(fid_temp;
% fprintf(fid, '%s', tline);
% fprintf(fid, '\n');
% end

catch err
disp('网络出错啦')
disp(tempstr)
end













filename_temp = 'upperair_temp.txt'; % 生成的FSL文件名
% 1. 根据自身情况修改 &&&&&&@@@@@@******
fid_temp = fopen(filename_temp, 'w');

cengshu=0;


try
content2 = webread(url2_last,options);
positions1 = strfind(content2, '<PRE>');
positions2 = strfind(content2, '</PRE>');
content21=content2(positions1+5+314:positions2-1);

date_num=num2str(str2num(tempstr(1:2)));
[long,long1]=size(date_num);


[long,long2]=size(content21);
if rem(long2,78)==0
hang=floor(long2/78);
else
hang=floor(long2/78)+1;
end



if long1==2
date_num=strcat(32,32,32,32,32,date_num);
else
date_num=strcat(32,32,32,32,32,32,date_num);
end
month_capital=upper(tempstr(4:6));

for j=1:hang-1

PRES_temp=content21((j-1)*78+1:(j-1)*78+7);
HGHT_temp=content21((j-1)*78+8:(j-1)*78+14);
TEMP_temp=content21((j-1)*78+15:(j-1)*78+21);
DWPT_temp=content21((j-1)*78+22:(j-1)*78+28);
DRCT_temp=content21((j-1)*78+43:(j-1)*78+49);
SPED_temp=content21((j-1)*78+50:(j-1)*78+56);

end

num2=7;
end
switch num2
case 5
fprintf(fid_temp, ' ');
fprintf(fid_temp, PRES_temp);
case 7
fprintf(fid_temp, PRES_temp);
end

if HGHT_temp== ' '

HGHT_temp='32727';
num2=5;
else
num2=7;
end
switch num2
case 5
fprintf(fid_temp, ' ');
fprintf(fid_temp, HGHT_temp);
case 7
fprintf(fid_temp, HGHT_temp);
end

[num1,num2]=size(str2num(TEMP_temp));
if num2
TEMP_temp=str2num(TEMP_temp)*10;
TEMP_temp = round(TEMP_temp, 1);

TEMP_temp=num2str(TEMP_temp);
[num1,num2]=size(TEMP_temp);
else
TEMP_temp='99999';
num2=5;
end

if DRCT_temp== ' '

DRCT_temp='99999';
num2=5;
else
num2=7;
end
switch num2
case 5
fprintf(fid_temp, ' ');
fprintf(fid_temp, DRCT_temp);
case 7
fprintf(fid_temp, DRCT_temp);
end

head21=strcat(' 254',' 12',date_num,32,32,32,32,32,32,month_capital,32,32,32,32,tempstr(8:11));

if haiba<10000
head22=strcat(' 1',' 3190',32,32,id,32,32,Latitude,Longitude,32,32,32,haiba,' 2303');%% 3190 此处随便4位？？
elseif haiba<1000
head22=strcat(' 1',' 3190',32,32,id,32,32,Latitude,Longitude,32,32,32,32,haiba,' 2303');%% 3190 此处随便4位？？
elseif haiba<100
head22=strcat(' 1',' 3190',32,32,id,32,32,Latitude,Longitude,32,32,32,32,32,haiba,' 2303');%% 3190 此处随便4位？？
elseif haiba<10
head22=strcat(' 1',' 3190',32,32,id,32,32,Latitude,Longitude,32,32,32,32,32,32,haiba,' 2303');%% 3190 此处随便4位？？
end


if buge==4
   head23=strcat(' 2 100 100 194',32,32,32,32,num2str(cengshu),' 32767 3');%% '44'为层数*****
elseif buge==5
   head23=strcat(' 2 100 100 194',32,32,32,32,32,num2str(cengshu),' 32767 3');%% '44'为层数*****
elseif buge==6
   head23=strcat(' 2 100 100 194',32,32,32,32,32,32,num2str(cengshu),' 32767 3');%% '44'为层数*****
end

head24=' 3 32767 kt';




fprintf(fid, head21);
fprintf(fid, '\n');
fprintf(fid, head22);
fprintf(fid, '\n');
fprintf(fid, head23);
fprintf(fid, '\n');
fprintf(fid, head24);
fprintf(fid, '\n');


fid_temp = fopen('upperair_temp.txt', 'r');
fid_temp = fopen('upperair_temp.txt', 'r');
while ~feof(fid_temp)
tline=fgetl(fid_temp);
fprintf(fid, '%s', tline);
fprintf(fid, '\n');
end
catch err
disp('网络出错啦')
disp(tempstr)
end
fclose(fid_temp);