diff --git a/main_test.m b/main_test.m index 2cedeb443706fe31860c646214e4955805838681..741f00f0794694bedeb3e86fa4788cc17ba20bdf 100644 --- a/main_test.m +++ b/main_test.m @@ -36,9 +36,9 @@ zf0=csvread('H_ble.csv'); LAI0=csvread('LAI_ble.csv'); - LAIdate=LAI0(:,1); % dates de mesure de l'indice foliaire - gLAIval=LAI0(:,2); % valeur de l'indice foliaire à ces dates - LAIval=LAI0(:,3); % valeur de l'indice foliaire total à ces dates + LAIdate=LAI0(:,1); % Dates + gLAIval=LAI0(:,2); % Grean LAI at those dates + LAIval=LAI0(:,3); % Total LAI at those dates glai=interp1(LAIdate,gLAIval,doy,'pchip','extrap'); lai=interp1(LAIdate,LAIval,doy,'pchip','extrap'); @@ -94,7 +94,9 @@ betastress=0.01; end -% bounding / stress +% bounding (can be adusted according to pot and stress baselines) + +% bounding / stress rns(les<lesstress)=rnsstress(les<lesstress); g(les<lesstress)=gstress(les<lesstress); @@ -184,36 +186,23 @@ hpot2=hspot2+hvpot2; a=doy-floor(doy); time=find(a>0.485-(1/24) & a<0.515-(1/24)); time2=find(a>0.485-(1/24) & a<0.515-(1/24) & leobs >-100 & hobs>-100 & rnobs >-100 & gobs> - 100); -time3=time2; -RMSE_LE = sqrt(mean((leobs(time2) - le(time2)).^2)) % RMSE -RMSE_LE2 = sqrt(mean((leobs(time2) - le2(time2)).^2)) % RMSE 2 +% Rott Mean Square Error on Latent Heat LE, Sensible HEat H, Net radiation +% Rn and Soil Heat Flux G -RMSE_H = sqrt(mean((hobs(time2) - h(time2)).^2)) % RMSE -RMSE_H2 = sqrt(mean((hobs(time2) - h2(time2)).^2)) % RMSE 2 +RMSE_LE = sqrt(mean((leobs(time2) - le(time2)).^2)) % RMSE Layer +RMSE_LE2 = sqrt(mean((leobs(time2) - le2(time2)).^2)) % RMSE Patch -RMSE_RN = sqrt(mean((rnobs(time2) - rn(time2)).^2)) % RMSE -RMSE_RN2 = sqrt(mean((rnobs(time2) - rn2(time2)).^2)) % RMSE 2 +RMSE_H = sqrt(mean((hobs(time2) - h(time2)).^2)) % RMSE Layer +RMSE_H2 = sqrt(mean((hobs(time2) - h2(time2)).^2)) % RMSE Patch -RMSE_G = sqrt(mean((gobs(time2) - g(time2)).^2)) % RMSE -RMSE_G2 = sqrt(mean((gobs(time2) - g2(time2)).^2)) % RMSE 2 +RMSE_RN = sqrt(mean((rnobs(time2) - rn(time2)).^2)) % RMSE Layer +RMSE_RN2 = sqrt(mean((rnobs(time2) - rn2(time2)).^2)) % RMSE Patch -figure(1) -plot(doy(time),ts(time),'bo',doy(time),ts2(time),'bs',doy(time),tspot(time),'go',doy(time),tsstress(time),'ro',doy(time),tspot2(time),'gs',doy(time),tsstress2(time),'rs',doy(time),ta(time)-273.15,':k') -legend('simulé série','simulé parallèle','potentiel série','stress série','potentiel parallèle','stress parallèle','air') -title('température sol') +RMSE_G = sqrt(mean((gobs(time2) - g(time2)).^2)) % RMSE Layer +RMSE_G2 = sqrt(mean((gobs(time2) - g2(time2)).^2)) % RMSE Patch -figure(2) -plot(doy(time),tv(time),'bo',doy(time),tv2(time),'bs',doy(time),tvpot(time),'go',doy(time),tvstress(time),'ro',doy(time),tvpot2(time),'gs',doy(time),tvstress2(time),'rs',doy(time),ta(time)-273.15,':k') -legend('simulé série','simulé parallèle','potentiel série','stress série','potentiel parallèle','stress parallèle','air') -title('température végétation') - -figure(3) -plot(doy(time),trad(time),'bo',doy(time),trad2(time),'bs',doy(time),tradpot(time),'go',doy(time),tradstress(time),'ro',doy(time),tradpot2(time),'gs',doy(time),tradstress2(time),'rs',doy(time),ta(time)-273.15,':k',doy(time),tsobs(time),'ok') -legend('simulé série','simulé parallèle','potentiel série','stress série','potentiel parallèle','stress parallèle','air') -title('température radiative') - -figure(4) +figure(1) plot(leobs(time2),le(time2),'bo',leobs(time2),le2(time2),'rs') hold on lsline @@ -221,128 +210,51 @@ axis([-100 500 -100 500]); legend('Series model','Parallel model') xlabel('Observed Latent Heat Flux at midday [W/m2]') ylabel('Retrieved Latent Heat Flux at midday [W/m2]') - -% figure(12) -% %plot(doy(time),lev(time)./levpot(time),'ko',doy(time),lev2(time)./levpot(time),'bo',... -% % doy(time),les(time)./lespot(time),'k*',doy(time),les2(time)./lespot(time),'b*') -% % plot(doy(time),lev(time)./levpot(time),'go',doy(time),lev(time)./levstress(time),'gs',... -% % doy(time),les(time)./lespot(time),'ro',doy(time),les(time)./lesstress(time),'rs') -% plot(doy(time),Hv(time)./Hvpot(time),'go',doy(time),Hv(time)./Hvstress(time),'gs',... -% doy(time),Hs(time)./Hspot(time),'ro',doy(time),Hs(time)./Hsstress(time),'rs') -% axis([0 180 -1 2]); -figure(5) -plot(doy(time2),leobs(time2)./lepot(time2),'rs',doy(time2),le(time2)./lepot(time2),'bo',doy(time2),le2(time2)./lepot2(time2),'g*',doy(time2),lai(time2)/2,'g-') -%plot(doy(time2),leobs(time2).*2./(lepot(time2)+lepot2(time2)),'rs',doy(time2),lereal(time2)./lepot(time2),'ro',doy(time2),lereal2(time2)./lepot2(time2),'r*',doy(time2),le(time2)./lepot(time2),'bo',doy(time2),le2(time2)./lepot2(time2),'g*') -%axis([0 180 -0.2 1.2]); -legend('observed','series','parallel','LAI/2') -xlabel('Date [DOY]') -ylabel('Evaporative Efficiency LE/LEpot [-], LAI/2 [-]') - - figure(6) - %plot(1-leobs(time3)./lepot(time3),1-le(time3)./lepot(time3),'bo',1-leobs(time3)./lepot2(time3),1-le2(time3)./lepot2(time3),'rs') - - scatter(1-leobs(time3)./lepot(time3),1-le(time3)./lepot(time3),lepot(time3)/2,'bo') - - hold on - scatter(1-leobs(time3)./lepot2(time3),1-le2(time3)./lepot2(time3),lepot2(time3)/2,'rs') - %lsline - axis([-0.2 1.2 -0.2 1.2]); - legend('Series model','Parallel model') - xlabel('Observed surface stress at midday (-)') - ylabel('Simulated surface stress at midday (-)') +figure(2) +scatter(1-leobs(time2)./lepot(time2),1-le(time2)./lepot(time2),'bo') +hold on +scatter(1-leobs(time2)./lepot2(time2),1-le2(time2)./lepot2(time2),'rs') +axis([-0.2 1.2 -0.2 1.2]); +legend('Series model','Parallel model') +xlabel('Observed surface stress at midday (-)') +ylabel('Simulated surface stress at midday (-)') -figure(7) +figure(3) plot(doy(time2),leobs(time2),'rs',doy(time2),le(time2),'bo',doy(time2),le2(time2),'g*') axis([0 180 0 500]); legend('observé','série','parallèle') xlabel('Date [DOY]') ylabel('LE [W/m2]') -save('bounded.mat','doy','time2','le','le2') - - weight=0.2*ones(1,5); -t=doy(time3); -es=les(time3); -es2=les2(time3); -esp=lespot(time3); -esp2=lespot2(time3); +t=doy(time2); +es=les(time2); +es2=les2(time2); +esp=lespot(time2); +esp2=lespot2(time2); esf=filter(weight,1,es); es2f=filter(weight,1,es2); espf=filter(weight,1,esp); esp2f=filter(weight,1,esp2); - -figure(8) -%plot(doy(time2),les(time2),'bo',doy(time2),les2(time2),'g*',dateobs+0.5,pluie*10,'k-',doy(time2),lesobs(time2),'r:',doy(time2),leobs(time2),'rs') +figure(4) plot(t,esf,'bo:',t,es2f,'rs:',t,espf,'g:',t,esp2f,'g-') axis([0 180 0 250]); legend('Series model','Parallel model','Potential serie','potential parallel') xlabel('Date [DOY]') ylabel('Soil Evaporation, 5 days running average [W/m2]') -figure(12) +figure(5) plot(t,esf./espf,'bo:',t,es2f./esp2f,'rs:') axis([0 180 -0.1 1.1]); legend('Series model','Parallel model') xlabel('Date [DOY]') ylabel('Soil Evaporation Efficiency, 5 days running average [-]') -% figure(9) -% plot(doy(time2),lev2(time2),'bo',doy(time2),lev(time2),'b-',doy(time2),levpot(time2),'g:') -% legend('Old Srface Energy Balance Model TSEB','New Surface Energy Balance Model TSEB2','potentiel') -% xlabel('Date [DOY]') -% ylabel('LE [W/m2]') - -figure(9) - plot(doy(time2),(h(time2)-hpot(time2))./(hstress(time2)-hpot(time2)),'b-',doy(time2),(hobs(time2)-hpot(time2))./(hstress(time2)-hpot(time2)),'bo',... - doy(time2),1-le(time2)./lepot(time2),'r-',doy(time2),1-leobs(time2)./lepot(time2),'rs') - %axis([0 180 -0.2 1.2]); - legend('Stress H simulé','Stress H observé','Stress LE simulé','Stress LE observé') - xlabel('Date [DOY]') - ylabel('Stress [-]') - -figure(10) -plot(doy(time2),hv2(time2),'b-',doy(time2),hvpot2(time2),'r-',doy(time2),hvstress2(time2),'g-') -axis([0 180 -200 500]); - -figure(11) -plot(doy(time2),hv(time2),'b-',doy(time2),hvpot(time2),'r-',doy(time2),hvstress(time2),'g-') -axis([0 180 -200 500]); - -figure(13) - plot((hobs(time2)-hpot(time2))./(hstress(time2)-hpot(time2)),(h(time2)-hpot(time2))./(hstress(time2)-hpot(time2)),'bo',... - (hobs(time2)-hpot(time2))./(hstress(time2)-hpot(time2)),(h2(time2)-hpot(time2))./(hstress(time2)-hpot(time2)),'rs') - hold on - lsline - axis([-0.2 1.2 -0.2 1.2]); - legend('Series model','Parallel model') - xlabel('Observed surface stress at midday (-)') - ylabel('Simulated surface stress at midday (-)') - - - -% figure(13) -% plot(doy(time2),leobs(time2),'rs',doy(time2),lev(time2),'bo',doy(time2),lev2(time2),'g*',dateobs+0.5,pluie*10,'k:') -% axis([0 180 -200 500]); -% legend('observé total','végétation série','végétation parallèle','pluie') -% xlabel('Date [DOY]') -% ylabel('LE [W/m2], pluie (*0.1 mm)') - -% figure(14) -% plot(lesobs(time2),les(time2),'bo',lesobs(time2),les2(time2),'g*') -% hold on -% lsline -% axis([-100 400 -100 400]); -% title('évaporation du sol') -% legend('série','parallèle') -% xlabel('Flux observé à midi (W/m2)') -% ylabel('Flux simulé à midi (W/m2)') - -figure(15) +figure(6) plot(hobs(time2),h(time2),'bo',hobs(time2),h2(time2),'g*') hold on lsline @@ -351,7 +263,7 @@ legend('series','parallel') xlabel('Observed Sensible Heat Flux at Midday [W/m2]') ylabel('Simulated Sensible Heat Flux at Midday [W/m2]') -figure(16) +figure(7) plot(rnobs(time2),rn(time2),'bo',rnobs(time2),rn2(time2),'g*') hold on lsline @@ -359,7 +271,8 @@ axis([-100 700 -100 700]); legend('series','parallel') xlabel('Observed Net Radiation at Midday [W/m2]') ylabel('Simulated Net Radition at Midday [W/m2]') -figure(18) + +figure(8) plot(gobs(time2),g(time2),'bo',gobs(time2),g2(time2),'g*') hold on lsline @@ -368,33 +281,5 @@ legend('series','parallel') xlabel('Observed Soil Heat Flux at Midday [W/m2]') ylabel('Observed Soil Heat Flux at Midday [W/m2]') -figure(21) -plot(doy(time3),1-leobs(time3)./lepot(time3),'rs',doy(time3),(tsobs(time3)-tradpot(time3))./(tradstress(time3)-tradpot(time3)),'b-',... - doy(time3),(tsobs(time3)-tradpot2(time3))./(tradstress2(time3)-tradpot2(time3)),'g-') - axis([0 180 -0.2 1.2]); - legend('Stress observé','Stress tso-tsp/ts0-tsp serie','Stress tso-tsp/ts0-tsp parallel') - xlabel('Date [DOY]') - ylabel('Stress [-]') - -figure(22) -plot(1-leobs(time2)./lepot(time2),(tsobs(time2)-tradpot(time2))./(tradstress(time2)-tradpot(time2)),'bo',... - 1-leobs(time2)./lepot(time2),1-le(time2)./lepot(time2),'go') - hold on - lsline - axis([-0.2 1.2 -0.2 1.2]); - legend('Stress tso-tsp/ts0-tsp serie','Stress 1-LE/LEpot') - xlabel('Stress observé [-]') - ylabel('Stress simulé [-]') - - - figure(17) - plot(1-leobs(time2)./lepot(time2),(hobs(time2)-hpot(time2))./(hstress(time2)-hpot(time2)),'rs') - hold on - lsline - axis([-0.2 1.2 -0.2 1.2]); - legend('Stress LE','Stress H') - xlabel('Stress [-]') - ylabel('Stress [-]') - \ No newline at end of file