|||
活动温度总和(简称积温)是某一段时间内逐日平均气温≥10℃持续期间日平均气温的总和。是研究温度与生物有机体发育速度之间关系的一种指标,从强度和作用时间两个方面表示温度对生物有机体生长发育的影响。一般以度日(d·℃)为单位。
积温一般使用站点记录的日平均气温进行计算。我国建设有2000余个国家气象站点,最早观测时间至1951。观测数据可在国际气象信息中心申请(中国气象数据网)。这些站点的日值观测数据可以用于积温的计算,图1是站点日尺度气温数据的格式。温度数据共13列,第1列为站点编号,第2、3、4列站点经纬度和海拔高度,第5-7列为观测日期,第8-10列为日平均温、最高温和最低温,11-13列为数据观测状态。
图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<dateofYear
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<dateofYear-1
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
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
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
end
%% 关闭文本文件。
fclose(fileID);
%% 对无法导入的数据进行的后处理。
% 在导入过程中未应用无法导入的数据的规则,因此不包括后处理代码。要生成适用于无法导入的数据的代码,请在文件中选择无法导入的元胞,然后重新生成脚本。
%% 创建输出变量
temthis = [dataArray{1:end-1}];
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-11-23 09:29
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社