%% metwrite.m % making a check of meteorological data from Tarfala weather station % peter.jansson@2000-05-07 -> 2002-04-21 %% Define year and read data % Assembling this data is done in two parts and by inserting columns of % NaNs since the number of data columns change in September. Year = 1989; hourly = dlmread('TRSmet1989.csv'); indnan = hourly==-999; hourly(indnan) = NaN; clear indnan; indnan = hourly==6999; hourly(indnan) = NaN; clear indnan; indnan = hourly==-6999; hourly(indnan) = NaN; clear indnan; % Define columns of data to convert to a date string. % Note dummy seconds variable Hyr = hourly(:,1); Hjd = hourly(:,2); Hhhmm = hourly(:,3); Hs = zeros(length(hourly(:,1)),1); %Dyr = daily(:,2); %Djd = daily(:,3); %Dhhmm = daily(:,4); %Dyr = daily(:,2); %Djd = floor(daily(:,3)); %Dhhmm = zeros(length(daily(:,1)),1); %Ds = zeros(length(daily(:,1)),1); % Call the string funmction. [HourlyString] = jday2date(Hyr,Hjd,Hhhmm,Hs); %[DailyString] = jday2date(Dyr,Djd,Dhhmm,Ds); [nrows, ncols] = size(hourly); % Open files for printing fidT=fopen(['TRS_met_',num2str(Year),'_Temperature.csv'],'w'); fidRh=fopen(['TRS_met_',num2str(Year),'_Relative_humidity.csv'],'w'); fidW=fopen(['TRS_met_',num2str(Year),'_Wind.csv'],'w'); %fidB=fopen(['TRS_met_',num2str(Year),'_Barometric_pressure.csv'],'w'); fidP=fopen(['TRS_met_',num2str(Year),'_Precipitation.csv'],'w'); fidR=fopen(['TRS_met_',num2str(Year),'_Radiation.csv'],'w'); % print the different files for row=1:nrows fprintf(fidT,'%19s,%.2f\n',HourlyString(row,:),hourly(row,4)); fprintf(fidRh,'%19s,%.1f\n',HourlyString(row,:),hourly(row,5)); fprintf(fidW,'%19s,%.1f,%.1f,%.1f,%.1f\n',HourlyString(row,:),hourly(row,6:9)); % fprintf(fidB,'%19s,%.1f\n',HourlyString(row,:),hourly(row,26)); fprintf(fidP,'%19s,%.2f\n',HourlyString(row,:),hourly(row,13)); fprintf(fidR,'%19s,%.2f,%.2f,%.2f\n',HourlyString(row,:),hourly(row,10:12)); end %[mrows, mcols] = size(daily); % open file for daily output an dprint %fidD=fopen(['TRS_met_',num2str(Year),'_Daily_data.csv'],'w'); %for row=1:mrows % fprintf(fidD,'%19s,%.2f,%.2f,%.2f,%.1f,%.2f,%.0f,%.2f,%.0f,%.1f,%.0f,%.1f,%.1f,%.1f,%.1f,%.2f,%.2f\n',DailyString(row,:),daily(row,3:16)); %end fclose('all'); b=[[num2str(Year),'-01-01']; [num2str(Year),'-02-01']; [num2str(Year),'-03-01']; [num2str(Year),'-04-01']; [num2str(Year),'-05-01']; [num2str(Year),'-06-01']; [num2str(Year),'-07-01']; [num2str(Year),'-08-01']; [num2str(Year),'-09-01']; [num2str(Year),'-10-01']; [num2str(Year),'-11-01']; [num2str(Year),'-12-01']]; format short g H = [datenum(HourlyString) hourly(:,4:14)]; TickNum = datenum(b); TickDate = '1 Jan.|1 Feb.|1 Mar.|1 Apr.|1 May|1 Jun.|1 Jul.|1 Aug.|1 Sep.|1 Oct.|1 Nov.|1 Dec.'; %return StartDate = datenum([num2str(Year),'-01-01']); axismin = StartDate; axismax = datenum([num2str(Year+1),'-01-01']); textmin = axismin + 5; %% Figures % setting up constants for plotting screenrect = get(0,'screensize'); screenwidth = screenrect(3); screenheight = screenrect(4); figwidth = floor(0.9*screenwidth); figheight = floor(0.9*screenheight); figure1=figure('position', [(screenwidth/2-figwidth/2) (screenheight/2-figheight/2) figwidth figheight],... 'units','pixels','NumberTitle','off','Name',['Meteorological data from Tarfala ',num2str(Year)]); orient landscape set(gcf, 'PaperType', 'A4'); n=1; % for plotting selected points, n=4 is pseudo hourly data rows=6; subplot(rows,1,1) plot(H(1:n:end,1),H(1:n:end,2),'b',... [axismin axismax],[0 0],'k') axis([axismin axismax -20 20]); text(textmin,15,'Air temperature'); legend('Young (T/Rh)','Location','southwest') ylabel('(^\circ C)'); title(['Tarfala Research Station meteorological data for ',num2str(Year)]) set(gca,'XTick',[]) %return subplot(rows,1,2) plot(H(1:n:end,1),H(1:n:end,3),'b') axis([axismin axismax 0 100]); text(textmin,70,'Relative humidity'); legend('Toung (T/Rh)','Location','southwest') ylabel('(%)'); set(gca,'XTick',[]) subplot(rows,1,3) plot(H(1:n:end,1),H(1:n:end,8),'b',H(1:n:end,1),H(1:n:end,10),'r') axis([axismin axismax 0 1000]); text(textmin,800,'Incoming global radiation'); %legend('Incoming','Net','Location','southwest') ylabel('(W/m^2)'); set(gca,'XTick',[]) subplot(rows,1,4) %plot(M(1:n:end,1),M(1:n:end,29),'r.','MarkerSize',2) % hold on; plot(H(1:n:end,1),H(1:n:end,4),'b') axis([axismin axismax 0 20]); text(textmin,15,'Wind Speed'); ylabel('(m/s)'); set(gca,'XTick',[]) % hold off subplot(rows,1,5) plot(H(1:n:end,1),H(1:n:end,6),'b.') axis([axismin axismax 0 400]); text(textmin,300,'Wind direction'); ylabel('(^\circ)'); set(gca,'XTick',[]) subplot(rows,1,6) %bar([axismin axismax],[0 0],'b') bar(H(1:n:end,1),H(1:n:end,11),'b') axis([axismin axismax 0 5]); text(textmin,3,'Precipitation'); ylabel('(mm/15 min)'); xlabel (['Date (in ',num2str(Year),')']) set(gca,'XTick',TickNum,'XTickLabel',TickDate); %% Calculate monthly values of data including number of data %Find leap year if mod(Year, 400) == 0 LastFeb = [num2str(Year),'-02-29']; % leap year elseif mod(Year, 4) == 0 && mod(Year, 100) ~= 0 LastFeb = [num2str(Year),'-02-29']; % leap year else LastFeb = [num2str(Year),'-02-28']; % ordinary year end %Initialize summary Sum = zeros(12,16); % Temperature for k=1:12 switch k % select mont and criteria for extracting data from that month case 1 avgind = find(H(:,1) < datenum([num2str(Year),'-02-01'])); case 2 avgind = find(H(:,1) > datenum([num2str(Year),'-01-31']) & H(:,1) < datenum(LastFeb)); case 3 avgind = H(:,1) > datenum(LastFeb) & H(:,1) < datenum([num2str(Year),'-04-01']); case 4 avgind = H(:,1) > datenum([num2str(Year),'-03-31']) & H(:,1) < datenum([num2str(Year),'-05-01']); case 5 avgind = H(:,1) > datenum([num2str(Year),'-04-31']) & H(:,1) < datenum([num2str(Year),'-06-01']); case 6 avgind = H(:,1) > datenum([num2str(Year),'-05-31']) & H(:,1) < datenum([num2str(Year),'-07-01']); case 7 avgind = H(:,1) > datenum([num2str(Year),'-06-31']) & H(:,1) < datenum([num2str(Year),'-08-01']); case 8 avgind = H(:,1) > datenum([num2str(Year),'-07-31']) & H(:,1) < datenum([num2str(Year),'-09-01']); case 9 avgind = H(:,1) > datenum([num2str(Year),'-08-31']) & H(:,1) < datenum([num2str(Year),'-10-01']); case 10 avgind = H(:,1) > datenum([num2str(Year),'-09-30']) & H(:,1) < datenum([num2str(Year),'-11-01']); case 11 avgind = H(:,1) > datenum([num2str(Year),'-10-31']) & H(:,1) < datenum([num2str(Year),'-12-01']); case 12 avgind = H(:,1) > datenum([num2str(Year),'-11-30']) & H(:,1) < datenum([num2str(Year+1),'-01-01']); end temp = H(avgind,2); % monthly average temperature 1 Sum(k,1) = mean(temp(~isnan(temp))); Sum(k,2) = length(temp(~isnan(temp))); clear temp temp = H(avgind,3);% monthly average temperature 2 Sum(k,3) = mean(temp(~isnan(temp))); Sum(k,4) = length(temp(~isnan(temp))); clear temp temp = H(avgind,2); % monthly average temperature 3 Sum(k,5) = mean(temp(~isnan(temp))); Sum(k,6) = length(temp(~isnan(temp))); clear temp temp = H(avgind,2); % monthly totalized positive temperature 3 Sum(k,7) = sum(temp(temp(~isnan(temp))>0)); Sum(k,8) = length(temp(~isnan(temp(temp>0)))); clear temp temp = H(avgind,8); % monthly average radiation Sum(k,9) = mean(temp(~isnan(temp))); Sum(k,10) = length(temp(~isnan(temp))); clear temp temp = H(avgind,8); % monthly sum of radiation Sum(k,11) = sum(temp(temp(~isnan(temp))>0)); Sum(k,12) = length(temp(~isnan(temp(temp>0)))); clear temp temp = H(avgind,9); % monthly taverage net radiation Sum(k,13) = sum(temp(~isnan(temp))); Sum(k,14) = length(temp(~isnan(temp))); clear temp temp = H(avgind,11); % monthly totalized precipitation Sum(k,15) = sum(temp(~isnan(temp))); Sum(k,16) = length(temp(~isnan(temp))); clear temp temp = H(avgind,4); % monthly average wind speed Sum(k,17) = mean(temp(~isnan(temp))); Sum(k,18) = length(temp(~isnan(temp))); clear temp clear avgind end try fidSum = fopen(['TRS_met_summary_',num2str(Year),'.tex'], 'w'); catch ME fclose(fidSum); error('Error writing to output file'); end Units = ... ['(\si{\degreeCelsius}) '; '$n$ '; '(\si{\percent}) '; '$n$ '; '(\si{\degreeCelsius}) '; '$n$ '; '(\si{\degreeCelsius}) '; '$n$ '; '(\si{\watt\per\meter\squared}) '; '$n$ '; '(\si{\watt\per\meter\squared}) '; '$n$ '; '(\si{\watt\per\meter\squared}) '; '$n$ '; '(\si{\milli\meter}) '; '$n$ '; '(\si{\meter\per\second}) '; '$n$ ' ]; fprintf(fidSum,'\\begin{table}[ht]\n'); fprintf(fidSum,'\\caption{Monthly averages of meteorological parameters from the Tarfala Research Station automatic weather station \\Year.}\n'); fprintf(fidSum,'\\begin{center}{\\scriptsize\n'); fprintf(fidSum,'\\begin{tabular}{l d{0} d{0} d{0} d{0} d{0} d{0} d{0} d{0} d{0} d{0} d{0} d{0}}\n'); fprintf(fidSum,'\\toprule\n'); fprintf(fidSum,' &\\multicolumn{1}{c}{Jan.} &\\multicolumn{1}{c}{Feb.}&\\multicolumn{1}{c}{Mar.}&\\multicolumn{1}{c}{Apr.}&\\multicolumn{1}{c}{May }&\\multicolumn{1}{c}{Jun.}&\\multicolumn{1}{c}{Jul.}&\\multicolumn{1}{c}{Aug.}&\\multicolumn{1}{c}{Sep.}&\\multicolumn{1}{c}{Oct.}&\\multicolumn{1}{c}{Nov.}&\\multicolumn{1}{c}{Dec.}\\\\ \n'); l = 1;k = 1; %1 avg T fprintf(fidSum,'\\midrule\n \\multicolumn{13}{l}{Average air temperature} \\\\ \\\\ \n \\multicolumn{1}{r}{%s} &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f \\\\ \n \\multicolumn{1}{r}{%s} &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f \\\\', ... Units(k,1:37), Sum(1,k), Sum(2,k), Sum(3,k), Sum(4,k), Sum(5,k),... Sum(6,k), Sum(7,k), Sum(8,k), Sum(9,k), Sum(10,k),... Sum(11,k), Sum(12,k),... Units(k+1,1:37), Sum(1,k+1), Sum(2,k+1), Sum(3,k+1), Sum(4,k+1), Sum(5,k+1),... Sum(6,k+1), Sum(7,k+1), Sum(8,k+1), Sum(9,k+1), Sum(10,k+1),... Sum(11,k+1), Sum(12,k+1) ); l = l + 1; k = k + 2; %2 avg Rh fprintf(fidSum,'\\midrule\n \\multicolumn{13}{l}{Average relative humidity} \\\\ \\\\ \n \\multicolumn{1}{r}{%s} &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f \\\\ \n \\multicolumn{1}{r}{%s} &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f \\\\', ... Units(k,1:37), Sum(1,k), Sum(2,k), Sum(3,k), Sum(4,k), Sum(5,k),... Sum(6,k), Sum(7,k), Sum(8,k), Sum(9,k), Sum(10,k),... Sum(11,k), Sum(12,k),... Units(k+1,1:37), Sum(1,k+1), Sum(2,k+1), Sum(3,k+1), Sum(4,k+1), Sum(5,k+1),... Sum(6,k+1), Sum(7,k+1), Sum(8,k+1), Sum(9,k+1), Sum(10,k+1),... Sum(11,k+1), Sum(12,k+1) ); l = l + 1; k = k + 2; %3 avg T fprintf(fidSum,'\\midrule\n \\multicolumn{13}{l}{Average air temperature} \\\\ \\\\ \n \\multicolumn{1}{r}{%s} &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f \\\\ \n \\multicolumn{1}{r}{%s} &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f \\\\', ... Units(k,1:37), Sum(1,k), Sum(2,k), Sum(3,k), Sum(4,k), Sum(5,k),... Sum(6,k), Sum(7,k), Sum(8,k), Sum(9,k), Sum(10,k),... Sum(11,k), Sum(12,k),... Units(k+1,1:37), Sum(1,k+1), Sum(2,k+1), Sum(3,k+1), Sum(4,k+1), Sum(5,k+1),... Sum(6,k+1), Sum(7,k+1), Sum(8,k+1), Sum(9,k+1), Sum(10,k+1),... Sum(11,k+1), Sum(12,k+1) ); l = l + 1; k = k + 2; %4 sum pos T fprintf(fidSum,'\\midrule\n \\multicolumn{13}{l}{Positive degree sum} \\\\ \\\\ \n \\multicolumn{1}{r}{%s} &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f \\\\ \n \\multicolumn{1}{r}{%s} &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f \\\\', ... Units(k,1:37), Sum(1,k), Sum(2,k), Sum(3,k), Sum(4,k), Sum(5,k),... Sum(6,k), Sum(7,k), Sum(8,k), Sum(9,k), Sum(10,k),... Sum(11,k), Sum(12,k),... Units(k+1,1:37), Sum(1,k+1), Sum(2,k+1), Sum(3,k+1), Sum(4,k+1), Sum(5,k+1),... Sum(6,k+1), Sum(7,k+1), Sum(8,k+1), Sum(9,k+1), Sum(10,k+1),... Sum(11,k+1), Sum(12,k+1) ); l = l + 1; k = k + 2; %5 avg Ir fprintf(fidSum,'\\midrule\n \\multicolumn{13}{l}{Average incoming global radiation} \\\\ \\\\ \n \\multicolumn{1}{r}{%s} &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f \\\\ \n \\multicolumn{1}{r}{%s} &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f \\\\', ... Units(k,1:37), Sum(1,k), Sum(2,k), Sum(3,k), Sum(4,k), Sum(5,k),... Sum(6,k), Sum(7,k), Sum(8,k), Sum(9,k), Sum(10,k),... Sum(11,k), Sum(12,k),... Units(k+1,1:37), Sum(1,k+1), Sum(2,k+1), Sum(3,k+1), Sum(4,k+1), Sum(5,k+1),... Sum(6,k+1), Sum(7,k+1), Sum(8,k+1), Sum(9,k+1), Sum(10,k+1),... Sum(11,k+1), Sum(12,k+1) ); l = l + 1; k = k + 2; %6 Sum Ir fprintf(fidSum,'\\midrule\n \\multicolumn{13}{l}{Global incoming energy sum} \\\\ \\\\ \n \\multicolumn{1}{r}{%s} &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f \\\\ \n \\multicolumn{1}{r}{%s} &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f \\\\', ... Units(k,1:37), Sum(1,k), Sum(2,k), Sum(3,k), Sum(4,k), Sum(5,k),... Sum(6,k), Sum(7,k), Sum(8,k), Sum(9,k), Sum(10,k),... Sum(11,k), Sum(12,k),... Units(k+1,1:37), Sum(1,k+1), Sum(2,k+1), Sum(3,k+1), Sum(4,k+1), Sum(5,k+1),... Sum(6,k+1), Sum(7,k+1), Sum(8,k+1), Sum(9,k+1), Sum(10,k+1),... Sum(11,k+1), Sum(12,k+1) ); l = l + 1; k = k + 2; %7 Sum net fprintf(fidSum,'\\midrule\n \\multicolumn{13}{l}{Average net radiation} \\\\ \\\\ \n \\multicolumn{1}{r}{%s} &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f \\\\ \n \\multicolumn{1}{r}{%s} &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f \\\\', ... Units(k,1:37), Sum(1,k), Sum(2,k), Sum(3,k), Sum(4,k), Sum(5,k),... Sum(6,k), Sum(7,k), Sum(8,k), Sum(9,k), Sum(10,k),... Sum(11,k), Sum(12,k),... Units(k+1,1:37), Sum(1,k+1), Sum(2,k+1), Sum(3,k+1), Sum(4,k+1), Sum(5,k+1),... Sum(6,k+1), Sum(7,k+1), Sum(8,k+1), Sum(9,k+1), Sum(10,k+1),... Sum(11,k+1), Sum(12,k+1) ); l = l + 1; k = k + 2; %8 sum P fprintf(fidSum,'\\midrule\n \\multicolumn{13}{l}{Totalized precipitation} \\\\ \\\\ \n \\multicolumn{1}{r}{%s} &%.2f &%.2f &%.2f &%.2f &%.2f &%.2f &%.2f &%.2f &%.2f &%.2f &%.2f &%.2f \\\\ \n \\multicolumn{1}{r}{%s} &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f \\\\', ... Units(k,1:37), Sum(1,k), Sum(2,k), Sum(3,k), Sum(4,k), Sum(5,k),... Sum(6,k), Sum(7,k), Sum(8,k), Sum(9,k), Sum(10,k),... Sum(11,k), Sum(12,k),... Units(k+1,1:37), Sum(1,k+1), Sum(2,k+1), Sum(3,k+1), Sum(4,k+1), Sum(5,k+1),... Sum(6,k+1), Sum(7,k+1), Sum(8,k+1), Sum(9,k+1), Sum(10,k+1),... Sum(11,k+1), Sum(12,k+1) ); l = l + 1; k = k + 2; %9 avg wsp fprintf(fidSum,'\\midrule\n \\multicolumn{13}{l}{Average wind speed} \\\\ \\\\ \n \\multicolumn{1}{r}{%s} &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f &%.1f \\\\ \n \\multicolumn{1}{r}{%s} &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f &%.0f \\\\', ... Units(k,1:37), Sum(1,k), Sum(2,k), Sum(3,k), Sum(4,k), Sum(5,k),... Sum(6,k), Sum(7,k), Sum(8,k), Sum(9,k), Sum(10,k),... Sum(11,k), Sum(12,k),... Units(k+1,1:37), Sum(1,k+1), Sum(2,k+1), Sum(3,k+1), Sum(4,k+1), Sum(5,k+1),... Sum(6,k+1), Sum(7,k+1), Sum(8,k+1), Sum(9,k+1), Sum(10,k+1),... Sum(11,k+1), Sum(12,k+1) ); fprintf(fidSum,'\\bottomrule\n'); fprintf(fidSum,'\\end{tabular}}\n'); fprintf(fidSum,'\\end{center}\n'); fprintf(fidSum,'\\label{tab:Summary}\n'); fprintf(fidSum,'\\end{table}\n'); fclose(fidSum); % Pr(k) = nansum(P(avgind,10)); % monthly totalized precipitation % set(gca,'XTick',[]) %subplot(rows,1,7) %plot(M(1:n:end,1),M(1:n:end,35),'b') % axis([axismin axismax 0 5]); % text(textmin,3,'SR50A height'); % ylabel('(m)'); % set(gca,'XTick',TickNum,'XTickLabel',TickDate);