# INPUTS

1. sitename (text_string)

2. year (e.g. 2013)

3. doy (e.g. doy=2 is January 2)

4. latitude (decimal degrees)

5. longitude (decimal degrees)

6. tmin = minumum daily air temperature (degrees C)

7. tmax = maximum daily air temperature (degrees C)

8. sr_min_ndvi = 0.1

9. sr_min_evi = 0.05

10. sr_min_red = 0

11. sr_min_nir = 0

12. sr_min_mir = 0

13. sr_min_blue = 0

14. slope_thrshld = 0.05   #maximum change per day for acceptable quality  

15. nsmooth = 20   #number of smoothing iterations to perform

16. tc = 7   #timeconstant for running means used for smoothing environmental variables (days)

17. dt = 1   #time interval between observations (days)


# Gap-filling and smoothing

#Reliable high-quality MODIS observations of surface reflectance, albedo or vegetation indices are not available every day due to clouds and other atmospheric effects. The first step is to filter out the bad data points, the second step is to fill the gaps between the remaining good points using linear interpolation and the third step is to smooth the data to create a physically realistic time-series of daily values. A generic function could be created to apply the same processing steps to surface reflectance, albedo or vegetation indices. This gap-filling and smoothing approach would not appropriate, however, for thermal observations because land surface temperatures can be different under clear and cloudy conditions.

# Filtering

#For sr=MYD13Q1_250m_16_days_NIR_reflectance, MYD13Q1_250m_16_days_red_reflectance, MYD13Q1_250m_16_days_NDVI, etc.

#sr_min=minimum value of reasonable surface reflectance (0.1 for NDVI, 0.05 for EVI, ? for red reflectance, ? for NIR reflectance, etc.). Reasonable values for this is half of the global bare soil minimum for the respective surface reflectance (e.g. minimum NDVI~0.2, thus sr_min~0.1). The global minimum can be determined from a histogram of the surface reflectance.

#slope_thrshld = maximum acceptable rate of change in surface reflectance for filtering bad values

#For the first entry in the dataset, set sr_filtered = sr and sr_filtered_prev = sr, also set doy_prev = doy

if sr < sr_min or abs(sr - sr_prev)/(doy - doy_prev) > slope_thrshld then sr_filtered = sr_filtered_prev else sr_filtered = sr

sr_slope = (sr_filtered - sr_filtered_prev)/(doy - doy_prev)

sr_intercept = sr - (sr_slope)*doy

sr_filtered_prev = sr_filtered

doy_prev = doy

# Gap Filling

#Create a continuous time-series of daily surface reflectances such that there is a value for each day of the year (366 days for leap years). In Excel I would create a list in a column with 365 (or 366) rows long then assign surface reflectances to each day using a lookup function. The assigned values would either be the good surface reflectances or if there were no good values, calculated ones based on the above slope and intercepts for that day. The calculated value will equal the good values so one can simply use the calculated value.

if modulus of year = 4 then doy_max = 366 else doy_max=365

for i = 1 to doy_max

     sr_filled = sr_slope*doy + sr_intercept

next i 

# Smoothing

#Once a continuous gap-filled time-series of surface reflectances is created, the values need to be smoothed to create gradual and physically realistic transitions between points in the time-series. This is done use an iterative 3-point moving average applied to the time-series 20 times. This involves calculating the mean surface reflectance  using the previous, current and next surface reflectances, then calculating the mean using the previous mean values for each iteration. This provides smoothing while preserving the dynamic features of the time-series and without introducing time-shifts or lags.

sr_prev = sr at previous timestep

sr_next = sr at next timestep

for the first timestep, sr_prev = sr

nsmooth = 20   #number of smoothing iterations to perform 

for i = 1 to n_smooth

     sr_mean = (sr_prev + sr + sr_next)/3

     sr = sr_mean

next i


# LE Equations

tc = 7   #timeconstant for running means (days)

dt = 1   #time interval between observations (days)

w = exp(-dt/tc)   #weighting coefficient for running mean (0<w<1)

tday = 0.25*tmin + 0.75*tmax   #average daytime temperature (deg C)

decl = 23.45*sin((2*pi*(284 + latitude)/365))   #solar declination (degrees) 

hourangle = acos(-tan(latitude)*tan(decl))   #hour angle (degrees) 

dl = 2*hourangle/15   #daylength (hours)

mm2wm2=((2500-2.361*tday)*100*100*100)/(1000*60*60*dl)   #multiplier to convert mm of water into average daytime W/m2