<jupyter_start><jupyter_text>DALEC-GRASS tutorial <NAME> | 02.06.2022 University of Edinburgh | Global Change Ecology Lab  This is a tutorial on using DALEC-Grass. The tutorial has 3 sections 1. DALEC-Grass code2. DALEC-Grass inputs3. Model data fusion with DALEC-Grass DALEC-GRASS DALEC-Grass is a grassland specific version of DALEC that was developed from DALEC_GIS_DFOL_FR.f90. It is written in fortran and its code is available on github (https://github.com/GCEL/DALEC-Grass). Since DALEC-Grass is work in progress the "latest version" is not on the github repo.<jupyter_code># import necessary packages
%matplotlib inline
from wand. image import Image as WImage
import salem
from convertbng.util import convert_bng, convert_lonlat
import geopandas as gpd 
import pandas as pd 
import numpy as np
import matplotlib.pyplot as plt 

# Plot a schematic of the DALEC-Grass model
img = WImage(filename='dalec_grass.pdf')
img

# Show some basic info on DALEC-Grass pools, fluxes, parameters and units
img = WImage(filename='dalec_pars.png')
img<jupyter_output><empty_output><jupyter_text>I am using DALEC-Grass in python. The fortran code will be integrated into the DALEC/CARDAMOM structure soon. There are python scipts packaged and available from github that allow the user to :1. source climate and optical/radar EO data form Alaska Satellite Facility, AWS and ECMWF (user accounts needed)2. process them into time-series for weather variables and LAI observations (storage/RAM demands)3. implement different probabilistic parameter optimisation algorithms (with/without paralellisation)The MDF_DALEC_GRASS package can be installed and DALEC_GRASS.f90 can be compilled by running these 2 lines :<jupyter_code>! python setup.py install # install package
! f2py -c DALEC_GRASS.f90 -m DALEC_GRASS # compile the DALEC_GRASS.f90 into a python object 
# see how the MDF_DALEC_GRASS folder looks
! ls -lh<jupyter_output><empty_output><jupyter_text>DALEC-Grass inputs The standard DALEC inputs needed with the GSI module : min/max T, srad, photoperiod, VPD, atm CO2 PLUS information on what is happening to the vegetation of a grassland field i.e. weekly reduction in LAI. Vegetation volume change inputs are a time-series of week-to-week reduction in LAI (in total m2.m-2 per week, can be zero or positive).The input_data_procuction module of the MDF_DALEC_GRASS package handles input data creation. Provided with a polygon (field limits) and a start/end date it will: 1. source raw met (ECMWF) and EO (ASF,AWS) data 2. process them into met and veg_reduction time-series; and into3. observational LAI time-series, based on fusing optical (sentinel-2) and radar (sentinel-1) dataSteps 1-3 take time (Step1: 20%, Step2: 30%, Step3: 50% of total duration).<jupyter_code># Example site location and limits
gdf = gpd.read_file('greatfield.geojson') # a geojson or shapefile with the location of the site
gdf = gdf.to_crs(epsg=4326)
gdf[['minx','miny','maxx','maxy']] = gdf.geometry.bounds

## Get the google map basemap
xdist = (gdf.minx.min() - gdf.maxx.max())*1.5
ydist = (gdf.miny.min() - gdf.maxy.max())*1.5
g = salem.GoogleVisibleMap(x=[gdf.minx.min()-xdist, gdf.maxx.max()+xdist],
                           y=[gdf.miny.min()-ydist, gdf.maxy.max()+ydist],
                           maptype='satellite', scale=1,
                           size_x=400, size_y=400)
ggl_img = g.get_vardata()

sm = salem.Map(g.grid, factor=1, countries=False)
sm.set_rgb(ggl_img,)  # add the background rgb image
cmap = plt.get_cmap('Blues')
for g, geo in enumerate(gdf.geometry) : sm.set_geometry(geo,alpha=0.75,facecolor=cmap(gdf.area.iloc[g]/(gdf.area.max()-gdf.area.min())))

# plot!
f, ax = plt.subplots(figsize=(10, 8))
ax.set_position([0.05, 0.06, 0.7, 0.9])
sm.visualize(ax=ax)
plt.show()
drivers = np.load('greatfield_M.npy') # load met/veg_red time-series
drivers.shape # see above for basic info on DALEC-Grass pools, fluxes, parameters and units 
years = drivers.shape[1]/52
years # how many years simulated 
# preview met/management time-series
inputs = pd.DataFrame(drivers)
inputs = inputs.T
inputs.index = pd.date_range('2017-01-01', periods = drivers.shape[1], freq='W') 
inputs = inputs.rename(columns={inputs.columns[1]: 'minT'}) 
inputs = inputs.rename(columns={inputs.columns[2]: 'maxT'}) 
inputs = inputs.rename(columns={inputs.columns[3]: 'srad'}) 
inputs = inputs.rename(columns={inputs.columns[7]: 'veg_red'}) 

plt.figure(figsize=(16, 4))
plt.subplot(1, 3, 1)
plt.title('temperature min/max')
plt.ylabel('$^{o}$C')
plt.plot(inputs.minT,label='min',color='blue')
plt.plot(inputs.maxT,label='max',color='red')
plt.xticks(rotation=90)
plt.legend()
plt.subplot(1, 3, 2)
plt.title('srad')
plt.ylabel('Mj m$^{2}$ d${-1}$')
plt.plot(inputs.srad)
plt.xticks(rotation=90)
plt.subplot(1, 3, 3)
plt.title('vegetation reduction')
plt.ylabel('LAI m$^{2}$ m$^{-2}$ w${-1}$')
plt.plot(inputs.veg_red)
plt.xticks(rotation=90)
plt.tight_layout()
plt.show()
# preview observational LAI time-series
laiobs = np.load('greatfield_O.npy')
laiobs = pd.DataFrame(laiobs)
laiobs.index = pd.date_range('2017-01-01', periods = len(laiobs), freq='W') 
laiobs = laiobs.rename(columns={inputs.columns[0]: 'lai'})
plt.plot(laiobs.lai)
plt.xticks(rotation=90)
plt.title('Observational LAI')
plt.ylabel('m$^{2}$ m$^{-2}$')
plt.show()<jupyter_output><empty_output><jupyter_text>Model data fusion with DALEC-GrassThe MDF algorithm is implemented by calling the run() function of the MDF_DALEC_GRASS package.The run() function can be parallelised to speed up the assimilation duration. You can uncomment the following block of code and run it. MDF.run() can take between 15-30 mins to complete.• Note that you need to provide the full path to the MDF_DALEC_GRASS folder (workingdir)<jupyter_code># import MDF
# MDF.run(workingdir="/full/path/to/MDF_DALEC_GRASS",sitename='greatfield')<jupyter_output><empty_output><jupyter_text>When the assimilation is completed we can draw samples from the posterior distribution and run DALEC-Grass in forward mode to obtain estimates of C pools, fluxes and balance. The following block of code completes the forward runs.• Note you need to edit the workingdir string (full path to the MDF_DALEC_GRASS folder).<jupyter_code># forward runs 
import numpy as np 
import pandas as pd 
import DALEC_GRASS
from sklearn.metrics import mean_squared_error
from math import sqrt
import datetime
import matplotlib.pyplot as plt 
from scipy import stats

workingdir = "/full/path/to/MDF_DALEC_GRASS"
sitename = 'greatfield'

# load posteriors
mcmc_outs = pd.read_csv('MDF_outs_%s.csv'%sitename)
mcmc_outs = mcmc_outs.sort_values(by='like1',ascending=False)
mcmc_outs = mcmc_outs[:100]
    
## Load drivers 
met       = np.array(np.load('%s/%s_M.npy' %(workingdir,sitename)),order="F")   
met       = np.append(met[:,:52],met,axis=1) ## add 1 spinup year 
met[0,:]  = np.arange(1,len(met[0,:])+1) ## re-create index 

## Fill DALEC-Grass input variables 
deltat   = np.zeros([(met.shape[1])]) + 7 # weekly runs 
nodays   = met.shape[1]
noyears  = int(nodays/float(52)) - 1 # weekly runs 
start    = 1     
finish   = int(nodays)
nopools  = 6     
nofluxes = 21 
nopars   = 34
nomet    = met.shape[0]
lat      = 50.77
VC       = 1

#### Load LAI observations 
obs_lai = np.load("%s/%s_O.npy" %(workingdir,sitename)) 

### Complete forward runs 
nosamples  = len(mcmc_outs)
lai_mult   = np.zeros([(nosamples),nodays]) * np.nan
gpp_mult   = np.zeros([(nosamples),nodays]) * np.nan 
npp_mult   = np.zeros([(nosamples),nodays]) * np.nan 
nee_mult   = np.zeros([(nosamples),nodays]) * np.nan 
nbe_mult   = np.zeros([(nosamples),nodays]) * np.nan 
graze_mult = np.zeros([(nosamples),nodays]) * np.nan 
soilC_mult = np.zeros([(nosamples),nodays+1]) * np.nan 
cut_mult   = np.zeros([(nosamples),nodays]) * np.nan  
abgC_mult  = np.zeros([(nosamples),nodays+1]) * np.nan 
rootC_mult = np.zeros([(nosamples),nodays+1]) * np.nan 
litterC_mult = np.zeros([(nosamples),nodays+1]) * np.nan 
animalCgas_mult = np.zeros([(nosamples),nodays]) * np.nan 
animalCmanure_mult = np.zeros([(nosamples),nodays]) * np.nan 
flux_autresp_mult = np.zeros([(nosamples),nodays]) * np.nan 
flux_hetresp_mult = np.zeros([(nosamples),nodays]) * np.nan 
flux_leaf_mult = np.zeros([(nosamples),nodays])  * np.nan
flux_stem_mult = np.zeros([(nosamples),nodays]) * np.nan 
flux_root_mult = np.zeros([(nosamples),nodays]) * np.nan 
flux_litter_mult = np.zeros([(nosamples),nodays]) * np.nan 
flux_som_mult = np.zeros([(nosamples),nodays]) * np.nan 

unc_DF = pd.DataFrame()

for i in range(0,(nosamples)):

	pars = np.array(mcmc_outs[mcmc_outs.columns[1:(nopars+1)]].iloc[i],order="F")
	lai,gpp,nee,pools,fluxes,rem = DALEC_GRASS.carbon_model_mod.carbon_model(start,finish,deltat,lat,met,pars,nopools,nofluxes,VC,nodays,nopars,nomet)
		
	lai_mult[i,:] = lai
	gpp_mult[i,:] = gpp
	npp_mult[i,:] = gpp - fluxes[:,2] 
	nee_mult[i,:] = nee
	graze_mult[i,:] = rem[0] 
	cut_mult[i,:] = rem[1] 
	soilC_mult[i,:] = pools[:,5]
	abgC_mult[i,:] =  pools[:,0] + pools[:,1]
	rootC_mult[i,:] = pools[:,2]
	litterC_mult[i,:] = pools[:,4]
	nbe_mult[i,:] = nee*7 + (rem[0] + rem[1]) - fluxes[:,18]
	animalCgas_mult[i,:] = fluxes[:,19] + fluxes[:,20]
	animalCmanure_mult[i,:] = fluxes[:,18]
	flux_autresp_mult[i,:] = fluxes[:,2] 
	flux_hetresp_mult[i,:] =  fluxes[:,12] + fluxes[:,13]
	flux_leaf_mult[i,:] = fluxes[:,3] 
	flux_stem_mult[i,:] =  fluxes[:,4] 
	flux_root_mult[i,:] = fluxes[:,5] 
	flux_litter_mult[i,:] =  fluxes[:,9] + fluxes[:,10] + fluxes[:,11]
	flux_som_mult[i,:] =  fluxes[:,14] 

cut_mult[cut_mult==0] = np.nan

for i in range(nodays): 
	unc_DF = unc_DF.append({
			  'lai': stats.bayes_mvs(lai_mult[:,i][~np.isnan(lai_mult[:,i])],alpha=0.99)[0][0] , 
			  'lai_std': stats.bayes_mvs(lai_mult[:,i][~np.isnan(lai_mult[:,i])],alpha=0.99)[2][0],
			  'gpp': stats.bayes_mvs(gpp_mult[:,i][~np.isnan(gpp_mult[:,i])],alpha=0.99)[0][0] , 
			  'gpp_std': stats.bayes_mvs(gpp_mult[:,i][~np.isnan(gpp_mult[:,i])],alpha=0.99)[2][0],
			  'nee': stats.bayes_mvs(nee_mult[:,i][~np.isnan(nee_mult[:,i])],alpha=0.99)[0][0] , 
			  'nee_std': stats.bayes_mvs(nee_mult[:,i][~np.isnan(nee_mult[:,i])],alpha=0.99)[2][0],
			  'npp': stats.bayes_mvs(npp_mult[:,i][~np.isnan(npp_mult[:,i])],alpha=0.99)[0][0] , 
			  'npp_std': stats.bayes_mvs(npp_mult[:,i][~np.isnan(npp_mult[:,i])],alpha=0.99)[2][0],
			  'nbe': stats.bayes_mvs(nbe_mult[:,i][~np.isnan(nbe_mult[:,i])],alpha=0.99)[0][0] , 
			  'nbe_std': stats.bayes_mvs(nbe_mult[:,i][~np.isnan(nbe_mult[:,i])],alpha=0.99)[2][0],
			  'graze_prob': len(graze_mult[:,i][graze_mult[:,i]>0])/float(len(graze_mult[:,i])),
			  'graze': stats.bayes_mvs(graze_mult[:,i][~np.isnan(graze_mult[:,i])],alpha=0.99)[0][0] , 
			  'graze_std': stats.bayes_mvs(graze_mult[:,i][~np.isnan(graze_mult[:,i])],alpha=0.99)[2][0],
			  'cut_prob': len(cut_mult[:,i][cut_mult[:,i]>0])/float(len(cut_mult[:,i])),
			  'cut_mean' : np.nanmean(cut_mult[:,i]),
			  'cut_std' : np.nanstd(cut_mult[:,i]),	  
			  'soilC': stats.bayes_mvs(soilC_mult[:,i][~np.isnan(soilC_mult[:,i])],alpha=0.99)[0][0] , 
			  'soilC_std': stats.bayes_mvs(soilC_mult[:,i][~np.isnan(soilC_mult[:,i])],alpha=0.99)[2][0],
			  'abgC': stats.bayes_mvs(abgC_mult[:,i][~np.isnan(abgC_mult[:,i])],alpha=0.99)[0][0] , 
			  'abgC_std': stats.bayes_mvs(abgC_mult[:,i][~np.isnan(abgC_mult[:,i])],alpha=0.99)[2][0],
			  'rootC': stats.bayes_mvs(rootC_mult[:,i][~np.isnan(rootC_mult[:,i])],alpha=0.99)[0][0] , 
			  'rootC_std': stats.bayes_mvs(rootC_mult[:,i][~np.isnan(rootC_mult[:,i])],alpha=0.99)[2][0],
			  'litterC': stats.bayes_mvs(litterC_mult[:,i][~np.isnan(litterC_mult[:,i])],alpha=0.99)[0][0] , 
			  'litterC_std': stats.bayes_mvs(litterC_mult[:,i][~np.isnan(litterC_mult[:,i])],alpha=0.99)[2][0],
			  'animalCgas': stats.bayes_mvs(animalCgas_mult[:,i][~np.isnan(animalCgas_mult[:,i])],alpha=0.99)[0][0] , 
			  'animalCgas_std': stats.bayes_mvs(animalCgas_mult[:,i][~np.isnan(animalCgas_mult[:,i])],alpha=0.99)[2][0],
			  'animalCmanure': stats.bayes_mvs(animalCmanure_mult[:,i][~np.isnan(animalCmanure_mult[:,i])],alpha=0.99)[0][0] , 
			  'animalCmanure_std': stats.bayes_mvs(animalCmanure_mult[:,i][~np.isnan(animalCmanure_mult[:,i])],alpha=0.99)[2][0],
			  'autresp': stats.bayes_mvs(flux_autresp_mult[:,i][~np.isnan(flux_autresp_mult[:,i])],alpha=0.99)[0][0], 
			  'autresp_std': stats.bayes_mvs(flux_autresp_mult[:,i][~np.isnan(flux_autresp_mult[:,i])],alpha=0.99)[2][0],
			  'hetresp': stats.bayes_mvs(flux_hetresp_mult[:,i][~np.isnan(flux_hetresp_mult[:,i])],alpha=0.99)[0][0], 
			  'hetresp_std': stats.bayes_mvs(flux_hetresp_mult[:,i][~np.isnan(flux_hetresp_mult[:,i])],alpha=0.99)[2][0],			  
			  'F_leaf': stats.bayes_mvs(flux_leaf_mult[:,i][~np.isnan(flux_leaf_mult[:,i])],alpha=0.99)[0][0], 
			  'F_leaf_std': stats.bayes_mvs(flux_leaf_mult[:,i][~np.isnan(flux_leaf_mult[:,i])],alpha=0.99)[2][0],
			  'F_stem': stats.bayes_mvs(flux_stem_mult[:,i][~np.isnan(flux_stem_mult[:,i])],alpha=0.99)[0][0], 
			  'F_stem_std': stats.bayes_mvs(flux_stem_mult[:,i][~np.isnan(flux_stem_mult[:,i])],alpha=0.99)[2][0],
			  'F_root': stats.bayes_mvs(flux_root_mult[:,i][~np.isnan(flux_root_mult[:,i])],alpha=0.99)[0][0], 
			  'F_root_std': stats.bayes_mvs(flux_root_mult[:,i][~np.isnan(flux_root_mult[:,i])],alpha=0.99)[2][0],			  			  
			  'F_litter': stats.bayes_mvs(flux_litter_mult[:,i][~np.isnan(flux_litter_mult[:,i])],alpha=0.99)[0][0], 
			  'F_litter_std': stats.bayes_mvs(flux_litter_mult[:,i][~np.isnan(flux_litter_mult[:,i])],alpha=0.99)[2][0],
			  'F_som': stats.bayes_mvs(flux_som_mult[:,i][~np.isnan(flux_som_mult[:,i])],alpha=0.99)[0][0], 
			  'F_som_std': stats.bayes_mvs(flux_som_mult[:,i][~np.isnan(flux_som_mult[:,i])],alpha=0.99)[2][0]			  
			  },ignore_index=True)

unc_DF.index = pd.date_range('2016-01-08', periods = nodays, freq='7D')
unc_DF = unc_DF['2017':]

# what did the forward runs produce 
unc_DF # the dataframe holding DALEC-Grass predictions (mean and std) for C pools, fluxes and removals<jupyter_output><empty_output><jupyter_text>Plotting the results of forward runs<jupyter_code>## Leaf Area Index 
plt.figure(figsize=(8, 4))
plt.title('Leaf Area Index')
plt.ylabel('m$^{2}$ m$^{-2}$')
plt.plot(unc_DF.lai,color='green',label='sim')
plt.fill_between(unc_DF.index, unc_DF.lai + 2*unc_DF['lai_std'], unc_DF.lai - 2*unc_DF['lai_std'], color='green', alpha=0.33)
plt.scatter(unc_DF.index,obs_lai,marker='+',color='red',label='obs')
plt.xticks(rotation=90)
plt.legend()
plt.show()
## NEE 
plt.figure(figsize=(12, 4))
plt.subplot(1,2,1)
plt.title('NEE')
plt.ylabel('g C m$^{-2}$ d$^{-1}$')
plt.plot(unc_DF.nee,color='black')
plt.fill_between(unc_DF.index, unc_DF.nee + 2*unc_DF['nee_std'], unc_DF.nee - 2*unc_DF['nee_std'], color='black', alpha=0.33)
plt.xticks(rotation=90)
plt.subplot(1,2,2)
plt.plot(unc_DF.nee.cumsum(),color='black')
plt.ylabel('g C m$^{-2}$')
plt.xticks(rotation=90)
plt.show()
## AGB 
plt.figure(figsize=(8, 4))
plt.title('above/below - ground Biomass')
plt.fill_between(unc_DF.index, unc_DF.abgC, unc_DF.abgC*0 , color='black', alpha=0.33,label='AGB')
plt.fill_between(unc_DF.index, -unc_DF.rootC, unc_DF.rootC*0 , color='brown', alpha=0.33,label='roots')
plt.ylabel('g C m$^{-2}$')
plt.xticks(rotation=90)
plt.legend()
plt.show()
## Grazing 
plt.figure(figsize=(8, 4))
plt.title('Predicted livestock density')
plt.xlabel('month')
plt.ylabel('LSU per ha')
plt.bar(list((unc_DF.graze.resample('M').sum()).index),(unc_DF.graze.resample('M').sum()*21/float(600*pars[30]))/float(30),width=20,color='black')
plt.xticks(rotation=90)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.186846326596086, 'fraction_numerical': 0.030093038177735003, 'mean_word_length': 5.854001759014952, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 28, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '197834', 'license': 'CC0-1.0', 'n_tokens_mistral': 7002, 'n_tokens_neox': 6483, 'n_words': 1083, 'path': 'dalec_grass_tutorial.ipynb', 'repo_name': 'vmyrgiotis/MDF_DALEC_Grass', 'size': '15596'}