%tree girdle experiment clear load('girdle_data1.mat'); load('girdle_data2.mat'); allleafcodes = [Cam01conleafcodes' Cam01grdleafcodes' Cam02grdleafcodes' Cam02leafcodes' Cam03conspectracodes' Cam03grdspectracodes' Cam04conleafcodes' Cam04grdleafcodes' Cam05conleafcodescorrected' Cam05grdleafcodescorrected']; allspectrap = [Cam01conspectra' Cam01grdspectra' Cam02grdspectra' Cam02spectra' Cam03conspectra' Cam03grdspectra' Cam04conspectra' Cam04grdspectra' Cam05conspectracorrected' Cam05grdspectracorrected']; % allspectra = allspectrap(1:621,i) allspectrap(622:1411,i)+0.0034 allspectrap(1412,i)+0.006 allspectrap(1413:end,i)+0.0088 %correct spectra for minor offsets at 621 and 1411 bands for i = 1:length(allleafcodes) allspectra2(:,i) = allspectrap(:,i); ch1 = allspectrap(621,i)-allspectrap(622,i); ch2 = allspectrap(1411,i)-allspectrap(1413,i); ch3 = ch2/2; allspectra(:,i) = [allspectrap(1:621,i); allspectrap(622:1411,i)+ch1; allspectrap(1412,i)+ch3+ch1; allspectrap(1413:end,i)+ch2+ch1]; end for i = 1:length(allleafcodes) str1 = find(strcmp(allleafcodes(i),SpectraLeafcodes.Spectra_name_asgiven)==1); branchz(i) = SpectraLeafcodes.Code_Branch(str1); treez(i) = SpectraLeafcodes.Tree(str1); str3 = find(treez(i)==cencusgirdle.TagNumber); if ~isempty(str3) deathz(i) = cencusgirdle.GirdlingDeathDate(str3); treespz(i) = cencusgirdle.GenusSpecies(str3); wooddenz(i) = cencusgirdle.WoodDensity(str3); DBHz(i) = cencusgirdle.D_POM_cm(str3); heightz(i) = cencusgirdle.Height_m(str3); end Campz(i) = SpectraLeafcodes.Campaign(str1); %Campaign 1=21-Jan-16, 2=10-Feb-16, 3=01-Mar-16, 4=29-Mar-16, 5 %08-Jun-16, 6 = 06-Aug-16, 7 = 05-Oct-16, 8 = 08-Feb-17 treatz(i) = SpectraLeafcodes.Treatment(str1); leafnumz = string(SpectraLeafcodes.Code_Leaf(str1)); str2 = find(strcmp(leafnumz,leaftraitdata.Code_Leaf)==1); if ~isempty(str2) lph(i) = leaftraitdata.Photo_m(str2); lres(i) = leaftraitdata.Resp_m(str2); lLMA(i) = leaftraitdata.LMA(str2); lwatercon(i) = (leaftraitdata.Leaf_Fresh(str2)-leaftraitdata.Leaf_Dry(str2))./leaftraitdata.Leaf_Fresh(str2); ldate(i) = datenum(leaftraitdata.Sample_Date(str2)); end end strtrc = find(strcmp(string(treatz), 'CON')==1); strtrg = find(strcmp(string(treatz), 'GRD')==1); allspectragrd = allspectra(:,strtrg); deathzg = datenum(deathz(strtrg)); deathzg2 =deathzg-min(deathzg); ldateg = ldate(strtrg); treespzg = treespz(strtrg); Campzg = Campz(strtrg); treezg = treez(strtrg); diffdeathg = deathzg-ldateg; diffdeathg(diffdeathg>1000)=NaN; diffdeathg(diffdeathg<0)=NaN; wooddenzg = wooddenz(strtrg); DBHzg = DBHz(strtrg); heightzg = heightz(strtrg); lphg = lph(strtrg); lresg = lres(strtrg); lLMAg= lLMA(strtrg); lwatercong= lwatercon(strtrg); ldateg = ldate(strtrg); allspectracon = allspectra(:,strtrc); deathzc = datenum(deathz(strtrc)); treespzc = treespz(strtrc); Campzc = Campz(strtrc); treezc = treez(strtrc); wooddenzc = wooddenz(strtrc); DBHzc = DBHz(strtrc); heightzc = heightz(strtrc); lphc = lph(strtrc); lresc = lres(strtrc); lLMAc= lLMA(strtrc); lwaterconc= lwatercon(strtrc); ldatec = ldate(strtrc); %filter data sf1 = find(allspectracon(500,:)>.35 & allspectracon(500,:)<.7); sf2 = find(allspectragrd(500,:)>.37 & allspectragrd(650,:)<.62); % subplot(211) % plot(allspectragrd(:,sf2)) % subplot(212) % plot(allspectracon(:,sf1)) allspectragrd2 = allspectragrd(:,sf2); deathzg2c = deathzg2(sf2); diffdeathg2 = diffdeathg(sf2); id = find(diffdeathg2<50); nid = find(diffdeathg2>50); treespzg2 = treespzg(sf2); Campzg2 = string(Campzg(sf2)); treezg2 = treezg(sf2); wooddenzg2 = wooddenzg(sf2); DBHzg2 = DBHzg(sf2); heightzg2 = heightzg(sf2); lphg2 = lphg(sf2); lresg2 = lresg(sf2); lLMAg2= lLMAg(sf2); lwatercong2= lwatercong(sf2); ldateg2 = ldateg(sf2); allspectracon2 = allspectracon(:,sf1); deathzc2 = deathzc(sf1); treespzc2 = treespzc(sf1); Campzc2 = string(Campzc(sf1)); treezc2 = treezc(sf1); wooddenzc2 = wooddenzc(sf1); DBHzc2 = DBHzc(sf1); heightzc2 = heightzc(sf1); lphc2 = lphc(sf1); lresc2 = lresc(sf1); lLMAc2= lLMAc(sf1); lwaterconc2= lwaterconc(sf1); ldatec2 = ldatec(sf1); spec=350:2500; %plot campaign spectra differences figure(1) subplot(221) plot(spec,nanmean(allspectracon2(:,(find(strcmp('C1',Campzc2)==1)))')) hold on plot(spec,nanmean(allspectracon2(:,(find(strcmp('C2',Campzc2)==1)))'),'r') plot(spec,nanmean(allspectracon2(:,(find(strcmp('C3',Campzc2)==1)))'),'k') plot(spec,nanmean(allspectracon2(:,(find(strcmp('C4',Campzc2)==1)))'),'g') plot(spec,nanmean(allspectracon2(:,(find(strcmp('C5',Campzc2)==1)))'),'y') xlabel('Wavelength (nm)') ylabel('reflectance') xlim([400 2500]) title('Drought') %gtext('A') subplot(222) plot(spec,nanmean(allspectragrd2(:,(find(strcmp('C1',Campzg2)==1)))')) hold on plot(spec,nanmean(allspectragrd2(:,(find(strcmp('C2',Campzg2)==1)))'),'r') plot(spec,nanmean(allspectragrd2(:,(find(strcmp('C3',Campzg2)==1)))'),'k') plot(spec,nanmean(allspectragrd2(:,(find(strcmp('C4',Campzg2)==1)))'),'g') plot(spec,nanmean(allspectragrd2(:,(find(strcmp('C5',Campzg2)==1)))'),'y') legend('C1', 'C2','C3','C4','C5') xlabel('Wavelength (nm)') ylabel('reflectance') xlim([400 2500]) title('Girdle') %gtext('B') %difference subplot(223) plot(spec,(nanmean(allspectracon2(:,(find(strcmp('C1',Campzc2)==1)))'))-nanmean(allspectracon2(:,(find(strcmp('C2',Campzc2)==1)))'),'r') hold on plot(spec,(nanmean(allspectracon2(:,(find(strcmp('C1',Campzc2)==1)))'))-nanmean(allspectracon2(:,(find(strcmp('C3',Campzc2)==1)))'),'k') plot(spec,(nanmean(allspectracon2(:,(find(strcmp('C1',Campzc2)==1)))'))-nanmean(allspectracon2(:,(find(strcmp('C4',Campzc2)==1)))'),'g') plot(spec,(nanmean(allspectracon2(:,(find(strcmp('C1',Campzc2)==1)))'))-nanmean(allspectracon2(:,(find(strcmp('C5',Campzc2)==1)))'),'y') xlabel('Wavelength (nm)') ylabel('reflectance') xlim([400 2500]) title('Drought difference') %gtext('C') subplot(224) plot(spec,(nanmean(allspectragrd2(:,(find(strcmp('C1',Campzg2)==1)))'))-nanmean(allspectragrd2(:,(find(strcmp('C2',Campzg2)==1)))'),'r') hold on plot(spec,(nanmean(allspectragrd2(:,(find(strcmp('C1',Campzg2)==1)))'))-nanmean(allspectragrd2(:,(find(strcmp('C3',Campzg2)==1)))'),'k') plot(spec,(nanmean(allspectragrd2(:,(find(strcmp('C1',Campzg2)==1)))'))-nanmean(allspectragrd2(:,(find(strcmp('C4',Campzg2)==1)))'),'g') plot(spec,(nanmean(allspectragrd2(:,(find(strcmp('C1',Campzg2)==1)))'))-nanmean(allspectragrd2(:,(find(strcmp('C5',Campzg2)==1)))'),'y') xlabel('Wavelength (nm)') ylabel('reflectance') legend('C1-C2','C1-C3','C1-C4','C1-C5') xlim([400 2500]) title('Girdle difference') %gtext('D') %calculate NDVI spec1 = nanmean(allspectracon2(:,(find(strcmp('C1',Campzc2)==1)))'); spec2 = nanmean(allspectracon2(:,(find(strcmp('C2',Campzc2)==1)))'); spec3 = nanmean(allspectracon2(:,(find(strcmp('C3',Campzc2)==1)))'); spec4 = nanmean(allspectracon2(:,(find(strcmp('C4',Campzc2)==1)))'); spec5 = nanmean(allspectracon2(:,(find(strcmp('C5',Campzc2)==1)))'); spec1g = nanmean(allspectragrd2(:,(find(strcmp('C1',Campzg2)==1)))'); spec2g = nanmean(allspectragrd2(:,(find(strcmp('C2',Campzg2)==1)))'); spec3g = nanmean(allspectragrd2(:,(find(strcmp('C3',Campzg2)==1)))'); spec4g = nanmean(allspectragrd2(:,(find(strcmp('C4',Campzg2)==1)))'); spec5g = nanmean(allspectragrd2(:,(find(strcmp('C5',Campzg2)==1)))'); spec1s = nanstd(allspectracon2(:,(find(strcmp('C1',Campzc2)==1)))'); spec2s = nanstd(allspectracon2(:,(find(strcmp('C2',Campzc2)==1)))'); spec3s = nanstd(allspectracon2(:,(find(strcmp('C3',Campzc2)==1)))'); spec4s = nanstd(allspectracon2(:,(find(strcmp('C4',Campzc2)==1)))'); spec5s = nanstd(allspectracon2(:,(find(strcmp('C5',Campzc2)==1)))'); spec1gs = nanstd(allspectragrd2(:,(find(strcmp('C1',Campzg2)==1)))'); spec2gs = nanstd(allspectragrd2(:,(find(strcmp('C2',Campzg2)==1)))'); spec3gs = nanstd(allspectragrd2(:,(find(strcmp('C3',Campzg2)==1)))'); spec4gs = nanstd(allspectragrd2(:,(find(strcmp('C4',Campzg2)==1)))'); spec5gs = nanstd(allspectragrd2(:,(find(strcmp('C5',Campzg2)==1)))'); ndvi1 = (spec1(650)-spec1(305))./(spec1(650)+spec1(305)); ndvi2 = (spec2(650)-spec2(305))./(spec2(650)+spec2(305)); ndvi3 = (spec3(650)-spec3(305))./(spec3(650)+spec3(305)); ndvi4 = (spec4(650)-spec4(305))./(spec4(650)+spec4(305)); ndvi5 = (spec5(650)-spec5(305))./(spec5(650)+spec5(305)); ndvig1 = (spec1g(650)-spec1g(305))./(spec1g(650)+spec1g(305)); ndvig2 = (spec2g(650)-spec2g(305))./(spec2g(650)+spec2g(305)); ndvig3 = (spec3g(650)-spec3g(305))./(spec3g(650)+spec3g(305)); ndvig4 = (spec4g(650)-spec4g(305))./(spec4g(650)+spec4g(305)); ndvig5 = (spec5g(650)-spec5g(305))./(spec5g(650)+spec5g(305)); ndvi1s = (spec1s(650)-spec1s(305))./(spec1s(650)+spec1s(305)); ndvi2s = (spec2s(650)-spec2s(305))./(spec2s(650)+spec2s(305)); ndvi3s = (spec3s(650)-spec3s(305))./(spec3s(650)+spec3s(305)); ndvi4s = (spec4s(650)-spec4s(305))./(spec4s(650)+spec4s(305)); ndvi5s = (spec5s(650)-spec5s(305))./(spec5s(650)+spec5s(305)); ndvig1s = (spec1gs(650)-spec1gs(305))./(spec1gs(650)+spec1gs(305)); ndvig2s = (spec2gs(650)-spec2gs(305))./(spec2gs(650)+spec2gs(305)); ndvig3s = (spec3gs(650)-spec3gs(305))./(spec3gs(650)+spec3gs(305)); ndvig4s = (spec4gs(650)-spec4gs(305))./(spec4gs(650)+spec4gs(305)); ndvig5s = (spec5gs(650)-spec5gs(305))./(spec5gs(650)+spec5gs(305)); title('>0.5 WD Drought') %plot spectral changes by functional type wdcg = find(wooddenzc2>0.5); wdcl = find(wooddenzc2<0.5); allspectracon2wdg=allspectracon2(:,wdcg); allspectracon2wdl=allspectracon2(:,wdcl); Campzc2wdg = Campzc2(wdcg); Campzc2wdl = Campzc2(wdcl); wdgg = find(wooddenzg2>0.5); wdgl = find(wooddenzg2<0.5); allspectragrd2wdg=allspectragrd2(:,wdgg); allspectragrd2wdl=allspectragrd2(:,wdgl); Campzg2wdg = Campzg2(wdgg); Campzg2wdl = Campzg2(wdgl); figure(2) subplot(321) plot(spec,(nanmean(allspectracon2wdg(:,(find(strcmp('C1',Campzc2wdg)==1)))'))-nanmean(allspectracon2wdg(:,(find(strcmp('C2',Campzc2wdg)==1)))'),'r') hold on plot(spec,(nanmean(allspectracon2wdg(:,(find(strcmp('C1',Campzc2wdg)==1)))'))-nanmean(allspectracon2wdg(:,(find(strcmp('C3',Campzc2wdg)==1)))'),'k') plot(spec,(nanmean(allspectracon2wdg(:,(find(strcmp('C1',Campzc2wdg)==1)))'))-nanmean(allspectracon2wdg(:,(find(strcmp('C4',Campzc2wdg)==1)))'),'g') plot(spec,(nanmean(allspectracon2wdg(:,(find(strcmp('C1',Campzc2wdg)==1)))'))-nanmean(allspectracon2wdg(:,(find(strcmp('C5',Campzc2wdg)==1)))'),'y') xlim([400 2500]) xlabel('Wavelength (nm)') ylabel('reflectance') %gtext('A') title('>0.5 WD Drought') subplot(323) plot(spec,(nanmean(allspectracon2wdl(:,(find(strcmp('C1',Campzc2wdl)==1)))'))-nanmean(allspectracon2wdl(:,(find(strcmp('C2',Campzc2wdl)==1)))'),'r') hold on plot(spec,(nanmean(allspectracon2wdl(:,(find(strcmp('C1',Campzc2wdl)==1)))'))-nanmean(allspectracon2wdl(:,(find(strcmp('C3',Campzc2wdl)==1)))'),'k') plot(spec,(nanmean(allspectracon2wdl(:,(find(strcmp('C1',Campzc2wdl)==1)))'))-nanmean(allspectracon2wdl(:,(find(strcmp('C4',Campzc2wdl)==1)))'),'g') plot(spec,(nanmean(allspectracon2wdl(:,(find(strcmp('C1',Campzc2wdl)==1)))'))-nanmean(allspectracon2wdl(:,(find(strcmp('C5',Campzc2wdl)==1)))'),'y') xlim([400 2500]) xlabel('Wavelength (nm)') ylabel('reflectance') title('<0.5 WD Drought') %gtext('B') subplot(322) plot(spec,(nanmean(allspectragrd2wdg(:,(find(strcmp('C1',Campzg2wdg)==1)))'))-nanmean(allspectragrd2wdg(:,(find(strcmp('C2',Campzg2wdg)==1)))'),'r') hold on plot(spec,(nanmean(allspectragrd2wdg(:,(find(strcmp('C1',Campzg2wdg)==1)))'))-nanmean(allspectragrd2wdg(:,(find(strcmp('C3',Campzg2wdg)==1)))'),'k') plot(spec,(nanmean(allspectragrd2wdg(:,(find(strcmp('C1',Campzg2wdg)==1)))'))-nanmean(allspectragrd2wdg(:,(find(strcmp('C4',Campzg2wdg)==1)))'),'g') plot(spec,(nanmean(allspectragrd2wdg(:,(find(strcmp('C1',Campzg2wdg)==1)))'))-nanmean(allspectragrd2wdg(:,(find(strcmp('C5',Campzg2wdg)==1)))'),'y') xlim([400 2500]) xlabel('Wavelength (nm)') ylabel('reflectance') title('>0.5 WD Girdle ') %gtext('D') subplot(324) plot(spec,(nanmean(allspectragrd2wdl(:,(find(strcmp('C1',Campzg2wdl)==1)))'))-nanmean(allspectragrd2wdl(:,(find(strcmp('C2',Campzg2wdl)==1)))'),'r') hold on plot(spec,(nanmean(allspectragrd2wdl(:,(find(strcmp('C1',Campzg2wdl)==1)))'))-nanmean(allspectragrd2wdl(:,(find(strcmp('C3',Campzg2wdl)==1)))'),'k') plot(spec,(nanmean(allspectragrd2wdl(:,(find(strcmp('C1',Campzg2wdl)==1)))'))-nanmean(allspectragrd2wdl(:,(find(strcmp('C4',Campzg2wdl)==1)))'),'g') plot(spec,(nanmean(allspectragrd2wdl(:,(find(strcmp('C1',Campzg2wdl)==1)))'))-nanmean(allspectragrd2wdl(:,(find(strcmp('C5',Campzg2wdl)==1)))'),'y') xlim([400 2500]) xlabel('Wavelength (nm)') ylabel('reflectance') title('<0.5 WD Girdle') %gtext('E') legend('C1-C2','C1-C3','C1-C4','C1-C5') %high wood den control allspectracon2wdgc1 = allspectracon2wdg(:,(find(strcmp('C1',Campzc2wdg)==1))); allspectracon2wdgc2 = allspectracon2wdg(:,(find(strcmp('C2',Campzc2wdg)==1))); allspectracon2wdgc3 = allspectracon2wdg(:,(find(strcmp('C3',Campzc2wdg)==1))); allspectracon2wdgc4 = allspectracon2wdg(:,(find(strcmp('C4',Campzc2wdg)==1))); allspectracon2wdgc5 = allspectracon2wdg(:,(find(strcmp('C5',Campzc2wdg)==1))); %low wood den allspectracon2wdlc1 = allspectracon2wdl(:,(find(strcmp('C1',Campzc2wdl)==1))); allspectracon2wdlc2 = allspectracon2wdl(:,(find(strcmp('C2',Campzc2wdl)==1))); allspectracon2wdlc3 = allspectracon2wdl(:,(find(strcmp('C3',Campzc2wdl)==1))); allspectracon2wdlc4 = allspectracon2wdl(:,(find(strcmp('C4',Campzc2wdl)==1))); allspectracon2wdlc5 = allspectracon2wdl(:,(find(strcmp('C5',Campzc2wdl)==1))); %significant differences bb=1; for bb = 1:length(spec) [H1,P] = ttest2(allspectracon2wdgc1(bb,:),allspectracon2wdlc1(bb,:)); cH1(bb) = H1; [H2,P] = ttest2(allspectracon2wdgc2(bb,:),allspectracon2wdlc2(bb,:)); cH2(bb) = H2; [H3,P] = ttest2(allspectracon2wdgc3(bb,:),allspectracon2wdlc3(bb,:)); cH3(bb) = H3; [H4,P] = ttest2(allspectracon2wdgc4(bb,:),allspectracon2wdlc4(bb,:)); cH4(bb) = H4; [H5,P] = ttest2(allspectracon2wdgc5(bb,:),allspectracon2wdlc5(bb,:)); cH5(bb) = H5; end cH1(cH1==0)=NaN; cH2(cH2==0)=NaN; cH3(cH3==0)=NaN; cH4(cH4==0)=NaN; cH5(cH5==0)=NaN; %high wood den girdle allspectragrd2wdgc1 = allspectragrd2wdg(:,(find(strcmp('C1',Campzg2wdg)==1))); allspectragrd2wdgc2 = allspectragrd2wdg(:,(find(strcmp('C2',Campzg2wdg)==1))); allspectragrd2wdgc3 = allspectragrd2wdg(:,(find(strcmp('C3',Campzg2wdg)==1))); allspectragrd2wdgc4 = allspectragrd2wdg(:,(find(strcmp('C4',Campzg2wdg)==1))); allspectragrd2wdgc5 = allspectragrd2wdg(:,(find(strcmp('C5',Campzg2wdg)==1))); %low wood den allspectragrd2wdlc1 = allspectragrd2wdl(:,(find(strcmp('C1',Campzg2wdl)==1))); allspectragrd2wdlc2 = allspectragrd2wdl(:,(find(strcmp('C2',Campzg2wdl)==1))); allspectragrd2wdlc3 = allspectragrd2wdl(:,(find(strcmp('C3',Campzg2wdl)==1))); allspectragrd2wdlc4 = allspectragrd2wdl(:,(find(strcmp('C4',Campzg2wdl)==1))); allspectragrd2wdlc5 = allspectragrd2wdl(:,(find(strcmp('C5',Campzg2wdl)==1))); %significant differences bb=1; for bb = 1:length(spec) [H1,P] = ttest2(allspectragrd2wdgc1(bb,:),allspectragrd2wdlc1(bb,:)); gH1(bb) = H1; [H2,P] = ttest2(allspectragrd2wdgc2(bb,:),allspectragrd2wdlc2(bb,:)); gH2(bb) = H2; [H3,P] = ttest2(allspectragrd2wdgc3(bb,:),allspectragrd2wdlc3(bb,:)); gH3(bb) = H3; [H4,P] = ttest2(allspectragrd2wdgc4(bb,:),allspectragrd2wdlc4(bb,:)); gH4(bb) = H4; [H5,P] = ttest2(allspectragrd2wdgc5(bb,:),allspectragrd2wdlc5(bb,:)); gH5(bb) = H5; end gH1(gH1==0)=NaN; gH2(gH2==0)=NaN; gH3(gH3==0)=NaN; gH4(gH4==0)=NaN; gH5(gH5==0)=NaN; %plot spectral changes by functional type high wood density minus low wood %density subplot(325) plot(spec,((nanmean(allspectracon2wdg(:,(find(strcmp('C1',Campzc2wdg)==1)))'))-nanmean(allspectracon2wdl(:,(find(strcmp('C1',Campzc2wdg)==1)))')).*cH1,'b') hold on plot(spec,((nanmean(allspectracon2wdg(:,(find(strcmp('C2',Campzc2wdg)==1)))'))-nanmean(allspectracon2wdl(:,(find(strcmp('C2',Campzc2wdg)==1)))')).*cH2,'r') plot(spec,((nanmean(allspectracon2wdg(:,(find(strcmp('C3',Campzc2wdg)==1)))'))-nanmean(allspectracon2wdl(:,(find(strcmp('C3',Campzc2wdg)==1)))')).*cH3,'k') plot(spec,((nanmean(allspectracon2wdg(:,(find(strcmp('C4',Campzc2wdg)==1)))'))-nanmean(allspectracon2wdl(:,(find(strcmp('C4',Campzc2wdg)==1)))')).*cH4,'g') plot(spec,((nanmean(allspectracon2wdg(:,(find(strcmp('C5',Campzc2wdg)==1)))'))-nanmean(allspectracon2wdl(:,(find(strcmp('C5',Campzc2wdg)==1)))')).*cH5,'y') xlim([400 2500]) xlabel('Wavelength (nm)') ylabel('reflectance') legend('C1', 'C2','C3','C4','C5') title('high density - low density drought') %gtext('C') subplot(326) plot(spec,((nanmean(allspectragrd2wdg(:,(find(strcmp('C1',Campzg2wdg)==1)))'))-nanmean(allspectragrd2wdl(:,(find(strcmp('C1',Campzg2wdl)==1)))')).*gH1,'b') hold on plot(spec,((nanmean(allspectragrd2wdg(:,(find(strcmp('C2',Campzg2wdg)==1)))'))-nanmean(allspectragrd2wdl(:,(find(strcmp('C2',Campzg2wdl)==1)))')).*gH1,'r') plot(spec,((nanmean(allspectragrd2wdg(:,(find(strcmp('C3',Campzg2wdg)==1)))'))-nanmean(allspectragrd2wdl(:,(find(strcmp('C3',Campzg2wdl)==1)))')).*gH1,'k') plot(spec,((nanmean(allspectragrd2wdg(:,(find(strcmp('C4',Campzg2wdg)==1)))'))-nanmean(allspectragrd2wdl(:,(find(strcmp('C4',Campzg2wdl)==1)))')).*gH1,'g') plot(spec,((nanmean(allspectragrd2wdg(:,(find(strcmp('C5',Campzg2wdg)==1)))'))-nanmean(allspectragrd2wdl(:,(find(strcmp('C5',Campzg2wdl)==1)))')).*gH1,'y') xlim([400 2500]) xlabel('Wavelength (nm)') ylabel('reflectance') title('high density - low density Girdle ') %gtext('F') % %plot photo and res allcon2C1 = (find(strcmp('C1',Campzc2)==1)); allcon2C2 = (find(strcmp('C2',Campzc2)==1)); allcon2C3 = (find(strcmp('C3',Campzc2)==1)); allcon2C4 = (find(strcmp('C4',Campzc2)==1)); allcon2C5 = (find(strcmp('C5',Campzc2)==1)); allgrd2C1 = (find(strcmp('C1',Campzg2)==1)); allgrd2C2 = (find(strcmp('C2',Campzg2)==1)); allgrd2C3 = (find(strcmp('C3',Campzg2)==1)); allgrd2C4 = (find(strcmp('C4',Campzg2)==1)); allgrd2C5 = (find(strcmp('C5',Campzg2)==1)); % sc = [ nanmean(lphc2(allcon2C1)); nanmean(lphc2(allcon2C2)); nanmean(lphc2(allcon2C3)); nanmean(lphc2(allcon2C4)); nanmean(lphc2(allcon2C5))]; % scg = [ nanmean(lphg2(allgrd2C1)); nanmean(lphg2(allgrd2C2)); nanmean(lphg2(allgrd2C3)); nanmean(lphg2(allgrd2C4)); nanmean(lphg2(allgrd2C5))]; % difzp = (sc-scg)./sc; % difz = (sc-scg); % plot(sc-difz(1)) % hold on % plot(scg,'r') % [H,Pgd] = ttest2(lphg2(allgrd2C5), lphc2(allcon2C5)-difz(1)); %significant change in girdle [H,Ppg] = ttest2(lphg2(allgrd2C1), lphg2(allgrd2C3)); [H,Ppc] = ttest2(lphc2(allcon2C1), lphg2(allcon2C3)); [H,Prg] = ttest2(lresg2(allgrd2C1), lresg2(allgrd2C3)); [H,Prc] = ttest2(lresc2(allcon2C1), lresg2(allcon2C3)); % gtext([' P=', num2str(P)]); %remove seasonal cyce using control data figure(3) subplot(321) errorbar([nanmean(lphg2(allgrd2C1)), nanmean(lphg2(allgrd2C2)), nanmean(lphg2(allgrd2C3)), nanmean(lphg2(allgrd2C4)), nanmean(lphg2(allgrd2C5))],[nanstd(lphg2((allgrd2C1)))/sqrt(20), nanstd(lphg2((allgrd2C2)))/sqrt(30), nanstd(lphg2((allgrd2C3)))/sqrt(20), nanstd(lphg2((allgrd2C4)))/sqrt(20), nanstd(lphg2((allgrd2C5)))/sqrt(20)],'bo-') hold on errorbar([nanmean(lphc2(allcon2C1))-2.19, nanmean(lphc2(allcon2C2))-2.19, nanmean(lphc2(allcon2C3))-2.19, nanmean(lphc2(allcon2C4))-2.19, nanmean(lphc2(allcon2C5))-2.19],[nanstd(lphc2((allcon2C1)))/sqrt(20), nanstd(lphc2((allcon2C2)))/sqrt(20), nanstd(lphc2((allcon2C3)))/sqrt(20), nanstd(lphc2((allcon2C4)))/sqrt(20), nanstd(lphc2((allcon2C5)))/sqrt(20)],'ro-') %bar([nanmean(lphg2(allgrd2C1)) nanmean(lphc2(allcon2C1)); nanmean(lphg2(allgrd2C2)) nanmean(lphc2(allcon2C2)); nanmean(lphg2(allgrd2C3)) nanmean(lphc2(allcon2C3)); nanmean(lphg2(allgrd2C4)) nanmean(lphc2(allcon2C4)); nanmean(lphg2(allgrd2C5)) nanmean(lphc2(allcon2C5))]) xlabel('Campaign') ylabel('Asat umol m^-^2 s^-^1') set(gca, 'XTick', 0:5) ylim([0 10]) legend('girdle','drought') %gtext('A') subplot(322) errorbar([nanmean(lresg2(find(strcmp('C1',Campzg2)==1))), nanmean(lresg2(find(strcmp('C2',Campzg2)==1))), nanmean(lresg2(find(strcmp('C3',Campzg2)==1))), nanmean(lresg2(find(strcmp('C4',Campzg2)==1))), nanmean(lresg2(find(strcmp('C5',Campzg2)==1)))],[nanstd(lresg2(find(strcmp('C1',Campzg2)==1)))/sqrt(20), nanstd(lresg2(find(strcmp('C2',Campzg2)==1)))/sqrt(20), nanstd(lresg2(find(strcmp('C3',Campzg2)==1)))/sqrt(20), nanstd(lresg2(find(strcmp('C4',Campzg2)==1)))/sqrt(20), nanstd(lresg2(find(strcmp('C5',Campzg2)==1)))/sqrt(20)],'bo-') hold on %bar([nanmean(lresg2(allgrd2C1)) nanmean(lresc2(allcon2C1)); nanmean(lresg2(allgrd2C2)) nanmean(lresc2(allcon2C2)); nanmean(lresg2(allgrd2C3)) nanmean(lresc2(allcon2C3)); nanmean(lresg2(allgrd2C4)) nanmean(lresc2(allcon2C4)); nanmean(lresg2(allgrd2C5)) nanmean(lresc2(allcon2C5))]) errorbar([nanmean(lresc2(find(strcmp('C1',Campzc2)==1))), nanmean(lresc2(find(strcmp('C2',Campzc2)==1))), nanmean(lresc2(find(strcmp('C3',Campzc2)==1))), nanmean(lresc2(find(strcmp('C4',Campzc2)==1))), nanmean(lresc2(find(strcmp('C5',Campzc2)==1)))],[nanstd(lresc2(find(strcmp('C1',Campzc2)==1)))/sqrt(20), nanstd(lresc2(find(strcmp('C2',Campzc2)==1)))/sqrt(20), nanstd(lresc2(find(strcmp('C3',Campzc2)==1)))/sqrt(20), nanstd(lresc2(find(strcmp('C4',Campzc2)==1)))/sqrt(20), nanstd(lresc2(find(strcmp('C5',Campzc2)==1)))/sqrt(20)],'ro-') xlabel('Campaign') ylabel('Rdark umol m^-^2 s^-^1') set(gca, 'XTick', 0:5) %gtext('B') subplot(323) plot([nanmean(lresg2(find(strcmp('C1',Campzg2)==1))), nanmean(lresg2(find(strcmp('C2',Campzg2)==1))), nanmean(lresg2(find(strcmp('C3',Campzg2)==1))), nanmean(lresg2(find(strcmp('C4',Campzg2)==1))), nanmean(lresg2(find(strcmp('C5',Campzg2)==1)))]./[nanmean(lphg2(allgrd2C1)), nanmean(lphg2(allgrd2C2)), nanmean(lphg2(allgrd2C3)), nanmean(lphg2(allgrd2C4)), nanmean(lphg2(allgrd2C5))],'bo-') hold on plot([nanmean(lresc2(find(strcmp('C1',Campzc2)==1))), nanmean(lresc2(find(strcmp('C2',Campzc2)==1))), nanmean(lresc2(find(strcmp('C3',Campzc2)==1))), nanmean(lresc2(find(strcmp('C4',Campzc2)==1))), nanmean(lresc2(find(strcmp('C5',Campzc2)==1)))]./[nanmean(lphc2(allcon2C1))-2.19, nanmean(lphc2(allcon2C2))-2.19, nanmean(lphc2(allcon2C3))-2.19, nanmean(lphc2(allcon2C4))-2.19, nanmean(lphc2(allcon2C5))-2.19],'ro-') xlabel('Campaign') ylabel('Rdark/Asat') set(gca, 'XTick', 0:5) %gtext('C') subplot(324) errorbar([nanmean(lLMAg2(find(strcmp('C1',Campzg2)==1))), nanmean(lLMAg2(find(strcmp('C2',Campzg2)==1))), nanmean(lLMAg2(find(strcmp('C3',Campzg2)==1))), nanmean(lLMAg2(find(strcmp('C4',Campzg2)==1))), nanmean(lLMAg2(find(strcmp('C5',Campzg2)==1)))],[nanstd(lLMAg2(find(strcmp('C1',Campzg2)==1)))/sqrt(20), nanstd(lLMAg2(find(strcmp('C2',Campzg2)==1)))/sqrt(20), nanstd(lLMAg2(find(strcmp('C3',Campzg2)==1)))/sqrt(333), nanstd(lLMAg2(find(strcmp('C4',Campzg2)==1)))/sqrt(20), nanstd(lLMAg2(find(strcmp('C5',Campzg2)==1)))/sqrt(20)],'bo-') hold on errorbar([nanmean(lLMAc2(find(strcmp('C1',Campzc2)==1))), nanmean(lLMAc2(find(strcmp('C2',Campzc2)==1))), nanmean(lLMAc2(find(strcmp('C3',Campzc2)==1))), nanmean(lLMAc2(find(strcmp('C4',Campzc2)==1))), nanmean(lLMAc2(find(strcmp('C5',Campzc2)==1)))],[nanstd(lLMAc2(find(strcmp('C1',Campzc2)==1)))/sqrt(20), nanstd(lLMAc2(find(strcmp('C2',Campzc2)==1)))/sqrt(20), nanstd(lLMAc2(find(strcmp('C3',Campzc2)==1)))/sqrt(20), nanstd(lLMAc2(find(strcmp('C4',Campzc2)==1)))/sqrt(20), nanstd(lLMAc2(find(strcmp('C5',Campzc2)==1)))/sqrt(20)],'ro-') %bar([nanmean(lLMAg2(allgrd2C1)) nanmean(lLMAc2(allcon2C1)); nanmean(lLMAg2(allgrd2C2)) nanmean(lLMAc2(allcon2C2)); nanmean(lLMAg2(allgrd2C3)) nanmean(lLMAc2(allcon2C3)); nanmean(lLMAg2(allgrd2C4)) nanmean(lLMAc2(allcon2C4)); nanmean(lLMAg2(allgrd2C5)) nanmean(lLMAc2(allcon2C5))]) set(gca, 'XTick', 0:5) xlabel('Campaign') ylabel('LMA g m^-^2') %gtext('D') subplot(325) errorbar([nanmean(lwatercong2(find(strcmp('C1',Campzg2)==1))), nanmean(lwatercong2(find(strcmp('C2',Campzg2)==1))), nanmean(lwatercong2(find(strcmp('C3',Campzg2)==1))), nanmean(lwatercong2(find(strcmp('C4',Campzg2)==1))), nanmean(lwatercong2(find(strcmp('C5',Campzg2)==1)))],[nanstd(lwatercong2(find(strcmp('C1',Campzg2)==1)))/sqrt(20), nanstd(lwatercong2(find(strcmp('C2',Campzg2)==1)))/sqrt(20), nanstd(lwatercong2(find(strcmp('C3',Campzg2)==1)))/sqrt(333), nanstd(lwatercong2(find(strcmp('C4',Campzg2)==1)))/sqrt(20), nanstd(lwatercong2(find(strcmp('C5',Campzg2)==1)))/sqrt(20)],'bo-') hold on errorbar([nanmean(lwaterconc2(find(strcmp('C1',Campzc2)==1))), nanmean(lwaterconc2(find(strcmp('C2',Campzc2)==1))), nanmean(lwaterconc2(find(strcmp('C3',Campzc2)==1))), nanmean(lwaterconc2(find(strcmp('C4',Campzc2)==1))), nanmean(lwaterconc2(find(strcmp('C5',Campzc2)==1)))],[nanstd(lwaterconc2(find(strcmp('C1',Campzc2)==1)))/sqrt(20), nanstd(lwaterconc2(find(strcmp('C2',Campzc2)==1)))/sqrt(20), nanstd(lwaterconc2(find(strcmp('C3',Campzc2)==1)))/sqrt(20), nanstd(lwaterconc2(find(strcmp('C4',Campzc2)==1)))/sqrt(20), nanstd(lwaterconc2(find(strcmp('C5',Campzc2)==1)))/sqrt(20)],'ro-') %bar([nanmean(lLMAg2(allgrd2C1)) nanmean(lLMAc2(allcon2C1)); nanmean(lLMAg2(allgrd2C2)) nanmean(lLMAc2(allcon2C2)); nanmean(lLMAg2(allgrd2C3)) nanmean(lLMAc2(allcon2C3)); nanmean(lLMAg2(allgrd2C4)) nanmean(lLMAc2(allcon2C4)); nanmean(lLMAg2(allgrd2C5)) nanmean(lLMAc2(allcon2C5))]) set(gca, 'XTick', 0:5) xlabel('Campaign') ylabel('Water Content %') %gtext('E') subplot(326) errorbar([ndvi1 ndvi2 ndvi3 ndvi4 ndvi5],([ndvi1s ndvi2s ndvi3s ndvi4s ndvi5s])./sqrt(200),'bo-'); hold on errorbar([ndvig1 ndvig2 ndvig3 ndvig4 ndvig5],([ndvig1s ndvig2s ndvig3s ndvig4s ndvig5s])./sqrt(200),'ro-'); set(gca, 'XTick', 0:5) xlabel('Campaign') ylabel('NDVI') %gtext('F') %legend('Drought','Girdle') i % do a repeated measures anova on these results % rm = [(lresc2(find(strcmp('C1',Campzc2)==1))), (lresc2(find(strcmp('C2',Campzc2)==1))), (lresc2(find(strcmp('C3',Campzc2)==1))), (lresc2(find(strcmp('C4',Campzc2)==1))), (lresc2(find(strcmp('C5',Campzc2)==1)))] % ranovatbl = ranova(rm) %plot drought spectra % find statistical differences allspectracon2C1 = allspectracon2(:,(find(strcmp('C1',Campzc2)==1))); allspectracon2C3 = allspectracon2(:,(find(strcmp('C3',Campzc2)==1))); allspectracon2C5 = allspectracon2(:,(find(strcmp('C5',Campzc2)==1))); for bb = 1:length(spec) [H1,P] = ttest2(allspectracon2C1(bb,:),allspectracon2C3(bb,:)); sH1(bb) = H1; [H2,P] = ttest2(allspectracon2C3(bb,:),allspectracon2C5(bb,:)); sH2(bb) = H2; [H3,P] = ttest2(allspectracon2C1(bb,:),allspectracon2C5(bb,:)); sH3(bb) = H3; end sH1(sH1==0)=NaN; sH2(sH2==0)=NaN; sH3(sH3==0)=NaN; % plot(spec,((nanmean(allspectracon2C1')-nanmean(allspectracon2C3')).*sH1),'k') % hold on % plot(spec,((nanmean(allspectracon2C3')-nanmean(allspectracon2C5')).*sH2),'r') % plot(spec,((nanmean(allspectracon2C1')-nanmean(allspectracon2C5')).*sH3),'g') % xlim([400 2500]) % xlabel('spectra') % ylabel('reflectance') % title('Drought spectra ') %The traits tree that died in the control site is tag number 134. %find data for tree 134, Tree died on March 23rd %Campaign 1=21-Jan-16, 2=10-Feb-16, 3=01-Mar-16, 4=29-Mar-16 fdt = find(treezc2==134); allspectracon2t134 = allspectracon2(:,fdt); Campzc2t134 = Campzc2(fdt); lphc2t134 = lphc2(fdt); lresc2t134 = lresc2(fdt); lLMAc2t134 = lLMAc2(fdt); lwaterconc2t134 = lwaterconc2(fdt); allspectracon2t134C1 = allspectracon2t134(:,(find(strcmp('C1',Campzc2t134)==1))); allspectracon2t134C2 = allspectracon2t134(:,(find(strcmp('C2',Campzc2t134)==1))); allspectracon2t134C3 = allspectracon2t134(:,(find(strcmp('C3',Campzc2t134)==1))); allspectracon2t134C4 = allspectracon2t134(:,(find(strcmp('C4',Campzc2t134)==1))); cha = allspectracon2t134C4(970,1)-allspectracon2t134C4(971,1); for x=1:14 cha = allspectracon2t134C4(621,1)-allspectracon2t134C4(622,1); allspectracon2t134C4n(:,x)=[allspectracon2t134C4(1:621,x); allspectracon2t134C4(622:end,x)+cha]; end for bb = 1:length(spec) [H1,P] = ttest2(allspectracon2t134C1(bb,:),allspectracon2t134C3(bb,:)); sH1(bb) = H1; [H2,P] = ttest2(allspectracon2t134C1(bb,:),allspectracon2t134C4n(bb,:)); sH2(bb) = H2; end sH1(sH1==0)=NaN; sH2(sH2==0)=NaN; figure(4) %subplot(211) %plot(spec,((nanmean(allspectracon2t134C1')-nanmean(allspectracon2t134C3')).*sH1),'ro') %plot death spectra % find statistical differences for bb = 1:length(spec) [H,P] = ttest2(allspectragrd2(bb,id),allspectragrd2(bb,nid)); sH(bb) = H; end sH(sH==0)=NaN; % subplot(211) % plot(spec,nanmean(allspectragrd2(:,id)')) % hold on % plot(spec,nanmean(allspectragrd2(:,nid)'),'r') % % plot(nanmean(allspectracon2'),'g') % xlabel('spectra') % ylabel('reflectance') % title('Girdle') %nid >50days from death %subplot(212) plot(spec,(nanmean(allspectragrd2(:,nid)')-nanmean(allspectragrd2(:,id)')), 'k') hold on plot(spec,((nanmean(allspectracon2t134C1')-nanmean(allspectracon2t134C4n'))),'r') plot(spec,(nanmean(allspectragrd2(:,nid)')-nanmean(allspectragrd2(:,id)')).*sH, 'ko') plot(spec,((nanmean(allspectracon2t134C1')-nanmean(allspectracon2t134C4n')).*sH2),'ro') xlabel('Wavelength (nm)') ylabel('reflectance') %title('tree mortality') xlim([400 2500]) legend('girdle','drought') i figure(5) % % plot drought death physiology subplot(421) bar([nanmean(lphc2t134(find(strcmp('C1',Campzc2t134)==1))), nanmean(lphc2t134(find(strcmp('C3',Campzc2t134)==1)))]) [H,P] = ttest2((lphc2t134(find(strcmp('C1',Campzc2t134)==1))), (lphc2t134(find(strcmp('C3',Campzc2t134)==1)))); P=round(P,3); % gtext('A') % gtext([' P=', num2str(P)]); xlabel('drought death') ylabel('Asat umol m^-^2 s^-^1') subplot(423) bar([nanmean(lresc2t134(find(strcmp('C1',Campzc2t134)==1))), nanmean(lresc2t134(find(strcmp('C3',Campzc2t134)==1)))]) [H,P1] = ttest2((lresc2t134(find(strcmp('C1',Campzc2t134)==1))), (lresc2t134(find(strcmp('C3',Campzc2t134)==1)))); P1=round(P1,3); % gtext('C') % gtext([' P=', num2str(P1)]); xlabel('drought death') ylabel('Rdark umol m^-^2 s^-^1') subplot(425) bar([nanmean(lresc2t134(find(strcmp('C1',Campzc2t134)==1))), nanmean(lresc2t134(find(strcmp('C3',Campzc2t134)==1)))]./[nanmean(lphc2t134(find(strcmp('C1',Campzc2t134)==1))), nanmean(lphc2t134(find(strcmp('C3',Campzc2t134)==1)))]) dtl = lresc2t134./lphc2t134; [H,P1] = ttest2((dtl(find(strcmp('C1',Campzc2t134)==1))), (dtl(find(strcmp('C3',Campzc2t134)==1)))); P1=round(P1,3); %gtext('E') %gtext([' P=', num2str(P1)]); xlabel('drought death') ylabel('Rdark/Asat') subplot(427) bar([nanmean(lLMAc2t134(find(strcmp('C1',Campzc2t134)==1))), nanmean(lLMAc2t134(find(strcmp('C3',Campzc2t134)==1)))]) [H,P3] = ttest2((lLMAc2t134(find(strcmp('C1',Campzc2t134)==1))), (lLMAc2t134(find(strcmp('C3',Campzc2t134)==1)))); P3=round(P3,3); % gtext('G') % gtext([' P=', num2str(P3)]); xlabel('drought death') ylabel('LMA g m^-^2') subplot(422) bar([nanmean(lphg2(nid)), nanmean(lphg2(id))]) % hold on % errorbar([nanmean(lphg2(nid)), nanmean(lphg2(id))],[nanstd(lphg2(id))/sqrt(122), nanstd(lphg2(nid))/sqrt(675)], 'ro') [H,P] = ttest2(lphg2(id), lphg2(nid)); P=round(P,3); % gtext('B') % gtext([' P=', num2str(P)]); xlabel('girdle death') ylabel('Asat umol m^-^2 s^-^1') subplot(424) bar([nanmean(lresg2(nid)), nanmean(lresg2(id))]) % hold on % errorbar([nanmean(lresg2(id)), nanmean(lresg2(nid))],[nanstd(lresg2(id))/sqrt(122), nanstd(lresg2(nid))/sqrt(675)], 'ro') [H,P1] = ttest2(lresg2(id), lresg2(nid)); P1=round(P1,3); % gtext('D') % gtext([' P=', num2str(P1)]); xlabel('girdle death') ylabel('Rdark umol m^-^2 s^-^1') subplot(426) bar([nanmean(lresg2(nid)), nanmean(lresg2(id))]./[nanmean(lphg2(nid)), nanmean(lphg2(id))]) % hold on % errorbar([nanmean(lresg2(id)), nanmean(lresg2(nid))],[nanstd(lresg2(id))/sqrt(122), nanstd(lresg2(nid))/sqrt(675)], 'ro') adres = lresg2./lphg2; adres(find(adres>10))=NaN; adres(find(adres<-10))=NaN; [H,P1] = ttest2(adres(id), adres(nid)); P1=round(P1,3); % gtext('F') % gtext([' P=', num2str(P1)]); xlabel('girdle death') ylabel('Rdark/Asat') subplot(428) bar([nanmean(lLMAg2(nid)), nanmean(lLMAg2(id))]) % hold on % errorbar([nanmean(lLMAg2(id)), nanmean(lLMAg2(nid))],[nanstd(lLMAg2(id))/sqrt(122), nanstd(lLMAg2(nid))/sqrt(675)], 'ro') [H,P3] = ttest2(lLMAg2(id), lLMAg2(nid)); P3=round(P3,3); % gtext('H') % gtext([' P=', num2str(P3)]); xlabel('girdle death') ylabel('LMA g m^-^2') figure(6) subplot(121) bar([nanmean(lwaterconc2t134(find(strcmp('C1',Campzc2t134)==1))), nanmean(lwaterconc2t134(find(strcmp('C3',Campzc2t134)==1)))]) [H,P3] = ttest2((lwaterconc2t134(find(strcmp('C1',Campzc2t134)==1))), (lwaterconc2t134(find(strcmp('C3',Campzc2t134)==1)))); P3=round(P3,3); % gtext('I') % gtext([' P=', num2str(P3)]); xlabel('drought death') ylabel('% water') %plot near death physiology subplot(122) bar([nanmean(lwatercong2(nid)), nanmean(lwatercong2(id))]) % hold on % errorbar([nanmean(lwatercong2(id)), nanmean(lwatercong2(nid))],[nanstd(lwatercong2(id))/sqrt(122), nanstd(lwatercong2(nid))/sqrt(675)], 'ro') [H,P3] = ttest2(lwatercong2(id), lwatercong2(nid)); P3=round(P3,3); % gtext('J') % gtext([' P=', num2str(P3)]); xlabel('girdle death') ylabel('% water') i % % % plot drought physiology % C1cz = (find(strcmp('C1',Campzc2)==1)); % C3cz = (find(strcmp('C3',Campzc2)==1)); % C5cz = (find(strcmp('C5',Campzc2)==1)); % subplot(322) % bar([nanmean(lphc2(C1cz)), nanmean(lphc2(C3cz))]) % hold on % errorbar([nanmean(lphc2(C1cz)), nanmean(lphc2(C3cz))],[nanstd(lphc2(C1cz))/sqrt(122), nanstd(lphc2(C3cz))/sqrt(675)], 'ro') % [H,P] = ttest2(lphc2(C1cz), lphc2(C3cz)); % gtext([' P=', num2str(P)]); % xlabel('drought') % ylabel('photosynthesis') % subplot(324) % bar([nanmean(lresc2(C1cz)), nanmean(lresc2(C3cz))]) % hold on % errorbar([nanmean(lresc2(C1cz)), nanmean(lresc2(C3cz))],[nanstd(lresc2(C1cz))/sqrt(122), nanstd(lresc2(C3cz))/sqrt(675)], 'ro') % [H,P1] = ttest2(lresc2(C1cz), lresc2(C3cz)); % gtext([' P=', num2str(P1)]); % xlabel('drought') % ylabel('Respiration') % subplot(326) % bar([nanmean(lLMAc2(C1cz)), nanmean(lLMAc2(C3cz))]) % hold on % errorbar([nanmean(lLMAc2(C1cz)), nanmean(lLMAc2(C3cz))],[nanstd(lLMAc2(C1cz))/sqrt(122), nanstd(lLMAc2(C3cz))/sqrt(675)], 'ro') % [H,P3] = ttest2(lLMAc2(C1cz), lLMAc2(C3cz)); % gtext([' P=', num2str(P3)]); % xlabel('drought') % ylabel('LMA') i %plot(nanmean(allspectragrd2(:,id)')-nanmean(allspectracon2')) %%do the PLSRs spec2 = 401:2500; %%%water % %remove NaNs from variables vtzw = ~isnan(lwaterconc2); vtestw = lwaterconc2(vtzw); spectestw = allspectracon2(52:end,vtzw); %remove outliers ivws = find(vtestw>0.4); vtestw2s = vtestw(ivws); spectestw2s = spectestw(:,ivws); ivwx = find(vtestw2s<0.9); vtestw2 = vtestw2s(ivwx); spectestw2 = spectestw2s(:,ivwx); %divide data to have 30% for validation vv=1; vv2=1; for j=1:length(vtestw2) if rand(1)>0.7 vtestw3(vv)=vtestw2(j); spectestw3(:,vv) = spectestw2(:,j); vv=vv+1; else vtestw4(vv2)=vtestw2(j); spectestw4(:,vv2) = spectestw2(:,j); vv2=vv2+1; end end %do the PLSR len = length(vtestw4); len2 = 30; [xl,yl,xs,ys,beta,pctvar,mse,stats] = plsregress(spectestw4',vtestw4',len2,'CV',len); dmin = find(mse(2,:)==min(mse(2,:))) RMSEl = sqrt(mse(2,dmin)) RMSEmeanl = round(sqrt(mse(2,dmin))/mean(vtestw4),2) RMSEstdl = round(sqrt(mse(2,dmin))/std(vtestw4),2) [xl,yl,xs,ys,beta,pctvar,mse,statz] = plsregress(spectestw4',vtestw4',dmin,'CV',len); d = [ones(size(spectestw4',1),1) spectestw4']*beta; dt = [ones(size(spectestw3',1),1) spectestw3']*beta; W1 = statz.W; W1=W1'; W1w = W1; figure(7) subplot(5,1,1) plot(spec2,W1(1,:)*pctvar(2,1)) hold on %plot(W1(2,:)*pctvar(2,2),'r') xlim([400 2500]) xlabel('Wavelength (nm)') ylabel('PLS Weightings') gtext('B') R12 = corrcoef(d,vtestw4'); R12t = corrcoef(dt,vtestw3'); r12 = R12(1,2); % Correlation coefficient r12t = R12t(1,2); r12sq = round(r12^2,2) % Coefficient of determination r12sqt = round(r12t^2, 2) % %r2 = stats(1) figure(8) subplot(5,1,1) plot(vtestw4,d,'.') hold on plot(vtestw3,dt,'r.') plot([min(d) ; max(d) ], [min(d) ;max(d) ]) xlabel('measured') ylabel('predicted') title('% water' ) gtext('A') axis square % gtext(['r^2=', num2str(r12sq) ' RMSE=', num2str(RMSEmeanl)]); %%LMA % %remove NaNs from variables vtzl = ~isnan(lLMAc2); vtestl = lLMAc2(vtzl); spectestl = allspectracon2(52:end,vtzl); %remove outliers ivl = find(vtestl>20); vtestl2 = vtestl(ivl); spectestl2 = spectestl(:,ivl); %divide data to have 30% for validation vv=1; vv2=1; for j=1:length(vtestl2) if rand(1)>0.7 vtestl3(vv)=vtestl2(j); spectestl3(:,vv) = spectestl2(:,j); vv=vv+1; else vtestl4(vv2)=vtestl2(j); spectestl4(:,vv2) = spectestl2(:,j); vv2=vv2+1; end end %do the PLSR len = length(vtestl4); len2 = 30; [xl,yl,xs,ys,beta,pctvar,mse,stats] = plsregress(spectestl4',vtestl4',len2,'CV',len); dmin = find(mse(2,:)==min(mse(2,:))) RMSEl = sqrt(mse(2,dmin)) RMSEmeanl = round(sqrt(mse(2,dmin))/mean(vtestl4),2) RMSEstdl = round(sqrt(mse(2,dmin))/std(vtestl4),2) [xl,yl,xs,ys,beta,pctvar,mse,statz] = plsregress(spectestl4',vtestl4',dmin,'CV',len); d = [ones(size(spectestl4',1),1) spectestl4']*beta; dt = [ones(size(spectestl3',1),1) spectestl3']*beta; W1 = statz.W; W1=W1'; W1l = W1; figure(7) subplot(5,1,2) plot(spec2,W1(1,:)*pctvar(2,1)) hold on %plot(W1(2,:)*pctvar(2,2),'r') xlim([400 2500]) xlabel('Wavelength (nm)') ylabel('PLS Weightings') gtext('D') R12 = corrcoef(d,vtestl4'); R12t = corrcoef(dt,vtestl3'); r12 = R12(1,2); % Correlation coefficient r12t = R12t(1,2); r12sq = round(r12^2,2) % Coefficient of determination r12sqt = round(r12t^2, 2) % %r2 = stats(1) figure(8) subplot(5,1,2) plot(vtestl4,d,'.') hold on plot(vtestl3,dt,'r.') plot([min(d) ; max(d) ], [min(d) ;max(d) ]) xlabel('measured') ylabel('predicted') title('LMA g m^-^2' ) gtext('C') axis square %gtext(['r^2=', num2str(r12sq) ' RMSE=', num2str(RMSEmeanl)]); % %photo % %remove NaNs from variables vtzp = ~isnan(lphc2); vtestp = lphc2(vtzp); spectestp = allspectracon2(52:end,vtzp); %remove outliers ivp = find(vtestp>0); vtestp2 = vtestp(ivp); spectestp2 = spectestp(:,ivp); %divide data to have 30% for validation vv=1; vv2=1; for j=1:length(vtestp2) if rand(1)>0.7 vtestp3(vv)=vtestp2(j); spectestp3(:,vv) = spectestp2(:,j); vv=vv+1; else vtestp4(vv2)=vtestp2(j); spectestp4(:,vv2) = spectestp2(:,j); vv2=vv2+1; end end %do the PLSR len = length(vtestp4); len2 = 30; [xl,yl,xs,ys,beta,pctvar,mse,stats] = plsregress(spectestp4',vtestp4',len2,'CV',len); dmin = find(mse(2,:)==min(mse(2,:))) RMSEl = sqrt(mse(2,dmin)) RMSEmeanl = round(sqrt(mse(2,dmin))/mean(vtestp4),2) RMSEstdl = round(sqrt(mse(2,dmin))/std(vtestp4),2) [xl,yl,xs,ys,beta,pctvar,mse,statz] = plsregress(spectestp4',vtestp4',dmin,'CV',len); d = [ones(size(spectestp4',1),1) spectestp4']*beta; dt = [ones(size(spectestp3',1),1) spectestp3']*beta; W1 = statz.W; W1=W1'; W1p=W1; figure(7) subplot(5,1,3) plot(spec2,W1(1,:)*pctvar(2,1)) hold on xlim([400 2500]) %plot(W1(2,:)*pctvar(2,2),'r') xlabel('Wavelength (nm)') ylabel('PLS Weightings') gtext('F') R12 = corrcoef(d,vtestp4'); R12t = corrcoef(dt,vtestp3'); r12 = R12(1,2); % Correlation coefficient r12t = R12t(1,2); r12sq = round(r12^2,2) % Coefficient of determination r12sqt = round(r12t^2, 2) % %r2 = stats(1) figure(8) subplot(5,1,3) plot(vtestp4,d,'.') hold on plot(vtestp3,dt,'r.') plot([min(d) ; max(d) ], [min(d) ;max(d) ]) xlabel('measured') ylabel('predicted') title('A umol m^-^2 s^-^1' ) gtext('E') axis square %gtext(['r^2=', num2str(r12sq) ' RMSE=', num2str(RMSEmeanl)]); % %wood den % %remove NaNs from variables vtzw = ~isnan(wooddenzc2); vtestw = wooddenzc2(vtzp); spectestw = allspectracon2(52:end,vtzw); %remove outliers ivw = find(vtestw>0.2); vtestw2 = vtestw(ivw); spectestw2 = spectestw(:,ivw); %divide data to have 30% for validation vv=1; vv2=1; for j=1:length(vtestw2) if rand(1)>0.7 vtestw3(vv)=vtestw2(j); spectestw3(:,vv) = spectestw2(:,j); vv=vv+1; else vtestw4(vv2)=vtestw2(j); spectestw4(:,vv2) = spectestw2(:,j); vv2=vv2+1; end end %do the PLSR len = length(vtestw4); len2 = 30; [xl,yl,xs,ys,beta,pctvar,mse,stats] = plsregress(spectestw4',vtestw4',len2,'CV',len); dmin = find(mse(2,:)==min(mse(2,:))) RMSEl = sqrt(mse(2,dmin)) RMSEmeanl = round(sqrt(mse(2,dmin))/mean(vtestw4),2) RMSEstdl = round(sqrt(mse(2,dmin))/std(vtestw4),2) [xl,yl,xs,ys,beta,pctvar,mse,statz] = plsregress(spectestw4',vtestw4',dmin,'CV',len); d = [ones(size(spectestw4',1),1) spectestw4']*beta; dt = [ones(size(spectestw3',1),1) spectestw3']*beta; W1 = statz.W; W1=W1'; W1w=W1; figure(7) subplot(5,1,4) plot(spec2,W1(1,:)*pctvar(2,1)) hold on %plot(W1(2,:)*pctvar(2,2),'r') xlim([400 2500]) xlabel('Wavelength (nm)') ylabel('PLS Weightings') gtext('H') R12 = corrcoef(d,vtestw4'); R12t = corrcoef(dt,vtestw3'); r12 = R12(1,2); % Correlation coefficient r12t = R12t(1,2); r12sq = round(r12^2,2) % Coefficient of determination r12sqt = round(r12t^2, 2) % %r2 = stats(1) figure(8) subplot(5,1,4) plot(vtestw4,d,'.') hold on plot(vtestw3,dt,'r.') plot([min(d) ; max(d) ], [min(d) ;max(d) ]) xlabel('measured') ylabel('predicted') title(' wood density (g cm^-^3)' ) gtext('G') axis square %gtext(['r^2=', num2str(r12sq) ' RMSE=', num2str(RMSEmeanl)]); %death %remove NaNs from variables vtzd = ~isnan(diffdeathg2); vtestd = diffdeathg2(vtzd); spectestd = allspectracon2(52:end,vtzd); %remove outliers ivd = find(vtestd>0); vtestd2 = vtestd(ivd); spectestd2 = spectestd(:,ivd); %divide data to have 30% for validation vv=1; vv2=1; for j=1:length(vtestd2) if rand(1)>0.7 vtestd3(vv)=vtestd2(j); spectestd3(:,vv) = spectestd2(:,j); vv=vv+1; else vtestd4(vv2)=vtestd2(j); spectestd4(:,vv2) = spectestd2(:,j); vv2=vv2+1; end end %do the PLSR len = length(vtestd4); len2 = 30; [xl,yl,xs,ys,beta,pctvar,mse,stats] = plsregress(spectestd4',vtestd4',len2,'CV',len); dmin = find(mse(2,:)==min(mse(2,:))) RMSEl = sqrt(mse(2,dmin)) RMSEmeanl = round(sqrt(mse(2,dmin))/mean(vtestd4),2) RMSEstdl = round(sqrt(mse(2,dmin))/std(vtestd4),2) [xl,yl,xs,ys,beta,pctvar,mse,statz] = plsregress(spectestd4',vtestd4',dmin,'CV',len); d = [ones(size(spectestd4',1),1) spectestd4']*beta; dt = [ones(size(spectestd3',1),1) spectestd3']*beta; W1 = statz.W; W1=W1'; W1d=W1; figure(7) subplot(5,1,5) plot(spec2, W1(1,:)*pctvar(2,1)) hold on %plot(W1(2,:)*pctvar(2,2),'r') xlim([400 2500]) xlabel('Wavelength (nm)') ylabel('PLS Weightings') gtext('J') R12 = corrcoef(d,vtestd4'); R12t = corrcoef(dt,vtestd3'); r12 = R12(1,2); % Correlation coefficient r12t = R12t(1,2); r12sq = round(r12^2,2) % Coefficient of determination r12sqt = round(r12t^2, 2) % %r2 = stats(1) figure(8) subplot(5,1,5) plot(vtestd4,d,'.') hold on plot(vtestd3,dt,'r.') plot([min(d) ; max(d) ], [min(d) ;max(d) ]) xlabel('measured') ylabel('predicted') title('time to mortality (days)' ) xlim([0 400]) ylim([0 400]) axis square gtext('I') %gtext(['r^2=', num2str(r12sq) ' RMSE=', num2str(RMSEmeanl)]);