matlab提取温度,科学网—站点气温数据的积温计算（含Matlab程序实现） - 朱永超的博文...

一时失言乱红尘 2023-01-13 10:29 66阅读 0赞

活动温度总和(简称积温)是某一段时间内逐日平均气温≥10℃持续期间日平均气温的总和。是研究温度与生物有机体发育速度之间关系的一种指标，从强度和作用时间两个方面表示温度对生物有机体生长发育的影响。一般以度日(d·℃)为单位。

积温一般使用站点记录的日平均气温进行计算。我国建设有2000余个国家气象站点，最早观测时间至1951。观测数据可在国际气象信息中心申请(中国气象数据网)。这些站点的日值观测数据可以用于积温的计算，图1是站点日尺度气温数据的格式。温度数据共13列，第1列为站点编号，第2、3、4列站点经纬度和海拔高度，第5-7列为观测日期，第8-10列为日平均温、最高温和最低温，11-13列为数据观测状态。

![89dd85de21621b3d6b05cf45c4aada03.png][]

图1. 日观测气温数据格式

积温计算代码

1. 将所有站点数据导入到Matlab中

%% 将所有的温度数据导入到matlab中

%% 选择数据所在的文件夹，获取所有文件信息

clear,clc

% 添加当前路径

addpath(genpath(pwd))

filedir = 'F:\\indices\_wangjing\\tem\\\*txt';

filenames = dir(filedir);

% 文件保存目录

filedir\_save = 'F:\\indices\_wangjing\\tem\_mat\\';

% 获取文件长度

file\_len = length(filenames);

% 逐个文件导入

for ifile = 1:file\_len

disp(ifile/file\_len);

% 当前文件

filename\_this = filenames(ifile).name;

% 获取观测时间

time\_obv = filename\_this(end-9:end-4);

% 读取当前文件(自定义函数参加第3部分)

sm\_singleMonth = importfile\_tem\_v1(filename\_this, 1, inf);

% 纬度转换

lats = sm\_singleMonth(:,2);

lats\_ = num2str(lats);

lats\_D = str2num(lats\_(:,1:2)); % degree

lats\_M = str2num(lats\_(:,3:end)); % minute

lats = lats\_D+lats\_M/60;

% 经度转换

lons = sm\_singleMonth(:,3);

lons\_ = num2str(lons);

lons\_D = str2num(lons\_(:,1:3)); % degree

lons\_M = str2num(lons\_(:,4:end)); % minute

lons = lons\_D+lons\_M/60;

% 转化为日期

date = datetime(sm\_singleMonth(:,5:7));

% 平均温，最高温和最低温

temAvg = sm\_singleMonth(:,8);

temMax = sm\_singleMonth(:,9);

temMin = sm\_singleMonth(:,10);

% 整合数据

temthis = \[sm\_singleMonth(:,1),lats,lons,datenum(date),temAvg,temMax,temMin\];

save(\[filedir\_save,'TEM',time\_obv,'.mat'\],'temthis')

end

2. 计算日平均温大于10度的天数和总积温(1980-2016)

clear, clcaddpath(genpath(pwd))

% load所有的气象观测站信息，包含站点编号、精度和纬度

load('station\_info.mat')stalen = length(stainfo);

% 预定于数据，用于存贮结果

datesGT10 = nan(2419,37);

accumulatedT = nan(2419,37);

allyears = 1980:1:2016;

for iyear = 1:37

tic

% this year

disp(iyear);

thisyear = allyears(iyear);

dateofYear = 365; % 1年中的天数

if rem(thisyear,4)==0

dateofYear = 366;

end

% load所有的温度数据

files = dir(\['G:\\indices\_wangjing\\tem\_mat\\TEM' num2str(thisyear) '\*.mat'\]);

% 将一年中的数据连在一起

temEachyear = \[\];

for imonth = 1:12

% load each month data

load(files(imonth).name)

% cat each month

temEachyear = cat(1,temEachyear,temthis);

end

%% 计算各站点的积温天数和总积温

for ista=1:2419

staNO = stainfo(ista,1);

temeachstation = temEachyear(temEachyear(:,1)==staNO,:);

% temperature=temperature/10

temeachstation(:,5:7) = temeachstation(:,5:7)/10;

%% 剔除异常值

% if the temperature gt 1000, turn in to nan

temeachstation(temeachstation(:,5)>10000,:)=nan;

% delete nan data

temeachstation = temeachstation(~isnan(temeachstation(:,5)),:);

tem\_avg = temeachstation(:,5);

dates\_mea = temeachstation(:,4);

% 如果这个站点一年中的观测数据少于300天，则不计算积温

lendays = length(temeachstation);

if lendays<300

continue;

end

datesta = datenum(thisyear,1,1);

dateend = datenum(thisyear,12,31);

% 如果观测数据少于该年的天数，则对观测数据进行插值

if lendays

tem\_interp = interp1(dates\_mea,tem\_avg,(datesta:1:dateend));

else

tem\_interp = tem\_avg;

end

%% 对插值后的气温数据进行5日滑动平均

for idate = 1:dateofYear

if idate>2 && idate

tem\_interp(idate) = (tem\_interp(idate-2)+tem\_interp(idate-1)+...

tem\_interp(idate)+tem\_interp(idate+1)+tem\_interp(idate+2))/5;

else

continue

end

%% ！！！积温计算算法

% 找到日平均气温大于10的日期

Iaa = find(tem\_interp>10);

len1 = zeros(366,1);

j=1;

for i=1:length(Iaa)-1

if Iaa(i+1)-Iaa(i)==1

len1(j) = len1(j)+1;

else

j=j+1;

end

% a代表气温连续超过10的最长天数

\[a,b\] = max(len1);

c = sum(len1(1:b-1));

% if no temperature gt 10, continue

if a<1

continue

end

% d代表积温计算的起始日期

d=Iaa(b+c);

% 保存数据

datesGT10(ista,iyear) = a;

accumulatedT(ista,iyear) = sum(tem\_interp(d:d+a),'omitnan');

end

toc

end

%% 将结果保存到excel表格中

xlswrite('大于10度的积温天数.xlsx',datesGT10)

xlswrite('大于10度的年积温.xlsx',accumulatedT)

3.自定义函数importfile\_tem\_v1.m

function temthis = importfile\_tem\_v1(filename, startRow, endRow)

%IMPORTFILE 将文本文件中的数值数据作为矩阵导入。

%   SURFCLICHNMULDAYTEM12001198001 = IMPORTFILE(FILENAME) 读取文本文件 FILENAME

%   中默认选定范围的数据。

%   SURFCLICHNMULDAYTEM12001198001 = IMPORTFILE(FILENAME, STARTROW, ENDROW)

%   读取文本文件 FILENAME 的 STARTROW 行到 ENDROW 行中的数据。

% Example:

%   temthis = importfile('SURF\_CLI\_CHN\_MUL\_DAY-TEM-12001-198001.TXT', 1, 74679);

%    另请参阅 TEXTSCAN。

% 由 MATLAB 自动生成于 2019/01/12 15:20:54

%% 初始化变量。

if nargin<=2

startRow = 1;

endRow = inf;

end

%% 每个文本行的格式字符串:

%   列1: 双精度值 (%f)

%   列2: 双精度值 (%f)

%   列3: 双精度值 (%f)

%   列4: 双精度值 (%f)

%   列5: 双精度值 (%f)

%   列6: 双精度值 (%f)

%   列7: 双精度值 (%f)

%   列8: 双精度值 (%f)

%   列9: 双精度值 (%f)

%   列10: 双精度值 (%f)

% 有关详细信息，请参阅 TEXTSCAN 文档。

formatSpec = '%5f%5f%6f%7f%5f%3f%3f%7f%7f%7f%\[^\\n\\r\]';

%% 打开文本文件。

fileID = fopen(filename,'r');

%% 根据格式字符串读取数据列。

% 该调用基于生成此代码所用的文件的结构。如果其他文件出现错误，请尝试通过导入工具重新生成代码。

dataArray = textscan(fileID, formatSpec, endRow(1)-startRow(1)+1, 'Delimiter', '', 'WhiteSpace', '', 'EmptyValue' ,NaN,'HeaderLines', startRow(1)-1, 'ReturnOnError', false);

for block=2:length(startRow)

frewind(fileID);

dataArrayBlock = textscan(fileID, formatSpec, endRow(block)-startRow(block)+1, 'Delimiter', '', 'WhiteSpace', '', 'EmptyValue' ,NaN,'HeaderLines', startRow(block)-1, 'ReturnOnError', false);

for col=1:length(dataArray)

dataArray\{col\} = \[dataArray\{col\};dataArrayBlock\{col\}\];

end

%% 关闭文本文件。

fclose(fileID);

%% 对无法导入的数据进行的后处理。

% 在导入过程中未应用无法导入的数据的规则，因此不包括后处理代码。要生成适用于无法导入的数据的代码，请在文件中选择无法导入的元胞，然后重新生成脚本。

%% 创建输出变量

temthis= \[dataArray\{1:end-1\}\];

转载本文请联系原作者获取授权，同时请注明本文来自朱永超科学网博客。

链接地址：http://blog.sciencenet.cn/blog-3367669-1161863.html

上一篇：风云3C土壤水分产品在河南省农业区的验证