;Code for the WLEF photosynthesis paper PRO wlef_smget,lat,lon,name IF n_elements(lat) EQ 0 THEN BEGIN lat = 45.94587778 lon = -90.27230417 name = 'wlef' ENDIF ;lon -180 to 180 (1440) ;lat -90 to 90 (720) ;sm(lat,lon) ;0.25 degree indir = '/sap/projects/psyn/soilmoist/' smdata = make_Array(366,32,/float,value=nan()) FOR yr = 1979,2010 DO BEGIN print,yr ystr = string(yr,format='(i4.4)') thedir = indir + ystr + '/' doy = 1l FOR mo = 1,12 DO BEGIN mstr = string(mo,format='(i2.2)') FOR day = 1,days_in_month(mo,y=yr) DO BEGIN dstr = string(day,format='(i2.2)') fname = thedir + 'ESACCI-L3S_SOILMOISTURE-SSMV-MERGED-'+ystr+mstr+dstr+'000000-fv00.1.nc' nc = ncdf_open(fname) IF yr EQ 1979 AND doy EQ 1 THEN BEGIN ncdf_varget,nc,'lat',lats ncdf_varget,nc,'lon',lons offset = [closest(lons,lon),closest(lats,lat)] ENDIF ncdf_varget1,nc,'sm',sm,offset=offset ncdf_close,nc IF sm LT 0 THEN sm = nan() ELSE sm *= 0.0001 smdata[doy-1l,yr-1979]=sm doy++ ENDFOR ENDFOR ENDFOR save,filename=indir+'wlef.sav',smdata END PRO wlef_readmodis ;products LST_day LST_night EVI, ;scale factor for daily LST is 0.02 (in Kelvin), fill is 0 ;valid 7500-65535 ;products EST_proc ;scale factor is 0.0001, fill is -3000 (valid -2000,10000) fname_lstd = '/sap/projects/psyn/modis/LST_day.csv' fname_lstn = '/sap/projects/psyn/modis/LST_night.csv' fname_evi = '/sap/projects/psyn/modis/EVI_proc.csv' lstd = (read_ascii(fname_lstd,delim=',',data_start=1)).(0) lstn = (read_ascii(fname_lstn,delim=',',data_start=1)).(0) evi = (read_ascii(fname_evi,delim=',',data_start=1)).(0) lstd_tm = lstd[0,*] + ((lstd[1,*]-1.0)/365.) lstn_tm = lstn[0,*] + ((lstn[1,*]-1.0)/365.) evi_tm = evi[0,*] + evi[1,*]/366. lstd = lstd[2:*,*] lstn = lstn[2:*,*] evi = evi[2:*,*] lstd[where(lstd LT 7500)]=nan() lstn[where(lstn LT 7500)]=nan() evi[where(evi LT -2000)]=nan() lstd*=0.02 lstn*=0.02 evi*=0.0001 lstd_avg = total(lstd,1,/nan)/total(finite(lstd),1) lstn_avg = total(lstn,1,/nan)/total(finite(lstn),1) evi_avg = total(evi,1,/nan)/total(finite(evi),1) lst_avg = (lstd_avg+lstn_avg)/2.0 tm_interp = replicate_arr(numgen(0,365,8)/365.,12) tm_interp += expand_arr(findgen(12)+2000.,46) wlef_lstd = reform(interpol(lstd_avg,lstd_tm,tm_interp),46,12) wlef_lstn = reform(interpol(lstn_avg,lstn_tm,tm_interp),46,12) wlef_lst = reform(interpol(lst_avg,lstd_tm,tm_interp),46,12) wlef_evi = reform(interpol(evi_avg,evi_tm,tm_interp),46,12) save,wlef_lstd,wlef_lstn,wlef_lst,wlef_evi,tm_interp,filename='/sap/projects/psyn/modis/wlef_modis.sav' stop END PRO wlef_readpsyn ;read fillenee for all years ;get pref par too ;gap fill ET ;store fillNEE, but flag by prefNEE indir = mydocs()+'cheas/wlef/filled/' indir2 = mydocs()+'cheas/wlef/flux/' yr = findgen(16)+1997 yr = yr[where(yr NE 2002)] nee = make_array(8784,15,/float,value=nan()) goodnee = make_array(8784,15,/byte,value=0b) et = nee par = nee temp = nee FOR i = 0,n_elements(Yr)-1 DO BEGIN yrstr = string(yr[i],format='(i4.4)') IF yr[i] LE 2005 THEN BEGIN fname = indir+'neefilled'+yrstr+'.txt' print,'Reading ',fname dta = (read_ascii(fname)).(0) bval = where(dta LE -999,nbd) IF nbd GT 0 THEN dta[bval] = nan() nee[0:n_elements(dta[16,*])-1,i]=dta[16,*] par[0:n_elements(dta[16,*])-1,i]=dta[12,*] temp[0:n_elements(dta[16,*])-1,i]=dta[7,*] goodnee[0:n_elements(dta[16,*])-1,i]=finite(dta[14,*]) AND (dta[16,*] EQ dta[14,*]) ;et! from flux file (timezone?) fname = indir2 + 'WLEF' + yrstr + '_flux.txt' print,'Reading ',fname dta = (read_ascii(fname)).(0) bval = where(dta LE -999,nbd) IF nbd GT 0 THEN dta[bval] = nan() et[0:n_elements(dta[0,*])-1,i] = dta[26,*] ENDIF ELSE BEGIN fname = indir+'prefnee_'+yrstr+'.txt' print,'Reading ',fname dta = (read_ascii(fname,data_Start=1)).(0) bval = where(dta LE -999,nbd) IF nbd GT 0 THEN dta[bval] = nan() nee[0:n_elements(dta[10,*])-1,i]=dta[10,*] et[0:n_elements(dta[10,*])-1,i]=dta[7,*] goodnee[0:n_elements(dta[10,*])-1,i]=finite(dta[6,*]) AND (dta[10,*] EQ dta[6,*]) par[0:n_elements(dta[16,*])-1,i]=dta[15,*] temp[0:n_elements(dta[16,*])-1,i]=dta[14,*] ENDELSE ENDFOR ;also precip and PAR - need to create pref precip/PAR from wcreek + lter save,filename=indir+'wlef19972012.sav',nee,goodnee,et,par,temp stop ;save the data ;processing (read save file, process) ;for each LT day, look for peak ER from previous or next night ;if >= 4 hours nighttime good, then take average (or max) - try both ;if >= 4 hours in daytime good, then take 10-2 filled average ;drawdown = daytime - nighttime ;anomaly analyis (separate function) - ;remove seasonal cycle replicate_arr(total(nee,2,/nan)/total(finite(nee),2),15) ;option 1: EMD, reconstruct minus first mode ;do same for ET (total in day) ;PAR, TEMP, other vars END PRO wlef_psyn_process,new=new indir = mydocs()+'cheas/wlef/filled/' IF keyword_set(new) THEN restore,indir+'wlef19972012-new.sav'ELSE restore,indir+'wlef19972012.sav' restore,mydocs()+'psyn/narr/narr_precip_wlef_daily.sav' restore,mydocs()+'psyn/modis/wlef_modis.sav' yr = findgen(16)+1997 yr = yr[where(yr NE 2002)] dn = intarr(8784,15) FOR i = 0,14 DO dn[*,i] = daynight(1,366,yr=yr[i],interval=24,lat=45.945,lon=-90.272) night = dn EQ 0 day = dn NE 0 noon = dn EQ 2 goodnee_night = average_Arr(night*goodnee,24,/tot) goodnee_day = average_Arr(day*goodnee,24,/tot) nightnee = nee nightnee[where(day)]=nan() nightnee[where(~goodnee)]=nan() daynee = nee daynee[where(night)]=nan() ; gpp = (366l*15l) neenight_mean = make_array(366,15,/float,value=nan()) neenight_max = neenight_mean needay_mean = neenight_mean FOR i = 0l,n_elements(neenight_mean)-1l DO IF goodnee_night[i] GE 4 THEN neenight_mean[i] = mean(nightnee[i*24l:(i*24l)+23l],/nan) FOR i = 0l,n_elements(neenight_max)-1l DO IF goodnee_night[i] GE 4 THEN neenight_max[i] = max(nightnee[i*24l:(i*24l)+23l],/nan) FOR i = 0l,n_elements(needay_mean)-1l DO IF goodnee_day[i] GE 4 THEN needay_mean[i] = mean(daynee[(i*24l)+10l:(i*24l)+13l],/nan) gpp = neenight_max-needay_mean gpp[where(gpp LT 0)]=nan() gpp[where(gpp GT 35)]=nan() meangpp = mean(gpp,dim=2,/nan) meangppr = replicate_arr(meangpp,15) sdgpp = stddev(gpp,dim=2,/nan) sdgppr = replicate_arr(sdgpp,15) gpp_an = gpp - meangppr gpp_rel = gpp_an/sdgppr ;normalized anomaly of daily PAR (daily tot) par_d = reform(average_arr(par,24,/nan),366,15) meanpar = mean(par_d,dim=2,/nan) meanparr = replicate_arr(meanpar,15) sdpar = stddev(par_d,dim=2,/nan) sdparr = replicate_arr(sdpar,15) par_an = par_d - meanparr par_rel = par_an/sdparr ;normalized anomaly of daily TEMP, Temp_max, temp_min, DTR temp_d = reform(average_arr(temp,24,/nan),366,15) meantemp = mean(temp_d,dim=2,/nan) meantempr = replicate_arr(meantemp,15) sdtemp = stddev(temp_d,dim=2,/nan) sdtempr = replicate_arr(sdtemp,15) temp_an = temp_d - meantempr temp_rel = temp_an/sdtempr mintemp_d = reform(average_arr(temp,24,/nan,/min),366,15) meanmintemp = mean(mintemp_d,dim=2,/nan) meanmintempr = replicate_arr(meanmintemp,15) sdmintemp = stddev(mintemp_d,dim=2,/nan) sdmintempr = replicate_arr(sdmintemp,15) mintemp_an = mintemp_d - meanmintempr mintemp_rel = mintemp_an/sdmintempr maxtemp_d = reform(average_arr(temp,24,/nan,/max),366,15) meanmaxtemp = mean(maxtemp_d,dim=2,/nan) meanmaxtempr = replicate_arr(meanmaxtemp,15) sdmaxtemp = stddev(maxtemp_d,dim=2,/nan) sdmaxtempr = replicate_arr(sdmaxtemp,15) maxtemp_an = maxtemp_d - meanmaxtempr maxtemp_rel = maxtemp_an/sdmaxtempr dtemp_d = maxtemp_d-mintemp_d meandtemp = mean(dtemp_d,dim=2,/nan) meandtempr = replicate_arr(meandtemp,15) sddtemp = stddev(dtemp_d,dim=2,/nan) sddtempr = replicate_arr(sddtemp,15) dtemp_an = dtemp_d - meandtempr dtemp_rel = dtemp_an/sddtempr ;normalized anomaly of daily total ET et_d = make_array(366,15,/float,value=nan()) goodet_d = average_arr(finite(et) AND (et GE 0),24,/tot) badet_d = where(~finite(et) OR (et LT 0) OR (et GT 700)) et_c = et et_c[badet_d] = nan() et_c[0:2000,5]=nan() FOR i = 0l,n_elements(et_d)-1l DO IF goodet_d[i] GT 12 THEN et_d[i] = mean(et_c[i*24l:(i*24l)+23l],/nan) meanet = mean(et_d,dim=2,/nan) meanetr = replicate_arr(meanet,15) sdet = stddev(et_d,dim=2,/nan) sdetr = replicate_arr(sdet,15) et_an = et_d - meanetr et_rel = et_an/sdetr ;normalized anomaly of daily SoilL meansoill = mean(soill,dim=2,/nan) meansoillr = replicate_arr(meansoill,15) sdsoill = stddev(soill,dim=2,/nan) sdsoillr = replicate_arr(sdsoill,15) soill_an = soill - meansoillr soill_rel = soill_an/sdsoillr ;normalized anomaly of daily PRECIP meanprecip = mean(precip,dim=2,/nan) meanprecipr = replicate_arr(meanprecip,15) sdprecip = stddev(precip,dim=2,/nan) sdprecipr = replicate_arr(sdprecip,15) precip_an = precip - meanprecipr precip_rel = precip_an/sdprecipr ;normalized anomaly of WUE wue = gpp/et_d meanwue = mean(wue,dim=2,/nan) meanwuer = replicate_arr(meanwue,15) sdwue = stddev(wue,dim=2,/nan) sdwuer = replicate_arr(sdwue,15) wue_an = wue - meanwuer wue_rel = wue_an/sdwuer ;normalized LST lst = make_Array(46,15,/float,value=nan()) wlef_lst[0:11]=nan() lst[*,3:4] = wlef_lst[*,0:1] lst[*,5:13] = wlef_lst[*,3:*] meanlst = mean(lst,dim=2,/nan) meanlstr = replicate_arr(meanlst,15) sdlst = stddev(lst,dim=2,/nan) sdlstr = replicate_arr(sdlst,15) lst_an = lst - meanlstr lst_rel = lst_an/sdlstr ;normalized EVI evi = make_Array(46,15,/float,value=nan()) wlef_evi[0:11]=nan() evi[*,3:4] = wlef_evi[*,0:1] evi[*,5:13] = wlef_evi[*,3:*] meanevi = mean(evi,dim=2,/nan) meanevir = replicate_arr(meanevi,15) sdevi = stddev(evi,dim=2,/nan) sdevir = replicate_arr(sdevi,15) evi_an = evi - meanevir evi_rel = evi_an/sdevir save,filename=mydocs()+'psyn/anomaly.sav',wue_rel,wue_an,wue,precip_rel,precip_an,precip,soill_rel,soill_an,soill,et_rel,et_an,et_d,$ dtemp_rel,dtemp_an,dtemp_d,mintemp_an,mintemp_rel,mintemp_d,maxtemp_an,maxtemp_rel,maxtemp_d,temp_an,temp_rel,temp_d,$ gpp_rel,gpp_an,gpp,neenight_mean,neenight_max,needay_mean,evi,evi_an,evi_rel,lst,lst_an,lst_rel stop ;for lag correl with relative anomaly - will have to start with _an, ;average ; recompute _rel END PRO wlef_psyn_emd ;do an emd of GPP ;do an emd of ET ;do a combined emd of GPP_an*ET_an ? ; restore,mydocs()+'psyn/anomaly.sav' restore,'/sap/projects/psyn/anomaly.sav' ; g = gpp[0:5310] ; e = et_d[0:5310] g = gpp[0:5489] e = et_d[0:5489] print,'Computing GPP emd' emd_gpp = demd(g,/nodel) print,'Computing ET emd' emd_et = demd(e,/nodel) save,emd_gpp,emd_et,filename=mydocs()+'psyn/emd.sav' hil_gpp = hilbert_spec(emd_gpp) hil_et = hilbert_spec(emd_et) han = hanning(8,8) hil_gpp2 = convol(hil_gpp,han)/16. hil_et2 = convol(hil_et,han)/16. m_gpp = n_elements(hil_gpp2[*,0])*total(hil_gpp2,2,/nan)/total(finite(hil_gpp2),2) m_wr = n_elements(hil_et2[*,0])*total(hil_et2,2,/nan)/total(finite(hil_et2),2) save,emd_gpp,emd_et,hil_gpp,hil_et,hil_gpp2,hil_et2,m_gpp,m_wr,filename=mydocs()+'psyn/emd.sav' ;hilbert spectrum analysis: ;one run all 5490 points - show marginal power spectrum (sum of total ; hilbert,2) ;emd of ;hanning filter = convol(array,hanning(8,8))/16. ;plot log of hilbert spectrum, freq ;co-plot of marginal power spectrum of GPP and ET END PRO wlef_psyn_emd_plot restore,mydocs()+'psyn/emd.sav' restore,mydocs()+'psyn/anomaly.sav' loadct,27 !p.color=255 !p.background=0 xr = (lindgen(5311)/360)+(findgen(5311) MOD 360)/360. !p.multi = [0,1,2] tvplot,transpose((hil_gpp2[0:400,*]))/7.31,/order,zrange=[0,1],title='P Power',ytitle='Wavelength (Days)',xtitle='Time (Years)',xrange=[0,14.75] tvplot,transpose((hil_et2[0:400,*]))/50.65,/order,zrange=[0,1],title='ET Power',ytitle='Wavelength (Days)',xtitle='Time (Years)',xrange=[0,14.75] ; tvplot,transpose((hil_gpp2[0:460,*])),/order,zrange=[0,1],title='P Power',ytitle='Wavelength (Days)',xtitle='Time (Years)',xrange=[0,14.75] ; tvplot,transpose((hil_et2[0:460,*])),/order,zrange=[0,1],title='ET Power',ytitle='Wavelength (Days)',xtitle='Time (Years)',xrange=[0,14.75] stop loadct,0 !p.color=0 !p.background=255 !p.multi = [0,1,4] plot,xr,gpp,ytitle='P (umol m-2 s-1)' plot,xr,et_d,ytitle='ET (W m-2)' stop loadct,0 !p.multi = [0,4,2,0,1] !p.color=0 !p.background=255 plot,m_gpp[0:400]/total(m_gpp),findgen(400),thick=2,title='P spectrum' plot,m_wr[0:400]/total(m_wr),findgen(400),thick=2,title='ET spectrum' stop END PRO wlef_psyn_lagtest restore,mydocs()+'psyn/anomaly.sav' sig = 0.1 lags = [0,1,3,8,15,30,60,90,180,360,720,1440] avg = [1,3,8,15,30,90,180,360,720,1440] glag = [1,3,8,15,30,60,90,180,360] ;randomly remove 6 days (first and last three) to get 360 data points gpp_an = gpp_an[3:362,*] et_an = et_an[3:362,*] wue_an = wue_an[3:362,*] precip_an = precip_an[3:362,*] soill_an = soill_an[3:362,*] dtemp_an = dtemp_an[3:362,*] mintemp_an = mintemp_an[3:362,*] maxtemp_an = maxtemp_an[3:362,*] temp_an = temp_an[3:362,*] evi_an2 = (congrid(evi_an[*,3:13],368,11,/center,/interp))[3:362,*] lst_an2 = (congrid(lst_an[*,3:13],368,11,/center,/interp))[3:362,*] gpp_an2 = gpp_an[*,3:13] lc_gpp = lagcorr(gpp_an,gpp_an,lags=lags,avg=avg,pval=pv_gpp,dof=df_gpp,/normalize,sig=sig) lc_et = lagcorr(et_an,gpp_an,lags=lags,avg=avg,pval=pv_et,dof=df_et,/normalize,/grange,glag=glag,cause=gr_et,sig=sig) lc_wue = lagcorr(wue_an,gpp_an,lags=lags,avg=avg,pval=pv_wue,dof=df_wue,/normalize,/grange,glag=glag,cause=gr_wue,sig=sig) lc_precip = lagcorr(precip_an,gpp_an,lags=lags,avg=avg,pval=pv_precip,dof=df_precip,/normalize,/grange,glag=glag,cause=gr_precip,sig=sig) lc_soill = lagcorr(soill_an,gpp_an,lags=lags,avg=avg,pval=pv_soill,dof=df_soill,/normalize,/grange,glag=glag,cause=gr_soill,sig=sig) lc_dtemp = lagcorr(dtemp_an,gpp_an,lags=lags,avg=avg,pval=pv_dtemp,dof=df_dtemp,/normalize,/grange,glag=glag,cause=gr_dtemp,sig=sig) lc_mintemp = lagcorr(mintemp_an,gpp_an,lags=lags,avg=avg,pval=pv_mintemp,dof=df_mintemp,/normalize,/grange,glag=glag,cause=gr_mintemp,sig=sig) lc_maxtemp = lagcorr(maxtemp_an,gpp_an,lags=lags,avg=avg,pval=pv_maxtemp,dof=df_maxtemp,/normalize,/grange,glag=glag,cause=gr_maxtemp,sig=sig) lc_temp = lagcorr(temp_an,gpp_an,lags=lags,avg=avg,pval=pv_temp,dof=df_temp,/normalize,/grange,glag=glag,cause=gr_temp,sig=sig) lc_evi = lagcorr(evi_an2,gpp_an2,lags=lags,avg=avg,pval=pv_evi,dof=df_evi,/normalize,/grange,glag=glag,cause=gr_evi,sig=sig) lc_lst = lagcorr(lst_an2,gpp_an2,lags=lags,avg=avg,pval=pv_lst,dof=df_lst,/normalize,/grange,glag=glag,cause=gr_lst,sig=sig) lc_gpp0 = lc_gpp[0,*] lc_et0 = lc_et[0,*] & lc_et0[where(pv_et[0,*] GT sig)] = 0 lc_wue0 = lc_wue[0,*] & lc_wue0[where(pv_wue[0,*] GT sig)] = 0 lc_precip0 = lc_precip[0,*] & lc_precip0[where(pv_precip[0,*] GT sig)] = 0 lc_soill0 = lc_soill[0,*] & lc_soill0[where(pv_soill[0,*] GT sig)] = 0 lc_dtemp0 = lc_dtemp[0,*] & lc_dtemp0[where(pv_dtemp[0,*] GT sig)] = 0 lc_mintemp0 = lc_mintemp[0,*] & lc_mintemp0[where(pv_mintemp[0,*] GT sig)] = 0 lc_maxtemp0 = lc_maxtemp[0,*] & lc_maxtemp0[where(pv_maxtemp[0,*] GT sig)] = 0 lc_temp0 = lc_temp[0,*] & lc_temp0[where(pv_temp[0,*] GT sig)] = 0 lc_evi0 = lc_evi[0,*] & lc_evi0[where(pv_evi[0,*] GT sig)] = 0 lc_lst0 = lc_lst[0,*] & lc_lst0[where(pv_lst[0,*] GT sig)] = 0 ;build into a 2D array lc0 = [lc_gpp0,lc_et0,lc_wue0,lc_precip0,lc_soill0,lc_temp0,lc_mintemp0,lc_maxtemp0,lc_dtemp0,lc_lst0,lc_evi0] lc_gpp2 = lc_gpp[1:*,*] & lc_gpp2[where(~finite(lc_gpp2) OR (pv_gpp[1:*,*] GT sig))]=0 lc_et2 = lc_et[1:*,*] & lc_et2[where(~finite(lc_et2) OR (pv_et[1:*,*] GT sig))]=0 lc_wue2 = lc_wue[1:*,*] & lc_wue2[where(~finite(lc_wue2) OR (pv_wue[1:*,*] GT sig))]=0 lc_precip2 = lc_precip[1:*,*] & lc_precip2[where(~finite(lc_precip2) OR (pv_precip[1:*,*] GT sig))]=0 lc_soill2 = lc_soill[1:*,*] & lc_soill2[where(~finite(lc_soill2) OR (pv_soill[1:*,*] GT sig))]=0 lc_dtemp2 = lc_dtemp[1:*,*] & lc_dtemp2[where(~finite(lc_dtemp2) OR (pv_dtemp[1:*,*] GT sig))]=0 lc_mintemp2 = lc_mintemp[1:*,*] & lc_mintemp2[where(~finite(lc_mintemp2) OR (pv_mintemp[1:*,*] GT sig))]=0 lc_maxtemp2 = lc_maxtemp[1:*,*] & lc_maxtemp2[where(~finite(lc_maxtemp2) OR (pv_maxtemp[1:*,*] GT sig))]=0 lc_temp2 = lc_temp[1:*,*] & lc_temp2[where(~finite(lc_temp2) OR (pv_temp[1:*,*] GT sig))]=0 lc_evi2 = lc_evi[1:*,*] & lc_evi2[where(~finite(lc_evi2) OR (pv_evi[1:*,*] GT sig))]=0 lc_lst2 = lc_lst[1:*,*] & lc_lst2[where(~finite(lc_lst2) OR (pv_lst[1:*,*] GT sig))]=0 ;build into 3d array lc2 = [[[lc_gpp2]],[[lc_et2]],[[lc_wue2]],[[lc_precip2]],[[lc_soill2]],[[lc_temp2]],[[lc_mintemp2]],[[lc_maxtemp2]],[[lc_dtemp2]],[[lc_lst2]],[[lc_evi2]]] lc_et3 = lc_et2 & lc_et3[where(abs(lc_gpp2) GT abs(lc_Et3))] = 0 lc_wue3 = lc_wue2 & lc_wue3[where(abs(lc_gpp2) GT abs(lc_Wue3))] = 0 lc_precip3 = lc_precip2 & lc_precip3[where(abs(lc_gpp2) GT abs(lc_Precip3))] = 0 lc_soill3 = lc_soill2 & lc_soill3[where(abs(lc_gpp2) GT abs(lc_Soill3))] = 0 lc_dtemp3 = lc_dtemp2 & lc_dtemp3[where(abs(lc_gpp2) GT abs(lc_Dtemp3))] = 0 lc_mintemp3 = lc_mintemp2 & lc_mintemp3[where(abs(lc_gpp2) GT abs(lc_Mintemp3))] = 0 lc_maxtemp3 = lc_maxtemp2 & lc_maxtemp3[where(abs(lc_gpp2) GT abs(lc_Maxtemp3))] = 0 lc_temp3 = lc_temp2 & lc_temp3[where(abs(lc_gpp2) GT abs(lc_Temp3))] = 0 lc_evi3 = lc_evi2 & lc_evi3[where(abs(lc_gpp2) GT abs(lc_Evi3))] = 0 lc_lst3 = lc_lst2 & lc_lst3[where(abs(lc_gpp2) GT abs(lc_Lst3))] = 0 lc3 = [[[abs(lc_gpp2)]],[[lc_et3]],[[lc_wue3]],[[lc_precip3]],[[lc_soill3]],[[lc_temp3]],[[lc_mintemp3]],[[lc_maxtemp3]],[[lc_dtemp3]],[[lc_lst3]],[[lc_evi3]]] save,/variables,filename=mydocs()+'psyn/correlations.sav',/compress stop ;figures LC0 as a tvplot ;tvplots of all LC2 (red/blue scale) ;diagnostics of LC3 END PRO wlef_psyn_figures,new=new indir = mydocs()+'cheas/wlef/filled/' ; restore,indir+'wlef19972012.sav' IF keyword_set(new) THEN restore,indir+'wlef19972012-new.sav'ELSE restore,indir+'wlef19972012.sav' restore,mydocs()+'psyn/correlations.sav' ;Figure 0 ;Conceptual figure of leaf and forest ;Figure 1 ;Hourly NEE and ET (two panels) ;NEE and ET fig1tm = 1997.0D + (dindgen(8784l*16l)/8784.0D) fig1tm = fig1tm[where(long(fig1tm) NE 2002)] et[where(et LT 0)]=nan() et[where(et GT 700)]=nan() et[where(fig1tm GT 2003.0 AND fig1tm LT 2003.1)]=nan() !p.multi = [0,1,2] loadct,0 red = fsc_color('red',200) plot,fig1tm,nee,psym=1,yrange=[-40,40],xticks=16,xtickformat='(i4.4)',xrange=[1997,2013],ytitle='umol m-2 s-1' oplot,fig1tm,smooth(nee[*],366,/nan,/edge),thick=3,color=200 plot,fig1tm,et,psym=1,yrange=[0,700],xticks=16,xtickformat='(i4.4)',xrange=[1997,2013],ytitle='W m-2' oplot,fig1tm,smooth(et[*],366,/nan,/edge),thick=3,color=200 stop ;NEW FIG X - also daily all others ;TEMP_D, evi, lst, precip (max 40) , -1, soill, !p.multi = [0,1,5] figAtm = 1997.0D + (dindgen(366l*16l)/366.0D) figAtm = figAtm[where(long(figAtm) NE 2002)] figBtm = 1997.0D + (dindgen(46l*16l)/46.0D) figBtm = figBtm[where(long(figBtm) NE 2002)] plot,figAtm,temp_d,psym=1,xticks=16,xtickformat='(i4.4)',xrange=[1997,2013],ytitle='degrees C',charsize=2 plot,figAtm,precip,max=50,psym=-1,xticks=16,xtickformat='(i4.4)',xrange=[1997,2013],ytitle='mm',charsize=2 plot,figAtm,soill,psym=1,xticks=16,xtickformat='(i4.4)',xrange=[1997,2013],ytitle='percent',charsize=2 plot,figBtm,evi,psym=1,xticks=16,xtickformat='(i4.4)',xrange=[1997,2013],ytitle='unitless',charsize=2 plot,figBtm,lst,psym=1,xticks=16,xtickformat='(i4.4)',xrange=[1997,2013],ytitle='degrees C',charsize=2 ;Figure 2 ;PD time series (360 day), anomaly time series, ;weekly average, relative anomaly fig2tm = 1997.0D + (dindgen(360l*16l)/360.0D) fig2tm = fig2tm[where(long(fig2tm) NE 2002)] sx = make_array(n_elements(gpp_an),/float,value=nan()) av = 30 FOR k = 0l,n_elements(gpp_an)-av DO sx[k] = mean(gpp_an[k:(av+k-1)],/nan) ndy = 360l nyr = 15l mn = mean(reform(sx,ndy,nyr),dim=2,/nan) mnr = replicate_arr(mn,nyr) sd = stddev(reform(sx,ndy,nyr),dim=2,/nan) sdr = replicate_arr(sd,nyr) sx_an = (sx-mnr) sx = (sx-mnr) / sdr !p.multi = [0,1,3] red = fsc_color('red',200) plot,fig2tm,gpp[3:362,*],xticks=16,xtickformat='(i4.4)',xrange=[1997,2013],ytitle='umol m-2 s-1',yrange=[0,35] plot,fig2tm,gpp_an,xticks=16,xtickformat='(i4.4)',xrange=[1997,2013],ytitle='umol m-2 s-1',yrange=[-15,15] oplot,fig2tm,sx_an,color=200,thick=3 oplot,[1997,2013],[0,0] plot,fig2tm,sx,yrange=[-3,3],xticks=16,xtickformat='(i4.4)',xrange=[1997,2013],ytitle='unitless' oplot,[1997,2013],[0,0] stop ;Figure 3 ;DOF plot df = df_gpp[0,*] !p.multi=[0,2,2] plot,avg,df,/xlog,xtitle='Averaging period (days)',ytitle='DOF (days)',thick=2,title='DOF GPP (n=5490)' stop ;Figure 4 ;Lag 0 show - all vars (lc0) ;[lc_gpp0,lc_et0,lc_wue0,lc_precip0,lc_soill0,lc_temp0,lc_mintemp0,lc_maxtemp0,lc_dtemp0,lc_lst0,lc_evi0] ;order: EVI, ET, WUE, PRECIP, SOILL, TEMPs lc0[9:10,0:1] = 0.0 ;remove 1 and 3 day avg from MODIS !p.multi = [0,2,2] loadct,33 tvlct,rd,gr,bl,/get rd[127] = 255 gr[127] = 255 bl[127] = 255 tvlct,rd,gr,bl !p.background = 127 tvplot,transpose(lc0[[10,1,2,3,4,5,6,7,8,9],*]),zrange=[-0.6,0.6],xticks=10,xtickv=numgen(0.0,9.0,0.9)+0.45,xtickn=['1','3','8','15','30','90','180','360','720','1440'],xtitle='Averaging Period (Days)',yticks=10,ytickv=numgen(0.0,9.0,0.9)+0.45,ytickn=['LST','Trange','Tmax','Tmin','Tavg','Qsoil','Precip','WUE','ET','EVI'],ntext=9 stop ;[lc_gpp0,lc_et0,lc_wue0,lc_precip0,lc_soill0,lc_temp0,lc_mintemp0,lc_maxtemp0,lc_dtemp0,lc_lst0,lc_evi0] ;transpose(lc0) - colors - no negative values - so just show white to ; black ,/creverse ;Figure 5 ;show LC2 for GPP, EVI, Tmean, LST ;cols are lags[1:*], rows are avg, z is variable lc2[*,0:1,9:10] = 0.0 !p.multi = [0,2,2] loadct,33 tvlct,rd,gr,bl,/get rd[127] = 255 gr[127] = 255 bl[127] = 255 tvlct,rd,gr,bl !p.background = 127 tvplot,transpose(lc2[*,*,0]),zrange=[-1,1],xticks=10,xtickv=numgen(0.0,9.0,0.9)+0.45,xtickn=['1','3','8','15','30','90','180','360','720','1440'],xtitle='Averaging Period (Days)',yticks=10,ytickv=numgen(0.0,9.0,0.9)+0.45,ytickn=['1440','720','360','180','90','60','30','15','8','3','1'],ntext=9,ytitle='Lag (days)',title='a) Pd' tvplot,transpose(lc2[*,*,10]),zrange=[-1,1],xticks=10,xtickv=numgen(0.0,9.0,0.9)+0.45,xtickn=['1','3','8','15','30','90','180','360','720','1440'],xtitle='Averaging Period (Days)',yticks=10,ytickv=numgen(0.0,9.0,0.9)+0.45,ytickn=['1440','720','360','180','90','60','30','15','8','3','1'],ntext=9,ytitle='Lag (days)',title='b) EVI' FOR x = 0,9 DO FOR y = 0,10 DO IF (lc2[y,x,10] EQ lc3[y,x,10]) AND (lc2[y,x,10] NE 0.0) THEN box,x*0.9,10-(y*0.9),(x*0.9)+0.9,9.1-(y*0.9),/nofill,/outline,outcolor=255,outthick=2 tvplot,transpose(lc2[*,*,5]),zrange=[-1,1],xticks=10,xtickv=numgen(0.0,9.0,0.9)+0.45,xtickn=['1','3','8','15','30','90','180','360','720','1440'],xtitle='Averaging Period (Days)',yticks=10,ytickv=numgen(0.0,9.0,0.9)+0.45,ytickn=['1440','720','360','180','90','60','30','15','8','3','1'],ntext=9,ytitle='Lag (days)',title='c) Tmean' FOR x = 0,9 DO FOR y = 0,10 DO IF (lc2[y,x,5] EQ lc3[y,x,5]) AND (lc2[y,x,5] NE 0.0) THEN box,x*0.9,10-(y*0.9),(x*0.9)+0.9,9.1-(y*0.9),/nofill,/outline,outcolor=255,outthick=2 tvplot,transpose(lc2[*,*,9]),zrange=[-1,1],xticks=10,xtickv=numgen(0.0,9.0,0.9)+0.45,xtickn=['1','3','8','15','30','90','180','360','720','1440'],xtitle='Averaging Period (Days)',yticks=10,ytickv=numgen(0.0,9.0,0.9)+0.45,ytickn=['1440','720','360','180','90','60','30','15','8','3','1'],ntext=9,ytitle='Lag (days)',title='d) LST' FOR x = 0,9 DO FOR y = 0,10 DO IF (lc2[y,x,9] EQ lc3[y,x,9]) AND (lc2[y,x,9] NE 0.0) THEN box,x*0.9,10-(y*0.9),(x*0.9)+0.9,9.1-(y*0.9),/nofill,/outline,outcolor=255,outthick=2 stop ;Figure 6 ;show LC3 for ET(1), WUE(2), Precip(3), SoilL(4) !p.multi = [0,2,2] loadct,33 tvlct,rd,gr,bl,/get rd[127] = 255 gr[127] = 255 bl[127] = 255 tvlct,rd,gr,bl !p.background = 127 tvplot,transpose(lc2[*,*,1]),zrange=[-1,1],xticks=10,xtickv=numgen(0.0,9.0,0.9)+0.45,xtickn=['1','3','8','15','30','90','180','360','720','1440'],xtitle='Averaging Period (Days)',yticks=10,ytickv=numgen(0.0,9.0,0.9)+0.45,ytickn=['1440','720','360','180','90','60','30','15','8','3','1'],ntext=9,ytitle='Lag (days)',title='ET' FOR x = 0,9 DO FOR y = 0,10 DO IF (lc2[y,x,1] EQ lc3[y,x,1]) AND (lc2[y,x,1] NE 0.0) THEN box,x*0.9,10-(y*0.9),(x*0.9)+0.9,9.1-(y*0.9),/nofill,/outline,outcolor=255,outthick=2 tvplot,transpose(lc2[*,*,2]),zrange=[-1,1],xticks=10,xtickv=numgen(0.0,9.0,0.9)+0.45,xtickn=['1','3','8','15','30','90','180','360','720','1440'],xtitle='Averaging Period (Days)',yticks=10,ytickv=numgen(0.0,9.0,0.9)+0.45,ytickn=['1440','720','360','180','90','60','30','15','8','3','1'],ntext=9,ytitle='Lag (days)',title='b) WUE' FOR x = 0,9 DO FOR y = 0,10 DO IF (lc2[y,x,2] EQ lc3[y,x,2]) AND (lc2[y,x,2] NE 0.0) THEN box,x*0.9,10-(y*0.9),(x*0.9)+0.9,9.1-(y*0.9),/nofill,/outline,outcolor=255,outthick=2 tvplot,transpose(lc2[*,*,3]),zrange=[-1,1],xticks=10,xtickv=numgen(0.0,9.0,0.9)+0.45,xtickn=['1','3','8','15','30','90','180','360','720','1440'],xtitle='Averaging Period (Days)',yticks=10,ytickv=numgen(0.0,9.0,0.9)+0.45,ytickn=['1440','720','360','180','90','60','30','15','8','3','1'],ntext=9,ytitle='Lag (days)',title='c) Precip' FOR x = 0,9 DO FOR y = 0,10 DO IF (lc2[y,x,3] EQ lc3[y,x,3]) AND (lc2[y,x,3] NE 0.0) THEN box,x*0.9,10-(y*0.9),(x*0.9)+0.9,9.1-(y*0.9),/nofill,/outline,outcolor=255,outthick=2 tvplot,transpose(lc2[*,*,4]),zrange=[-1,1],xticks=10,xtickv=numgen(0.0,9.0,0.9)+0.45,xtickn=['1','3','8','15','30','90','180','360','720','1440'],xtitle='Averaging Period (Days)',yticks=10,ytickv=numgen(0.0,9.0,0.9)+0.45,ytickn=['1440','720','360','180','90','60','30','15','8','3','1'],ntext=9,ytitle='Lag (days)',title='d) Qsoil' FOR x = 0,9 DO FOR y = 0,10 DO IF (lc2[y,x,4] EQ lc3[y,x,4]) AND (lc2[y,x,4] NE 0.0) THEN box,x*0.9,10-(y*0.9),(x*0.9)+0.9,9.1-(y*0.9),/nofill,/outline,outcolor=255,outthick=2 stop ;Panels for lc2 ;draw squares around lc3s? ;show only for GPP, EVI, ET, WUE, Precip, SOIll, meanTEMP, LST, (2x4) abcdefgh ;Figure 7 EMD stuff (hold) ;Table 1 ;Variables ;Table 2 ;Granger Cause ;in excel END ;steps ;1. Pref par for 2011 and 2012 (fix 2012 manually?) - got 2011, fix - DONE ; 2012 man (get ; ratio 2011) ;2. fix/get precip - DONE ;3. anomaly code (see notes above) - DONE ; need anomaly of MODIS wlef_lst and wlef_evi (46,12) 8 day - DONE ;4. lag/correl code ; when averaging - need to padd to 368 for 4,8,16, then to 384 for 32, ; 64 128. beyond 128 - need to account for ; interannual ;DONE - lagcorr.pro and granger.pro ;NOW: run through all explanatory variables for gpp ;restandardize? ; take two signals ; decide whether to detrend or not ; 2D output - correlation, signifiance, autocorrelation, reduced r^2 (DOF) ; smoother 0, 2, 4, 8, 16, ... until ; lags 0, -1, -2, -4, -8, -16, -32, -64, .... (shift/chop (make nan)) ; ;5. EMD/hilbert code ;1. gpp_an = EMD? - improve computation ;2. cross spectra EMD of most promising correlation (et?) ;lots of fun! ;LAG analysis plans: ;randomly remove 6 days (first and last three) to get 360 data points ;compute averages at 1,3,8,15,30,90,180,360,720,1440 ;compate lags at same 0,1,3,8,15,30,60,90,180,360,720,1440 ;take anomalies - average, then standardize, then lag (lagcorr ; normalize) 3:362 ;do for: ;gpp_an vs. ;gpp_an ;wue_an, precip_an, soill_an, et_an, dtemp_an, mintemp_an, maxtemp_an, ;temp_an, evi_an, lst_an (latter two require subset of gpp_an) ;take lcorr, remove where pval > 0.05 (set to 0, set nan to 0) ;compare gpp_an to others ;granger across each ;plot zero lags - for direct correlation by averaging ;then tvplot all lags by all averages, zero out ;REDO ;archive vars ;WLEF process with new data (.sav file) ;nee (fillednee) ;temp - air temp ;par ;et ;precip (daily) - compare to NARR ;redo anomaly analyis ;redo EMD analysis ;replot / retable PRO wlef_psyn_redo restore,'/data/incoming/wlef/alldata19952012-processed.sav' nee_new = fillnee[*,[2,3,4,5,6,8,9,10,11,12,13,14,15,16,17]] nn = replicate_arr(total(nee_new,2,/nan)/total(finite(nee_new),2),15) nee = merge_Array(nee_new,nn) et = q_pref[*,[2,3,4,5,6,8,9,10,11,12,13,14,15,16,17]] temp = tair_30_filled[*,[2,3,4,5,6,8,9,10,11,12,13,14,15,16,17]] par = par_pref[*,[2,3,4,5,6,8,9,10,11,12,13,14,15,16,17]] cp = c_prefflag[*,[2,3,4,5,6,8,9,10,11,12,13,14,15,16,17]] goodnee = (cp GE 0) ; precip_new = average_arr(precip[*,[2,3,4,5,6,8,9,10,11,12,13,14,15,16,17]],24,/tot,/nan) indir = mydocs()+'cheas/wlef/filled/' save,filename=indir+'wlef19972012-new.sav',nee,et,temp,par,goodnee ; restore,indir+'wlef19972012.sav' ; restore,mydocs()+'psyn/narr/narr_precip_wlef_daily.sav' ; restore,mydocs()+'psyn/modis/wlef_modis.sav' ;compare ;data has 2002 removed ;save the four vars stop END ;modify and run wlef_psyn_process ;modify and run wlef_psyn_emd, wlef_psyn_emd_plot, wlef_psyn_lagtest, wlef_psyn_figures ;Review: hepnology, PRO wlef_psyn_gsl ;for reviewer 1 indir = mydocs()+'cheas/wlef/filled/' restore,indir+'wlef19972012-new.sav' restore,mydocs()+'psyn/correlations.sav' gth = 5 ;look at gpp - for each yar gpp2 = fltarr(368,15) FOR i = 0,14 DO gpp2[*,i] = ts_smooth(congrid(zapbadval(average_arr(gpp[*,i],8,/nan)),368,/center,/interp),30) gpp2[266:*,7] = nan() gpp2[266,7] = 5 gpp2[267:*,7] = 0 cup = fltarr(15) FOR i = 0,14 DO cup[i] = n_elements(where(gpp2[*,i] GT gth)) evi2 = fltarr(368,15) FOR i = 0,14 DO evi2[*,i] = congrid(evi[*,i],368,/center,/interp) sday = fltarr(15) FOR i = 0,14 DO sday[i] = (where(gpp2[*,i] GT gth))[0] eday = fltarr(15) FOR i = 0,14 DO eday[i] = 368-(where(reverse(gpp2[*,i]) GT gth))[0] ev_av = total(evi2,1) ev_c = fltarr(15) FOR i = 0,14 DO ev_c[i] = mean(evi2[sday[i]:eday[i],i],/nan) gppe = fltarr(46,15) FOR i = 0,14 DO gppe[*,i] = [average_Arr(gpp[*,i],8,/nan),mean(gpp[360:*,i],/nan)] ev_sday = fltarr(15) ev_eday = fltarr(15) de = deriv(evi2) de[0:100,*] = nan() de[300:*,*] = nan() FOR i = 0,14 DO ev_sday[i] = (where(de[*,i] EQ max(de[*,i],/nan)))[0] FOR i = 0,14 DO ev_eday[i] = (where(de[*,i] EQ min(de[*,i],/nan)))[0] ev_sday[where(ev_sday LE 101)]=nan() ev_Eday[where(ev_eday GE 299)]=nan() ev_Eday[where(ev_eday LE 101)]=nan() evi3 = evi2 evi3[0:100,*] = nan() evi3[305:*,*]=nan() eth = 0.31 cup_E = cup FOR i = 0,14 DO cup_e[i] = n_elements(where(evi3[*,i] GT eth)) esday = fltarr(15) FOR i = 0,14 DO esday[i] = (where(evi3[*,i] GT eth))[0] eeday = fltarr(15) FOR i = 0,14 DO eeday[i] = 368-(where(reverse(evi3[*,i]) GT eth))[0] eeday[where(eeday GT 305)]=nan() esday[where(esday LT 100)]=nan() cup_e[where(cup_E LT 10)]=nan() ; ;CORRELATIONS: ;CUP to EV_C: long growing seasons have low EVI , but effect is small ;0.48 to 0.44, correlation is -0.88 ;strongest is to end of CUP period, -0.92 ;make like EVI - what is direct corr? it is 0.87 - explains 75% of ; biweekly variation ;correlation disappears with averaging, spurious ;cite HEinsch paper ;EVI may ben a detector of when leaves come out, but tests here show ;little relationship, -0.12 for end of season, -0.014 for start of ;season, -0.29 for carbon uptake period ;EVI during carbon uptake has a negative correlation ;doesn't invalidate it, just shows complexity of evaluation at site. END