%BAOPLT300_loop % Plot data from GMD 300m BAO data loggers % Dave has split these data into daily files and corrected the PRT and Wind Dir %You are correct on all counts. yes, m/s on the wind speed. %The PRT temperature is (accurate to within 0.1 deg C) is: a0 + %a1*(D3_PRTRES_Avg) + a2*(D3_PRTRES_Avg)^2 + a3*(D3_PRTRES_Avg)^3, %a0 = -246.058 %a1 = 2.36626 %a2 = 0.000923807 %a3 = 2.43906e-007 %The sonic anemometer is on the south boom % The north pole points to 322 degrees. %The the true wind direction should then be (SE2_SONICWD_Avg) - 38 degrees % % 300m %"TOA5","CR1000","CR1000","2846","CR1000.Std.11","CPU:060516b_Met.CR1","10044","wkt30m" %"TIMESTAMP","RECORD","SE1_SONICWS_Avg","SE2_SONICWD_Avg","SE2_SONICWD_Std","SE1_SONICWS_Std", %"SE3_SONICTEMP_Avg","D3_PRTRES_Avg","D3_PRTTEMPAPPROX_Avg","SE4_HMPTEMP_Avg","SE9_HMPRH_Avg", %"P1_AIRFLOW_Avg","SE10_ACPOWER","BT_BATVOLT_Min","PRT_AVG","SONICWDT_AVG" %"TS","RN","ms","Deg","Deg","ms","mV","","Deg C","Deg C","%","Counts","ms","Volts","DegC","Deg True" %"2007-06-14 14:22:00",0,3.37,230.5,0.502,0.012,13.42,105.8471,14.24,15.09,42.15,0,110.6,13.31,15.0,192.5 %"2007-06-14 14:22:30",1,3.418,231.1,0.269,0.026,13.38,105.8718,14.3,15.11,42.26,0,110.6,13.31,15.1,193.1 function [retval] = BAOPLT300_loop_test( byd, eyd, yr ) %clear pltflg = 1; %flag to plot the data: 1=yes, 0=no pngflg = 1; %flag to make a png of the plot: 1=yes, 0=no %jds = input('Input Julian day to plot ','s'); %jd = str2num(jds); inpath = 'z:\bao\Tower\Processed\daily\'; outpath = 'z:\bao\Tower\Processed_Images\daily\'; cd(inpath); if( nargin == 0 ) bds = input('Input the beginning day to process (165) ','s'); eds = input('Input end day to process (168) ','s'); ym = input('Input year to process (2007) ','s'); year = str2double( ym ); bd = str2double(bds); ed = str2double(eds); else year = yr; bd = byd; ed = eyd; end nmax = 3000; year = zeros(1,nmax); month = zeros(1,nmax); day = zeros(1,nmax); hour = zeros(1,nmax); minute = zeros(1,nmax); sec = zeros(1,nmax); yday = zeros(1,nmax); bao300 = zeros( nmax, 15); % 76 days 69 = 233 for i = bd:ed name300 = sprintf('BAO_300_%4d%03d.dat', year, i ); fprintf(1,'*** Processing %s ***\n', name10 ); disp(['Reading ',m(ii).name]); f300 = fopen(m(ii).name,'r'); %load f:\BAO_CR10\data\bao10.dat %read in 4 header lines tmp=fgetl(f300); tmp=fgetl(f300); tmp=fgetl(f300); tmp=fgetl(f300); bao300 = zeros( 3000, 15); n = 0; while ~feof(f300) n = n + 1; temp=fgetl(f300); year(n)= str2double(temp(2:5)); month(n)= str2double(temp(7:8)); day(n)= str2double(temp(10:11)); hour(n)= str2double(temp(13:14)); minute(n)= str2double(temp(16:17)); sec(n)= str2double(temp(19:20)); lt= length(temp); % strip off date/time info from first part of each line s1 = temp(23:lt); % convert each , to a blank space s1 = strrep( s1, ',', ' ' ); [data]=sscanf(s1,'%f') ; bao300(n,1) = data(2,1); % WS ms Sonic bao300(n,2) = data(3,1); % WD deg Sonic bao300(n,3) = data(4,1); % WD std deg Sonic bao300(n,4) = data(5,1); % WS STD Sonic bao300(n,5) = data(6,1); % Temp deg C Sonic bao300(n,6) = data(7,1); % PRT temp resistance bao300(n,7) = data(8,1); % PRT temp deg C bao300(n,8) = data(9,1); % HMP Temp Deg C bao300(n,9) = data(10,1); % RH bao300(n,10) = data(11,1); % airflow bao300(n,11) = data(12,1); % VAC bao300(n,12) = data(13,1); % Bat VDC bao300(n,13) = data(14,1); % PRT corrected deg C bao300(n,14) = data(15,1); % Dir true end %while ~feof(f300) yday = datenum(year,month,day,hour,minute,sec)-datenum(year-1,12,31); bao300(1:n,15) = yday(1:n); % Year Day jdb = fix(bao300(1,15)); jde = fix(bao300(n,15)) + 1; if pltflg % do plots if pltflg=1 figure(1) clf subplot(2,1,1),plot(bao300(1:n,15),bao300(1:n,5),'k') hold %subplot(2,1,1),plot(bao300(1:n,15),bao300(1:n,7),'r') subplot(2,1,1),plot(bao300(1:n,15),bao300(1:n,8),'g') subplot(2,1,1),plot(bao300(1:n,15),bao300(1:n,13),'b') %xlabel('Year Day') ylabel('T (C)') stitle=sprintf('BAO %d 300m', year(1)); title(stitle) ylim([0,40]) %axis([110,116,-10,30]) legend('Sonic T','HMP','PRTc','Location','Best') ylim([min(bao300(1:n,13))-3,max(bao300(1:n,13))+3]) xlim([jdb,jde]) subplot(2,1,2),plot(bao300(1:n,15),bao300(1:n,9)) %hold %subplot(2,1,2),plot(data300(1:n,3),data300(1:n,6),'r') xlabel('Year Day') ylabel('Rh (%)') %axis([110,116,0,100]) %legend('300m') %legend('10m','50m') ylim([0,109]) xlim([jdb,jde]) if pngflg print_buffer = [outpath, m(ii).name(1:15),'_TRH.png']; print('-dpng ', print_buffer); end figure(2) clf subplot(2,1,1),plot(bao300(1:n,15),bao300(1:n,1)) %hold %subplot(2,1,1),plot(data300(1:n,3),data300(1:n,7),'r') %xlabel('Year Day') ylabel('Spd (m/s)') title(stitle) %axis([110,116,0,30]) ylim([0,max(bao300(1:n,1))+5]) xlim([jdb,jde]) %legend('300m') %legend('10m','50m') subplot(2,1,2),plot(bao300(1:n,15),bao300(1:n,14)) %hold %subplot(2,1,2),plot(data300(1:n,3),data300(1:n,8),'r') xlabel('Year Day') ylabel('Dir (true)') xlim([jdb,jde]) %axis([110,116,0,360]) %legend('300m') %legend('10m','50m') xlim([jdb,jde]) if pngflg print_buffer = [outpath, m(ii).name(1:15),'_WS.png']; print('-dpng ', print_buffer); end figure(3) clf plot(bao300(1:n,15),bao300(1:n,5),'k') hold plot(bao300(1:n,15),bao300(1:n,7),'r') plot(bao300(1:n,15),bao300(1:n,8),'b') plot(bao300(1:n,15),bao300(1:n,13),'g') xlabel('Year Day (UTC)') ylabel('T (C)') title(stitle) ylim([min(bao300(1:n,13))-3,max(bao300(1:n,13))+3]) %axis([110,116,-10,30]) legend('Sonic T','PRT est','HMP T','PRT corrected','Location','Best') xlim([jdb,jde]) if pngflg print_buffer = [outpath, m(ii).name(1:15),'_Tall.png']; print('-dpng ', print_buffer); end % Scatter plot of HMP vs corrected PRT at 300m figure(4) clf plot(bao300(1:n,8),bao300(1:n,13),'k.') xlabel('HMP (C)') ylabel('PRT (C)') title(stitle) axis('square') %ylim([min(bao300(1:n,15))-3,max(bao300(1:n,15))+3]) %axis([110,116,-10,30]) %xlim([jdb,jde]) if pngflg print_buffer = [outpath, m(ii).name(1:15),'_HMPPRT.png']; print('-dpng ', print_buffer); end end % if doing plots fclose(f300); end % loop though each day