Interannual variability

[dlebauer]

from Tao:

I think it would be very interesting to quantify the inter-annual climate changes on perennial crop production and bioenergy supply chain development.

The key issues we need to work on further would be:

1. How to select random climate patterns based on historical climate data sets and how to assign probability for each climate scenario?
2. Develop probability distribution functions of accumulated perennial corp yield (y axis: probability, x axis: stepwise accumulated yield ( the step size depends on the variations of yield results)
3. Select representative yield samples (or mean yield in each yield range) with its probability density for optimization model inputs

Please let me know if you have any comments and any help on the model runs and data analysis. I would like to focus on this study for the following month.

---

references:

Wolfram Schlenker| 2006 Inter-annual Weather Variation and Crop Yields S. Asseng et al. 2013 Uncertainty in simulating wheat yields under climate change

[dlebauer]

There are two very different analyses in the papers you sent : the effects of climate change and the effects of interannual variability.

The scope of the Agmip Climate change scenarios (Asseng et al, predicting to 2100) is beyond the scope of this paper. I do have the met data, we just need an author

The analyses in the paper by Schlenker is basically a regression of yield on climate variables (and location). We can similarly calculate the relationship between yield and climate (e.g. the slope of temperature and precipitation responses).

e.g.

lm(yield ~ temperature + precipitation + random(site|year))

I think this will be more meaningful that running the model under 5 global climate models x 4 CO2 emissions scenarios, especially because the global climate models used to generate these analyses are explicitly intended to predict mean fields and decadal to century scale trends rather than inter-annual variability.

1. How to select random climate patterns based on historical climate data sets and how to assign probability for each climate scenario?

I'll take a closer look, but to a first appoximation,

• p(any year) = 1 / (total # years)
• p(any set of 15 years) ~ 1 / (choose 15 from 100 years)

To select the years, I will put the following at line 23 of the run.biocro function:

if(met.uncertainty = TRUE){
   met <- load.cfmet(met.nc, lat = lat, lon = lon, start.date = "1979-01-01", 
        end.date = "2010-12-31")
  year.sample <- met[,list(year = sample(unique(year), size = 15))] 
  met <- met[year %in% year.sample$year]
}

Then to loop over this, say 500 times, add the following to regionalbiocro.Rscript (and perhaps rename to regional_metuncertainty_biocro.Rscript

starting at line 42

for(point in 1:nrow(latlon)){

add

for (i in 1:500){
  for(point in 1:nrow(latlon)){
    set.seed(i)
    ....
    ## change def. of point.outdir
    point.outdir <- file.path(outdir, i, paste0(lat, 'x', lon)

1. Develop probability distribution functions of accumulated perennial corp yield (y axis: probability, x axis: stepwise accumulated yield ( the step size depends on the variations of yield results)

## compute cumulative probability as fn of yield
ecdf(yield)
## find a distribution that fits simulated yield
lapply(c('gamma', 'lnorm', 'weibull', 'norm') function(x) fitdist(yield, x)

1. Select representative yield samples (or mean yield in each yield range) with its probability density for optimization model inputs

quantile(yield, c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9)
## probability density: ## or best fit distribution / parameters from #2
dnorm(quantile(yield, c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9), mu_hat, sigma_hat)

[dlebauer]

Please let me know if you have any comments and any help on the model runs and data analysis. I would like to focus on this study for the following month.

Would be great if you or Hao wanted to do the above proposed computations / changes.

[vnbisn]

Stochastic weather generators can be used to generate artificial weather. This is one such paper: http://onlinelibrary.wiley.com/doi/10.1029/2006WR005364/abstract

[haohuuiuc]

David, can you instruct us on how we can load the met data in R? I saw it in run.biocro function that you filtered met data based on sampling from different years.

[dlebauer]

Hao - how familiar are you with R?

Start with:

debugonce(run.biocro)

Then run the code in regional_pecan_workflow.Rmd

It will go through this:

met.nc <- nc_open(settings$run$inputs$met$path)

## equals 
# met.nc <- nc_open("/home/a-m/dlebauer/ebimodeling/met/narr/threehourly/illinois.nc")
lat <- 40
lon <- -88
start.date <- '1979-01-01'
end.date <- '1990-01-01'

Then Rstudio will launch you into a debugger environment to walk through the function. At this point I find it more useful to begin editing the source code file in the pecan/models/biocro/R/run.biocro.R directory.

In either case, you will end up inside run.biocro looking at code like this:

met <- load.cfmet(met.nc, lat = lat, lon = lon, start.date = start.date, end.date = end.date)
met.hr <- cfmet.downscale.time(cfmet = met, output.dt = 1)
biocro.met <- cf2biocro(met.hr)

Which is where the code to randomly sample years above would be added. If the illinois.nc file doesn't start at 1900, I can check my other code.

In this case, I don't think it is necessary to 'prepare' new met files, just to enable BioCro to be run on a random sample of years. By adding a few lines to the run.biocro function. (note the weather should be sub-sampled before downscaling). I think that 500 15 year runs would be sufficient.

---

@vnbisn the code you sent is indeed an exceptional weather generator. Just what we need - it generates hourly estimates for many variables including humidity, direct and diffuse solar radiation, PAR. However ... a bit more than I can take on at this time unless there is a logical flaw in sampling from the historical record to estimate climate uncertainty. However, it looks like a great method for handling the climate model forecasts that we have.

http://www-personal.umich.edu/~ivanov/HYDROWIT/Models.html

source: http://www-personal.umich.edu/~ivanov/HYDROWIT/Models_files/AWE_GEN_COMPLETE.zip

[haohuuiuc]

Thanks David! When I started with debugonce(run.biocro) and then ran the code in regional_pecan_workflow.Rmd, it did not bring me to a debugging environment, but directly submit jobs to ROGER. I figured out that the code reads the met data on ROGER because when I print out settings$run$inputs$met$path in R, it returns me /gpfs_scratch/haohu3//illinois.nc, it that part correct?

[dlebauer]

yep. then you can start from

met.nc <- nc_open("/gpfs_scratch/haohu3//illinois.nc")

On Thu, Sep 3, 2015 at 11:50 PM, Hao Hu notifications@github.com wrote:

Thanks David! When I started with debugonce(run.biocro) and then ran the code in regional_pecan_workflow.Rmd, it did not bring me to a debugging environment, but directly submit jobs to ROGER. I figured out that the code reads the met data on ROGER because when I print out settings$run$inputs$met$path in R, it returns me /gpfs_scratch/haohu3// illinois.nc, it that part correct?

— Reply to this email directly or view it on GitHub https://github.com/ebimodeling/biocro_regional/issues/15#issuecomment-137648823 .

[dlebauer]

I apologize. I confused the workflow script with the script that reads the configuration file, opens the weather + soil, and then runs BioCro

On Roger, you can start in this file:

> settings$model$binary
[1] "/home/a-m/dlebauer/dev/pecan/models/biocro/inst//regionalbiocro.Rscript"

and start stepping through the steps. the outdir and rundir are found in any of the run/<runid>/job.sh scripts

For example, see

/home/dlebauer/pecan_remote/mxg_calibration/Orr_Research/run/17203/job.sh

You can just run that script from anywhere. But you can also look in and see how it is calling biocro:

This line has three arguments: the r script, the input directory, and the output directory

/home/dlebauer/dev/pecan/models/biocro/inst//pointbiocro.Rscript /home/dlebauer/
pecan_remote/mxg_calibration/Orr_Research/run/17203 /home/dlebauer/pecan_remote/
mxg_calibration/Orr_Research/out/17203

So ... you need to open /home/dlebauer/dev/pecan/models/biocro/inst//pointbiocro.Rscript, set

rundir  <- "/home/dlebauer/
pecan_remote/mxg_calibration/Orr_Research/run/17203"

outdir  <- "/home/dlebauer/
pecan_remote/mxg_calibration/Orr_Research/out/17203"

then

debugonce(run.biocro)