How to superimpose a plot of observations with a plot of information from a model?

2 views (last 30 days)
I was given a program and one block of code included the maximum wind speed of a hurricane. When you ran the program, a plot of the maximum wind speed with respect to time showed up (observation data). I am asked to tweak the code so that the information from the hurricane models also show up on that same graph (superimpose). I was given a file with numbers and was told which column was the maximum wind speed for the hurricane model, Katia. He told me this:
  • * * * * * *
" Please find GFDL model's Katia(2011090112) run. You can find track and intensity from these files. All these files have track and intensity data , but different formats. Choose one easier for you to plot. * * * * * * * * * *
To plot model trak and intensity, you mat want to use katia12l.2011090112.stats.complete file.
cp katia12l.2011090112.stats.complete track.dat load track.dat plot(track(:,1),track(:,4)) winds plot(track(:,1),track(:,6)) min pressure plot(track(:,2),track(:,3)) track
And then you can mege these commands with scipts that plot observation (intentrackobs.m)"
  • *
  • *
Intentrackobs.m is the following program:
clear all; opengl software set(0,'defaultaxesfontweight','bold') set(0,'defaultaxesfontsize',12) set(0,'defaultaxesfontName','times')
Name='KATRINA'; startdate = '2005082400' ; output1=['Inten.',Name,'.',startdate,'.png']; output2=['Track.',Name,'.',startdate,'.png'];
[tm0,pmin0,vmx0,lon0,lat0]=read_message(Name,startdate(1:4)); vmx0=1.944*vmx0; % converting to Kt because message file winds are in m/s julday0=date2day00(startdate);
date0=startdate; year=str2num(date0(1:4)); month0=str2num(date0(5:6)); day0=str2num(date0(7:8)); hour0=str2num(date0(9:10));
MM=str2mat('Apr','May','Jun','Jul','Aug','Sep','Oct','Nov'); MMF=str2mat('April','May','June','July','August','September','October','November'); tit2=['Tropical Cyclone ',upper(Name),'(',year,')']; tit1=['INITIAL TIME: ',sprintf('%4.4i',hour0),' UTC, ',int2str(day0),' ',MMF(month0-3,:),', ',year];
tm=tm0; lon=lon0; lat=lat0; mnp=pmin0; vmx=vmx0;
tmmin=tm(1); tmmax=tm(length(tm)); juldaym=mean([tmmin,tmmax]); datem=day2date00(year,juldaym); monthm=str2num(datem(5:6)); daym=str2num(datem(7:8)); hourm=str2num(datem(9:10));
juldaye=tmmax; datee=day2date00(year,juldaye) monthe=str2num(datee(5:6)); daye=str2num(datee(7:8)); houre=str2num(datee(9:10));
xtick=str2mat([MM(month0-3,:),'.',int2str(day0),',',int2str(hour0),'Z'],... [MM(monthm-3,:),'.',int2str(daym),',',int2str(hourm),'Z'],... [MM(monthe-3,:),'.',int2str(daye),',',int2str(houre),'Z']);
% for Pmin %-----------
%f1=figure(1); f1=figure ; set(f1,'units','inches','position',[1,0.25,8,10]); set(f1,'paperposition',[0.25,0.5,8,10]); whitebg(f1,'w');
dj=4; subplot(2,1,1) ll=length(tm); p0=plot(tm,mnp);hold on set(p0,'linewidth',2.5,'color','k','linestyle','-'); p0a=plot(tm(1:dj:ll),mnp(1:dj:ll),'.');hold on set(p0a,'color','k','markersize',40)
set(gca,'position',[0.2,0.5,0.6,0.37]) set(gca,'xlim',[julday0,juldaye]) set(gca,'xtick',[julday0,juldaym,juldaye]) set(gca,'xticklabels',xtick)
tit=title({tit2,tit1});
%title(['INITIAL TIME: ',sprintf('%4.4i',hour0),' UTC, ',... %int2str(day0),' ',MMF(month0-3,:),' ',int2str(year)]); %ti=get(gca,'title'); %ti1=text(1,1,['Tropical Cyclone ',upper(Name),'(',int2str(year),')']); %set(ti1,'units','normalized','position',[0.15,1.2]); %set(ti1,'fontName','times','fontweight','bold','fontsize',18) ylabel('Minimum sea-level pressure (hPa)') grid
% for Vmax %----------- subplot(2,1,2) ll=length(tm); p0=plot(tm,vmx,'-');hold on set(p0,'linewidth',2.5,'color','k'); p0a=plot(tm(1:dj:ll),vmx(1:dj:ll),'.');hold on set(p0a,'color','k','markersize',40)
%%te0=text(1,1,'OBSR'); %%set(te0,'units','normalized','position',[.30,0.35]) %%set(te0,'fontName','times','fontweight','bold','fontsize',13,'units','data') %%legx=get(te0,'position'); %%line([legx(1)-.1,legx(1)-.7],[legx(2),legx(2)],'linewidth',2.5,'color','k','linestyle','-') %%te1=text(1,1,run_name1); %%set(te1,'units','normalized','position',[.30,0.3]) %%set(te1,'fontName','times','fontweight','bold','fontsize',13,'units','data') %%legx=get(te1,'position'); %%line([legx(1)-.1,legx(1)-.7],[legx(2),legx(2)],'linewidth',2.5,'color','b','linestyle','-') %%te2=text(1,1,txt2); %%set(te2,'units','normalized','position',[0.30,0.25]) %%set(te2,'fontName','times','fontweight','bold','fontsize',13,'units','data') %%legx=get(te2,'position'); %%line([legx(1),legx(1)-1],[legx(2),legx(2)],'linewidth',2.5,'color','r','linestyle','-') set(gca,'position',[0.2,0.07,0.6,0.37]) set(gca,'xlim',[julday0,juldaye]) %set(gca,'ylim',[70,150]) set(gca,'xtick',[julday0,juldaym,juldaye]) set(gca,'xticklabels',xtick) ylabel(['Maximum Wind (kt)']) grid set(gca,'position',[0.2 0.07 0.6 0.37])
drawnow print ('-dpng',output1) % ----------------------------------------------------------------------------- % for track %-----------
%f2=figure(2); f2=figure; set(f2,'units','inches','position',[3,0.25,8,10]); set(f2,'paperposition',[0.25,0.5,8,10]); whitebg(f2,'w');
subplot(2,1,1) load gebco.dat, g=gebco; clear gebco; cl=plot(g(:,3),g(:,2)); set(cl,'linewidth',1,'color','k') hold on ll=length(lon); p0=plot(lon,lat); p0a=plot(lon(1:dj:ll),lat(1:dj:ll),'.');hold on set(p0,'linewidth',2.5,'color','k'),set(p0a,'color','k','markersize',30)
ll=length(lon); hx12=[lon(1),lon(ll)]; hy12=[lat(1),lat(ll)]; nd=length(lon); tt=[tm(1),tm(1)+(nd-1)/dj]; dhx=hx12(length(hx12))-hx12(1); dhy=hy12(length(hx12))-hy12(1); nor=sqrt(dhx^2+dhy^2); tx=hx12+5*abs(dhy)/nor; ty=hy12+5*abs(dhx)/nor; for n=1:length(tx) fprintf(1,'%6.2f',tt(n)) date1=day2date00(year,tt(n)); month1=str2num(date1(5:6)); day1=str2num(date1(7:8)); hour1=str2num(date1(9:10)); txtstr=[MM(month1-3,:),'.',int2str(day1),',',int2str(hour1),'Z']; ttmp=text(tx(n),ty(n),txtstr); set(ttmp,'fontName','times','fontweight','bold','fontsize',13) p3=plot([tx(n),hx12(n)],[ty(n),hy12(n)]) set(p3,'linewidth',1.5,'color','g'); end
%axis equal tmplat=[lat]; tmplon=[lon];
latmin=round(mean(tmplat))-30;latmax=round(mean(tmplat))+30; longmin=round(mean(tmplon))-35;longmax=round(mean(tmplon))+35; set(gca,'xlim',[longmin,longmax]); set(gca,'ylim',[latmin,latmax]); %set(gca,'xlim',[-95,-80]); %set(gca,'ylim',[20,35]);
ti1=text(1,1,['Tropical Cyclone ',upper(Name),' (',int2str(year),')']); set(ti1,'units','normalized','position',[0.15,1.15]); set(ti1,'fontName','times','fontweight','bold','fontsize',13,'fontsize',18) % labels ytlb=[-15:5:60]'; ytlb=ytlb(ytlb>latmin&ytlb<latmax); ytlb1=[num2str(ytlb),setstr(176*ones(length(ytlb),1)),setstr(78*ones(length(ytlb),1))]; %ytlb1=sprintf(['%3i',setstr(176),'N'],abs(ytlb)); ytlb1=(reshape(ytlb1,[5,length(ytlb1)/5]))'; set(gca,'ytick',ytlb,'yticklabels',ytlb1); ytlb=[-150:5:0]'; ytlb=ytlb(ytlb>longmin&ytlb<longmax); if length(ytlb)>6, ytlb=[-150:10:0]'; ytlb=ytlb(ytlb>longmin&ytlb<longmax); end %ytlb1=sprintf(['%3i',setstr(176),'W'],abs(ytlb)); ytlb1=(reshape(ytlb1,[5,length(ytlb1)/5]))'; ytlb1=[num2str(abs(ytlb)),setstr(176*ones(length(ytlb),1)),setstr(87*ones(length(ytlb),1))]; %ytlb1=sprintf(['%3i',setstr(176),'W'],abs(ytlb)); set(gca,'xtick',ytlb,'xticklabels',ytlb1); grid
drawnow print ('-dpng',output2)

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!