# -*- coding: utf-8 -*- """ Created on Tue Apr 26 22:11:25 2016 @author: wxf The code was wrote by Xufeng Wang """ #comprison of phenology methods """ methods choosed: 1. order polynomial(order=6) include 4 steps: step 1: get multi-year average seasonal NDVI dynamic curve. step 2: get NDVI_ratio N_ratio(t)=[NDVI(t+1)-NDVI(t)]/NDVI(t). step 3: get the maximum NDVI_ratio and minimum NDVI_ratio time. using the corresponding NDVI value at that time as NDVI threshhold value to get start of greening and dormancy. step 4: curve fitting using the NDVI and Julian day, function: NDVI = a0 + a1*x + a2*x**2 + a3*x**3 ... + an*x**n, n = 6 January to September data used to fit the start of greening, July to December data used to fit the start of dormancy paper: Zhang G.L., Zhang Y.J., Dong J.W., Xiao X.M., Green-up date in the Tibetan Plateau have continuously advanced from 1982 to 2011, PNAS, 2013, 110(11): 4309-4314. 2. logistic function step1: divide a full year of vegetation index data into appropriate sections for a spring sigmoid and a fall sigmoid. step2: fit sigmoid function, y = c / (1 + exp(b*(x+a/b)))+d; t is time(doy), a and b are fitting prameters, c+d is the maximum VI value, d is the initial background VI value. step3: get extrema of first derived value of function in step1 Zhang X.Y., Friedl M.A., Schaaf C.B., et al., Monitoring vegetation phenology using MODIS. Remote Sensing of Environment 84 (2003) 471–475 """ import numpy as np import numpy.polynomial.polynomial as poly def dbl_logistic_model ( p, doys ): """Double logistic function, as in Zhang et al or Sobrino and Juliean""" y= p[0] + p[1]* ( 1./(1+np.exp(p[2]*(doys-p[3]))) + 1./(1+np.exp(-p[4]*(doys-p[5]))) - 1) return y def mismatch_function ( p, pheno_func, ndvi, days): '''This function is used to measure the mismatch between two datasets ''' # output stores the predictions out = [] # define funciton using lambda structure fitness = lambda pa, ndvi_in, days: ndvi_in - pheno_func ( pa, days ) oot = fitness ( p, ndvi, days ) [ out.append ( x ) for x in oot ] return np.array(out).squeeze() def fit_phenology_model_dbl ( ndvi_1year, t_axis, xinit=None ): """Fitting NDVI data with double logistic model""" from scipy.optimize import leastsq n_params = 6 # 6 terms if xinit is None: xinit = [.5,] * n_params xinit[0] = ndvi_all.min() xinit[1] = ndvi_all.max() - ndvi_all.min() xinit[2] = 0.19 xinit[3] = 120 xinit[4] = 0.13 xinit[5] = 260 ( xsol, msg ) = leastsq ( mismatch_function, xinit, \ args=( dbl_logistic_model, ndvi_1year, t_axis ), maxfev=100000 ) #print xsol,msg ax = dbl_logistic_model ( xsol, np.arange(1, 366)) return ( xsol, msg,ax ) def get_phen_date_model_dbl(daily_ndvi,days): import numpy as np try: daily_ndvi = daily_ndvi.tolist() max_ndvi= max(daily_ndvi) ind_max = daily_ndvi.index(max_ndvi) sep_day = days[ind_max] days = np.asarray(days) daily_ndvi = np.asarray(daily_ndvi) first_half_ind = days < sep_day #spr_days = days[first_half_ind] spr_ndvi = daily_ndvi[first_half_ind] second_half_ind = days > sep_day #fall_days = days[second_half_ind] fall_ndvi = daily_ndvi[second_half_ind] spr_first_der = np.diff(spr_ndvi) doy_onset = np.diff(spr_first_der).argmax() # I get 297 for this # Onset of greenness is the minimum of the derivative fall_first_der = np.diff(fall_ndvi) doy_dormacy = np.diff(fall_first_der).argmax() + sep_day # I get 109 for my example return (doy_onset,doy_dormacy) except: return (-9999,-9999) def get_onset_dormacy_ndvi(NDVI_data,ndata_year): nyear = len(NDVI_data)/ndata_year NDVI_data = np.reshape(NDVI_data,[nyear,ndata_year]) years_avg_ndvi = np.nanmean(NDVI_data,axis=0) # get the mean of each column ndvi_ratio = [] for i in range(ndata_year-1): ndvi_ratio.append((years_avg_ndvi[i+1]-years_avg_ndvi[i])/years_avg_ndvi[i]) max_ratio = np.nanmax(ndvi_ratio) min_ratio = np.nanmin(ndvi_ratio) ind_onset = ndvi_ratio.index(max_ratio) ind_dormacy = ndvi_ratio.index(min_ratio) + 1 onset_ndvi = years_avg_ndvi[ind_onset] dormacy_ndvi = years_avg_ndvi[ind_dormacy] return [onset_ndvi, dormacy_ndvi] def fit_phenology_model_poly (ndvi_pixel, doys,onset_ndvi,dormacy_ndvi): SOS = -9999 EOS = -9999 #divide the data into two parts # spring part data # fitting the curve doys = np.asarray(doys) ind_spr = doys < 256 spr_doy=doys[ind_spr] spr_ndvi=ndvi_pixel[ind_spr] p_spr=poly.polyfit(spr_doy,spr_ndvi,6) # get the phenology date first_half_days=np.arange(1,210) ndvi_first_half=poly.polyval(first_half_days,p_spr) ndvi_first_half = ndvi_first_half.tolist() max_ind = ndvi_first_half.index(np.nanmax(ndvi_first_half)) max_ndvi_day_sec_first = first_half_days[max_ind] for d in range(max_ndvi_day_sec_first,50,-1): ind_d = d -1 -1 if(ndvi_first_half[ind_d]<=onset_ndvi): SOS=d break # fall part data # fitting the curve ind_fall = doys > 180 fall_doy=doys[ind_fall] fall_ndvi=ndvi_pixel[ind_fall] p_fall=poly.polyfit(fall_doy,fall_ndvi,6) # get the phenology date second_half_days=np.arange(181,365) ndvi_second_half=poly.polyval(second_half_days,p_fall) ndvi_second_half = ndvi_second_half.tolist() max_ind = ndvi_second_half.index(np.nanmax(ndvi_second_half)) max_ndvi_day_second = second_half_days[max_ind] for d in range(max_ndvi_day_second,340): ind_d = d-181 if(ndvi_second_half[ind_d]<=dormacy_ndvi): EOS=d break return (SOS,EOS,first_half_days,ndvi_first_half,second_half_days,ndvi_second_half) # the main program ########################## import pandas as pd import numpy as np import matplotlib.pyplot as plt import savitzky_golay #### setting location for different run ########################### #set the input path' infile_gimms_new_v1 = 'test_data/GIMMS_new_v1_NDVI_flux_sites.csv' infile_gimms_old = 'test_data/GIMMS_old_NDVI_flux_sites.csv' infile_modis = 'test_data/MODIS_NDVI_flux_sites.csv' infile_VGT = 'test_data/SPOT_VEG_NDVI_flux_sites.csv' infile_SWF = 'test_data/SeaWiFs_NDVI_flux_sites.csv' outpath='D:/UNH_visiting/Phen_QT' sites_name = ['Arou', 'Ha2', 'HaM', 'Dan'] #### start setting here parameters ######################### dataset = 'SPOT_VEG' col = 0 # set index of sites in sites_name site = sites_name[col] start_year = 1999 end_year = 2013 years = range(start_year,end_year+1) ndata_year = 36 doy = range(5,365,10) df=pd.DataFrame.from_csv(infile_VGT,header=0) print df.columns.values #print df.index #### end setting parameters ######################### NDVI_data = df.as_matrix(columns=[site]) #pheno_model="dbl_polynomial" SOS_sg=[] EOS_sg=[] SOS_sg_poly=[] EOS_sg_poly=[] para_phen = [] [onset_ndvi,dormacy_ndvi] = get_onset_dormacy_ndvi(NDVI_data,ndata_year) NDVI_data = np.squeeze(NDVI_data) NDVI_data_sg=savitzky_golay.savitzky_golay(NDVI_data,11,4) [onset_ndvi_sg,dormacy_ndvi_sg] = get_onset_dormacy_ndvi(NDVI_data_sg,ndata_year) print onset_ndvi,dormacy_ndvi print onset_ndvi_sg,dormacy_ndvi_sg for yr in years: t_axis=doy start_ind=(yr-start_year)*ndata_year end_ind=(yr-start_year+1)*ndata_year ndvi_all=NDVI_data[start_ind:end_ind] ### choose the NDVI value column NDVI_sg=NDVI_data_sg[start_ind:end_ind] # using polynomial function to fit [sos_poly_sg, eos_poly_sg,spr_days,spr_ndvi_sg,fall_days,fall_ndvi_sg] = \ fit_phenology_model_poly (NDVI_sg, t_axis, onset_ndvi_sg, dormacy_ndvi_sg) SOS_sg_poly.append(sos_poly_sg) EOS_sg_poly.append(eos_poly_sg) # using sig function to fit xinit=None pheno_model="dbl_logistic" #### time reconstructed NDVI data [para, msg, result_sg]=fit_phenology_model_dbl ( NDVI_sg, t_axis, xinit ) para_phen.append(para[3]) para_phen.append(para[5]) doys = range(1,len(result_sg)+1) [s_sos_sg, s_eos_sg] = get_phen_date_model_dbl(result_sg,doys) SOS_sg.append(s_sos_sg) EOS_sg.append(s_eos_sg) plt.plot(t_axis,ndvi_all,'bo',label = 'NDVI') plt.plot(doys,result_sg,'g-',label = 'Log_fit') plt.plot(spr_days,spr_ndvi_sg,'k-',label = 'Poly_fit') plt.plot(fall_days,fall_ndvi_sg,'k-') plt.plot([s_sos_sg,s_sos_sg],[-0.002,0.98],'g--',label = 'Log_phen') plt.plot([s_eos_sg,s_eos_sg],[-0.002,0.98],'g--') plt.plot([sos_poly_sg,sos_poly_sg],[-0.002,0.98],'k--',label = 'Poly_phen') plt.plot([eos_poly_sg,eos_poly_sg],[-0.002,0.98],'k--') plt.xlim([1,365]) plt.grid(True) ''' plt.legend(['ndvi','ndvi_sg','Log_fit','Log_fit_sg','sos','eos','sos_sg','eos_sg', \ 'spr_poly','fall_poly','spr_poly_sg','fall_poly_sg','sos_poly', \ 'eos_poly','sos_sg_poly','eos_sg_poly','sos_slog','eos_slog',\ 'sos_slog_sg','eos_slog_sg'], bbox_to_anchor=(1.4, 1.2)) ''' plt.legend(bbox_to_anchor=(1.4, 1.2)) title='%d' %(yr) plt.title(title) plt.show() plt.plot(years,SOS_sg,'b-o') plt.plot(years,SOS_sg_poly,'g-o') plt.legend(["%s_SOS_sig_sg" %(sites[col]), \ "%s_SOS_poly_sg" %(sites[col])], \ bbox_to_anchor=(0.8, 1.0),loc='best') plt.grid() plt.ylim([50,200]) plt.show() plt.plot(years,EOS_sg,'b-o') plt.plot(years,EOS_sg_poly,'g-o') plt.legend(["%s_EOS_sig_sg" %(sites[col]), \ "%s_EOS_poly_sg" %(sites[col])], \ bbox_to_anchor=(0.8, 1.0),loc='best') plt.grid() plt.ylim([220,360]) plt.show() outdata_poly = np.vstack((SOS_sg,EOS_sg,SOS_sg_poly,EOS_sg_poly)) outdata_sig = np.vstack((SOS_sg,EOS_sg,SOS_sg_poly,EOS_sg_poly)) outdata_poly = np.transpose(outdata_poly) outdata_sig = np.transpose(outdata_sig) outfile = 'D:/UNH_visiting/Phen_QT/papers/JGR-B/Codes and data share/Python code for retrieving SOS from NDVI/phendate_%s_%s.csv' %(dataset,site) indx = years fld = ['SOS_sg_log','EOS_sg_log','SOS_sg_poly','EOS_sg_poly'] df_out=pd.DataFrame(data=outdata_poly,index = indx, columns=fld) df_out.to_csv(outfile)