This page was generated from source/Parameters/tutorials/param1-LAI.ipynb. Interactive online version: Slideshow:

10.3.1. Leaf Area Index related parameters

Necessary packages for analysis:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime
from pathlib import Path
import supy as sp
import os
from shutil import copyfile
import geopandas as gpd
import pickle
pd.options.mode.chained_assignment = None
from IPython.display import clear_output
import warnings
warnings.filterwarnings('ignore')

This is a custom package specifically designed for this analysis. It contains various functions for reading and computing and plotting.

from alb_LAI_util import generate_array_dif, generate_array_same, modify_attr,func_parse_date

10.3.1.1. Parameters needed to calcuate LAI for non-crop vegetation (see Background https://umep-workshop.readthedocs.io/en/latest/Parameters/CalcBG.html)

• \(BaseT\) -

• \(BaseTe\)

• \(LAI_{MAX}\)

• \(LAI_{Min}\)

10.3.1.2. Necessary functions for analysis

10.3.1.2.1. Read in the data (including LAI data)

def read_data(year, name):

    df = pd.read_csv('data/MODIS_LAI_AmeriFlux/statistics_Lai_500m-' + name +
                     '.csv')
    if name!='crop':
        df.columns = ['product'
                      ] + [i.split(' ')[1] for i in df.columns if i != 'product']

    df = df.filter(['modis_date', 'value_mean'])

    df_period = df[[i.startswith('A' + str(year)) for i in df.modis_date]]
    df_period.loc[:, 'DOY'] = [
        int(i.split('A' + str(year))[1]) for i in df_period.modis_date
    ]
    df_period = df_period.set_index('DOY')

    copyfile("./runs/data/" + name + "_" + str(year) + "_data_60.txt",
             "runs/run/Input/Kc_2012_data_60.txt")
    df_forcing = pd.read_csv('runs/run' + '/Input/' + 'Kc' + '_' + '2012' +
                             '_data_60.txt',
                             sep=' ',
                             parse_dates={'datetime': [0, 1, 2, 3]},
                             keep_date_col=True,
                             date_parser=func_parse_date)

    all_sites_info = pd.read_csv('site_info.csv')
    site_info = all_sites_info[all_sites_info['Site Id'] == name]
    df = pd.DataFrame({
        'Site': [name],
        'Latitude': [site_info['Latitude (degrees)']],
        'Longitude': [site_info['Longitude (degrees)']]
    })

    path_runcontrol = Path('runs/run' + '/') / 'RunControl.nml'
    df_state_init = sp.init_supy(path_runcontrol)

    df_state_init ,level= modify_attr(df_state_init, df, name)

    grid = df_state_init.index[0]
    df_forcing_run = sp.load_forcing_grid(path_runcontrol, grid)

    df_output, df_state_final = sp.run_supy(df_forcing_run,
                                            df_state_init,
                                            save_state=False)

    return df_output, df_state_final, df_state_init, df_period, grid, df_forcing_run,level

10.3.1.2.2. Calculating some of the variables of the models (LAI,GDD,SDD,Tair)

def calc_vars(df_output, grid, df_forcing_run,year,level):

    df_output_2 = df_output.loc[grid, :]
    df_output_2 = df_output_2[df_output_2.index.year >= year]

    lai_model = pd.DataFrame(df_output_2.SUEWS.LAI)
    a = lai_model.index.strftime('%j')
    lai_model['DOY'] = [int(b) for b in a]

    if(level==1):
        nameGDD='GDD_DecTr'
        nameSDD='SDD_DecTr'
    elif(level==0):
        nameGDD='GDD_EveTr'
        nameSDD='SDD_EveTr'
    if(level==2):
        nameGDD='GDD_Grass'
        nameSDD='SDD_Grass'

    GDD_model = pd.DataFrame(df_output_2.DailyState[nameGDD])
    a = GDD_model.index.strftime('%j')
    GDD_model['DOY'] = [int(b) for b in a]
    GDD_model = GDD_model.dropna()
    GDD_model = GDD_model.iloc[1:]

    SDD_model = pd.DataFrame(df_output_2.DailyState[nameSDD])
    print(SDD_model)
    a = SDD_model.index.strftime('%j')
    SDD_model['DOY'] = [int(b) for b in a]
    SDD_model = SDD_model.dropna()
    SDD_model = SDD_model.iloc[1:]

    Tair = df_forcing_run.Tair.resample('1d', closed='left',
                                        label='right').mean()
    Tair = pd.DataFrame(Tair)
    a = Tair.index.strftime('%j')
    Tair['DOY'] = [int(b) for b in a]

    return lai_model, GDD_model, SDD_model, Tair,nameGDD,nameSDD

10.3.1.2.3. Calculate and plot the LAI for the site being analysed. The data are saved to a pickle file.

def LAI_tune(year,name):
    df_output,df_state_final,df_state_init,df_period,grid,df_forcing_run,level = read_data(year,name)
    lai_model,GDD_model,SDD_model,Tair,nameGDD,nameSDD=calc_vars(df_output,grid,df_forcing_run,year,level)
    clear_output()

    with open('outputs/LAI/'+name+'-'+str(year)+'-MODIS.pkl','wb') as f:
        pickle.dump(df_period, f)
    with open('outputs/LAI/'+name+'-'+str(year)+'-Model.pkl','wb') as f:
        pickle.dump(lai_model, f)

    attrs=[
    df_state_init.loc[:,'laimin'].loc[grid][level],
    df_state_init.loc[:,'laimax'].loc[grid][level],
    df_state_init.loc[:,'basete'].loc[grid][level],
    df_state_init.loc[:,'baset'].loc[grid][level]
    ]

    with open('outputs/LAI/'+name+'-attrs.pkl','wb') as f:
            pickle.dump(attrs, f)

    plt.rcParams.update({'font.size': 15})
    fig,axs=plt.subplots(2,1,figsize=(20,10))
    ax=axs[0]
    ax.plot(lai_model.DOY,lai_model.LAI,color='b',label='SUEWS-LAI')
    df_period.value_mean.plot(color='r',label='MODIS-LAI',ax=ax)
    ax.set_ylabel('LAI')
    ax.legend()
    ax.set_title(name+'-'+str(year))
    ax.set_xlabel('')

    ax=axs[1]
    plt.scatter(Tair.DOY,Tair.Tair,color='k')



    try:
        max_y=30
        for gdd_day in [90,100,110,120,130,140,150,170,180,210,260]:
            a=GDD_model[GDD_model.DOY==gdd_day][nameGDD]
            plt.plot([gdd_day,gdd_day],[-15,max_y],color='r')
            plt.annotate(str(np.round(a.values[0],0)),(gdd_day-5,-14),color='r',rotation=90)


        for sdd_day in [150,170,180,200,250,270,300,310,320,330,340,350,360]:
            a=SDD_model[SDD_model.DOY==sdd_day][nameSDD]
            plt.plot([sdd_day,sdd_day],[-15,max_y],color='b')
            plt.annotate(str(np.round(a.values[0],0)),(sdd_day-5,-14),color='b',rotation=90)
    except:
        pass

    ax.plot([0,365],[df_state_init.basete.iloc[0][0],df_state_init.basete.iloc[0][0]],label='BaseTe',color='k')
    ax.plot([0,365],[df_state_init.baset.iloc[0][0],df_state_init.baset.iloc[0][0]],'--k',label='BaseT')
    ax.legend()
    plt.ylabel('Tair (C)')
    plt.xlabel('DOY')
    ax.set_xlim(left=0,right=365)

10.3.1.2.4. Now evaluate the LAI parameters at other ‘Test’ sites. Calculate and plot LAI (and save data to pickle files)

def LAI_test(years,name):

    plt.rcParams.update({'font.size': 12})
    fig,axs=plt.subplots(len(years),2,figsize=(20,13))

    counter=-1
    for year in years:
        print(name+'-'+str(year))
        counter=counter+1
        df_output,df_state_final,df_state_init,df_period,grid,df_forcing_run,level = read_data(year,name)
        lai_model,GDD_model,SDD_model,Tair,nameGDD,nameSDD=calc_vars(df_output,grid,df_forcing_run,year,level)
        clear_output()

        with open('outputs/LAI/'+name+'-'+str(year)+'-MODIS.pkl','wb') as f:
            pickle.dump(df_period, f)
        with open('outputs/LAI/'+name+'-'+str(year)+'-Model.pkl','wb') as f:
            pickle.dump(lai_model, f)


        ax=axs[counter][0]
        ax.plot(lai_model.DOY,lai_model.LAI,color='b',label='SUEWS-LAI')
        df_period.value_mean.plot(color='r',label='MODIS-LAI',ax=ax)
        ax.set_ylabel('LAI')
        if counter==0:
            ax.legend()
        if counter!=len(years)-1:
            ax.set_xlabel('')
        ax.set_title(name+'-'+str(year))




        ax=axs[counter][1]
        ax.scatter(Tair.DOY,Tair.Tair,color='k')
        try:
            max_y=30
            for gdd_day in [90,110,130,150]:
                a=GDD_model[GDD_model.DOY==gdd_day][nameGDD]
                ax.plot([gdd_day,gdd_day],[-15,max_y],color='r')
                ax.annotate(str(np.round(a.values[0],0)),(gdd_day-10,-14),color='r',rotation=90)


            for sdd_day in [300,320,340,360]:
                a=SDD_model[SDD_model.DOY==sdd_day][nameSDD]
                ax.plot([sdd_day,sdd_day],[-15,max_y],color='b')
                ax.annotate(str(np.round(a.values[0],0)),(sdd_day-10,-14),color='b',rotation=90)
        except:
            pass
        ax.plot([0,365],[df_state_init.basete.iloc[0][0],df_state_init.basete.iloc[0][0]],label='BaseTe',color='k')
        ax.plot([0,365],[df_state_init.baset.iloc[0][0],df_state_init.baset.iloc[0][0]],'--k',label='BaseT')
        if counter==0:
            ax.legend()
        ax.set_ylabel('Tair (C)')
        if counter==len(years)-1:
            ax.set_xlabel('DOY')
        ax.set_xlim(left=0,right=365)
        ax.set_title(name+'-'+str(year))

    plt.savefig('figs/'+name+'-LAI.png',dpi=300,bbox_inches = 'tight',pad_inches = 0.01)

10.3.1.3. Plot for all sites

10.3.1.3.1. DecTr

10.3.1.3.1.1. US-MMS

name='US-MMS'
year=2017
LAI_tune(year,name)

name='US-MMS'
years=[2010,2012,2016]
LAI_test(years,name)

10.3.1.3.2. EveTr

10.3.1.3.2.1. CA-Obs

name = 'CA-Obs'
year = 2008
LAI_tune(year, name)

name = 'CA-Obs'
years = [2003, 2005, 2006]
LAI_test(years, name)

10.3.1.3.3. Grass

10.3.1.3.3.1. US-AR1

name = 'US-AR1'
year = 2012
LAI_tune(year, name)

name = 'US-AR1'
years = [2010, 2011]
LAI_test(years, name)