Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
0 parents
commit f68ca27
Showing
73 changed files
with
9,087 additions
and
0 deletions.
There are no files selected for viewing
Binary file not shown.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,11 @@ | ||
Package: CropWatR | ||
Type: Package | ||
Title: Crop-water balances according to algorithms described in the FAO Irrigation and drainage paper 56 | ||
Version: 1.0 | ||
Date: 2015-09-06 | ||
Author: Jacob Teter | ||
Requires: rgdal, sp, raster | ||
Suggests: rgeos, rasterVis, ggplot2 | ||
Maintainer: Jacob Teter <jacob.teter@gmail.com> | ||
Description: Functions to enable calculation of FAO Penman-Monteith reference evapotranspiration, crop- and soil-water balances for agricultural crops. GIS-enabled implementation (using R's basic spatial packages, including rgeos, rgdal sp, raster) at a daily time-step, including calibration methods for irrigation scheduling and variable crop planting and harvesting dates. Crop-water models according to the FAO 56 document 'Crop evapotranspiration - Guidelines for computing crop water requirements' | ||
License: GPL-2 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1 @@ | ||
exportPattern("^[[:alpha:]]+") |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,51 @@ | ||
Calc.Basal.Crop.Coeff <- | ||
function(Croplayer){ | ||
load(paste0(Intermediates, paste('Growing.Season', Croplayer, 'U2.final_', 'Rdata', sep = '.'))); U2 <- Growing.Season; rm(Growing.Season) | ||
load(paste0(Intermediates, paste('Growing.Season', Croplayer, 'MNRH_', 'Rdata', sep = '.'))); RHmin <- Growing.Season; rm(Growing.Season) | ||
load(paste0(Intermediates,paste('CropsList', Croplayer, 'Rdata', sep = '.'))) | ||
# check x, y coordinates listing of the two files: # all.equal(lapply(U2, function(x) dim(x)), lapply(RHmin, function(x) dim(x))) | ||
all.equal(lapply(U2, function(x) x[,c('x','y')]), lapply(RHmin, function(x) x[,c('x','y')])) | ||
all.equal(names(U2), names(Daily.Crops.list)) | ||
|
||
GS.U2 <- lapply(U2, function(x) x[,(grep('layer', names(x)))]) | ||
GS.RHmin <- lapply(RHmin, function(x) x[,(grep('layer', names(x)))]) | ||
|
||
sapply(GS.U2, dim) | ||
sapply(GS.RHmin, dim) | ||
all.equal(sapply(GS.U2, dim), sapply(GS.RHmin, dim)) | ||
|
||
### Subset only mid and late growing season | ||
Wind_term <- lapply(GS.U2, function(x) 0.04*(x[]-2)) | ||
RH_term <- lapply(GS.RHmin, function(x) 0.004*(x[]-45)) | ||
all.equal(lapply(Wind_term, function(x) c(x$x, x$y)), lapply(RH_term, function(x) c(x$x, x$y))) | ||
|
||
Max.season.heights <- lapply(Daily.Crops.list, function(x) tapply(x$day_height, x$season.ID, max)) | ||
Season.IDs <- lapply(Daily.Crops.list, function(x) x$season.ID); Plant_heights <- Season.IDs | ||
for (i in 1:length(Max.season.heights)){ | ||
Plant_heights[[i]] <- Max.season.heights[[i]][match(Season.IDs[[i]], names(Max.season.heights[[i]]))] | ||
} | ||
Kcb <- lapply(Daily.Crops.list, function(x) x$daily_Kcb) # dim(Kcb[[1]]); length(Kcb[[1]]) | ||
print('done plant heights') | ||
height_term <- lapply(Plant_heights, function(x) (x[]/3)^0.3) | ||
Kcb.corrected <- Wind_term | ||
summary(Kcb.corrected[[1]]) | ||
for (i in 1:length(Wind_term)){ | ||
for (j in 1:length(height_term[[i]])){ | ||
Kcb.corrected[[i]][,j] <- Kcb[[i]][j] + (Wind_term[[i]][,j] - RH_term[[i]][,j]) * height_term[[i]][j] | ||
} | ||
} | ||
print('done correction term') | ||
Mid.Late.Season.cuts <- lapply(Daily.Crops.list, function(x) which(x$season.ID == 3 | x$season.ID == 4)) | ||
|
||
Before <- Kcb.corrected | ||
for (i in 1:length(Kcb.corrected)){ # lots of looping, so it's slow. It takes: ~30 seconds | ||
for (j in 1:length(Kcb.corrected[[i]])){ | ||
Kcb.corrected[[i]][,j][which(GS.RHmin[[i]][,j] > 20 & GS.RHmin[[i]][,j] < 80)] <- Kcb[[i]][j] | ||
Kcb.corrected[[i]][,j][which(GS.U2[[i]][,j] > 1 & GS.U2[[i]][,j] < 6)] <- Kcb[[i]][j] | ||
} | ||
} | ||
all.equal(Before, Kcb.corrected) # not equal, but minor differences, as is should be. | ||
all.equal(lapply(Before, function(x) dim(x)), lapply(Kcb.corrected, function(x) dim(x))) | ||
all.equal(names(Before), names(Kcb.corrected)) | ||
save(Kcb.corrected, file = paste0(Intermediates, paste('Kcb.corrected', Croplayer, 'Rdata', sep = '.'))) | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,19 @@ | ||
Calc.Fc.Few <- | ||
function(Croplayer){ | ||
load(paste0(Intermediates, paste('KcMax', Croplayer, 'Rdata', sep = '.'))) | ||
GR.KcMax <- lapply(KcMax, function(x) x[,(grep('layer', names(x)))]) # clip off the coordinates for analysis | ||
load(paste0(Intermediates, paste('CropsList', Croplayer, 'Rdata', sep = '.'))) | ||
# str(Daily.Crops.list[[1]]) | ||
DayHeight <- lapply(Daily.Crops.list, function(x) x$day_height) | ||
load(paste0(Intermediates, paste('Kcb.corrected', Croplayer, 'Rdata', sep = '.'))) # Kcb.corrected | ||
|
||
KcMin <- lapply(Kcb.corrected, function(x) c(rep(min(x)-0.01, times = length(x)))); Fc <- GR.KcMax | ||
for (i in 1:length(Fc)){ | ||
for (j in 1:length(DayHeight[[i]])){ | ||
Fc[[i]][,j] <- ((Kcb.corrected[[i]][,j] - KcMin[[i]][j])/(GR.KcMax[[i]][,j] - KcMin[[i]][j]))^(1+0.5*DayHeight[[i]][j]) | ||
} | ||
} | ||
Few <- Fc | ||
Few <- lapply(Fc, function(x) 1-x[]) # This suffices for RAINFED crops | ||
save(Few, file = paste0(Intermediates, paste('Few', Croplayer, 'Rdata', sep = '.'))) | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,112 @@ | ||
Calculate.ETo <- | ||
function(Elevation, MaxTemperature, MinTemperature, MeanTemperature, Precipitation, VP, MaxRH, MinRH, Wind, SolarRad, Filename){ | ||
# Inputs are file names of the input raster files (daily weather files, elevation, etc.) | ||
|
||
### I. Load CONSTANTS | ||
LHV <- 2.45 # latent heat of vaporization [MJ/kg] | ||
Cp <- 1.013*10^-3 # specific heat at constant pressure | ||
E <- 0.622 # P from above (1 millibar = 100 pascals = 0.1 kilopascals) | ||
Alpha <- 0.23 # albedo or canopy reflection coefficient, which is 0.23 for the hypothetical grass reference crop [dimensionless], | ||
G <-0 #### 10. Soil heat flux (G) # G [MJ/m^2*day], "may be ignored for daily timesteps". | ||
Gsc <- 0.0820 # Solar constant (MJ/m^2*min) | ||
As0 <- 0.75 # as+bs fraction of extraterrestrial radiation reaching the earth on clear days (n = N). | ||
Bs0 <- 2*10^-5 | ||
Sigma <- 4.903*10^-9 # Stefan-Boltzmann constant [ 4.903 10^-9 MJ K^-4 m^-2 day^-1], | ||
|
||
### II. Load & Mask INPUT files: | ||
|
||
# In the Daymet algorithm, spatially and temporally explicit empirical analyses of the relationships of temperature and precipitation to elevation are performed. | ||
Elev <- raster(Elevation) # elevation (meters) | ||
MaxTemp <- brick(MaxTemperature) # max temperature. degrees Celsius | ||
MinTemp <- brick(MinTemperature) # min temperature. degrees Celsius | ||
MeanTemp <- brick(MeanTemperature) # mean temperature. degrees Celsius | ||
Precip <- brick(Precipitation) # precipitation (mm) | ||
|
||
### ADDED vapor pressure - DayMet units for VP are Pascal - need to convert to [kPa] | ||
EaPascal <- brick(VP) # "daily average vapor pressure" - Vapor Pressure (Ea) | ||
Ea <- EaPascal/1000 | ||
|
||
MxRH <- brick(MaxRH) # max relative humidity. 0-100 (%) (verified) | ||
MnRH <- brick(MinRH) # min relative humidity. 0-100 (%) (verified) | ||
|
||
### Wind speed is commonly measured at 10 m height above the ground surface. | ||
U2 <- brick(Wind) # wind speed @ 2 meters height, assuming conversion from 10 meters. in m/s | ||
|
||
Sol_watts <- brick(SolarRad) # length(Sol_watts[is.na(Sol_watts)]) # [1] 55339 | ||
# Temp <- mask(Sol_watts, ETo); all.equal(Temp, Sol_watts); rm(Temp) | ||
Solar <- Sol_watts/11.6 # Rs - solar or shortwave radiation [MJ m-2 day-1] | ||
|
||
#### aeaMapping, Lat/Long rasters - see Appendix 1 | ||
Lat <- raster('Lat.values.grd') | ||
Long <- raster('Long.values.grd') | ||
Julian <- brick('Julian.values.grd') | ||
Rasters <- list(MaxTemp, MinTemp, MeanTemp, MxRH, MnRH, U2, Solar, Ea, Julian) | ||
print('do the rasters match?') | ||
print(sapply(Rasters, function(x) compareRaster(x, Elev))) | ||
|
||
# "Julian" is the rasterBrick with the values being the Julian day - see Appendix 2 | ||
|
||
##### III. Calculate Solar Parameters | ||
E0Max <- calc(MaxTemp, fun=function(x) {0.6108*exp(17.27*x/(x+273.3))}) | ||
Dr <- calc(Julian, fun=function(x) {1 + 0.033*cos(2*pi/365*x)}) # Inverse relative earth-sun distance | ||
Theta <- calc(Julian, fun=function(x) {0.409*sin((2*pi/365*x)-1.39)}) # Solar declination | ||
Lrad <- Lat*pi/180 # Latitude (radians) | ||
b <- calc(Julian, fun=function(x) {2*pi*(x-81)/364}) # | ||
Sc <- 0.1645*sin(2*b)-0.1255*cos(b)-0.025*sin(b) # seasonal correction for solar time [hour] | ||
Ws <- acos(-1*tan(Lrad)*tan(Theta)) # Sunset hour angle (radians) -- equation 25 | ||
|
||
Gsc <- 0.0820 # Solar constant (MJ/m^2*min) | ||
N <- 24/pi*Ws # Daylight hours (N) | ||
Multiply.Day <- calc(Dr, fun=function(x) {(24*60)/pi*Gsc*x}) # first part of equation 28 | ||
Ra <- Multiply.Day*(Ws*sin(Lrad)*sin(Theta)+cos(Lrad)*cos(Theta)*sin(Ws)) | ||
Rs0 <- (As0+Bs0*Elev)*Ra ###### Clear-sky radiation (Rso) | ||
|
||
#### Net longwave radiation - (detailed explanation on page 51): Rnl (MJ/meter^2*day) (Equation 39) | ||
Sigma <- 4.903*10^-9 # Stefan-Boltzmann constant [ 4.903 10^-9 MJ K^-4 m^-2 day^-1], | ||
KTmax <- MaxTemp + 273.16 # Max/Min daily temp in Kelvin: | ||
KTmin <- MinTemp + 273.16 | ||
|
||
###### Clear-sky radiation (Rs0) | ||
Bs0 <- Elev*10^-5 | ||
Rs0 <- (As0+Bs0*Elev)*Ra | ||
|
||
# IV. Calculate Intermediates: | ||
P <- calc(Elev, fun=function(x) {101.3*((293-0.0065*(x))/293)^5.26}) | ||
gamma <- Cp*P/E*LHV | ||
|
||
E0Max <- calc(MaxTemp, fun=function(x) {0.6108*exp(17.27*x/(x+273.3))}) # as does this one | ||
E0Min <- calc(MinTemp, fun=function(x) {0.6108*exp(17.27*x/(x+273.3))}) | ||
Es <- (E0Min+E0Max)/2 | ||
E0Mean <- calc(MeanTemp, fun=function(x) {0.6108*exp(17.27*x/(x+273.3))}) | ||
Delta <- 4098*(E0Mean)/((MeanTemp+273.3)^2) | ||
|
||
### Use Raster of VP - need units of kPa | ||
Ea1 <- (E0Min*MxRH/100+E0Max*MnRH/100)/2 # Actual Vapor Pressure [kPa] | ||
# writeRaster(Ea1, 'Derived.VP.kPa.NCAR.min.max.RH.grd', overwrite = TRUE) | ||
VPD <- Es - Ea # Vapor Pressure Deficit [kPa] | ||
|
||
#### Net longwave radiation - (detailed explanation on page 51): Rnl (MJ/meter^2*day) (Equation 39) | ||
Sigma <- 4.903*10^-9 # Stefan-Boltzmann constant [ 4.903 10^-9 MJ K^-4 m^-2 day^-1], | ||
KTmax <- MaxTemp + 273.16 # Max/Min daily temp in Kelvin: | ||
KTmin <- MinTemp + 273.16 | ||
Rnl <- Sigma*((KTmax^4+KTmin^4)/2)*(0.34-0.14*Ea^0.5)*(1.35*(Solar/Rs0)-0.35) | ||
# Final product: Rnl (net longwave Radiation) | ||
|
||
# Solar radiation (Rs) | ||
Rns <- calc(Solar, fun=function(x) {(1-Alpha)*x}) ##### Net solar or net shortwave radiation: Rns (MJ/meter^2*day) | ||
Rn <- Rns - Rnl # Net longwave Radiation | ||
Rn <- dropLayer(Rn, 366) | ||
# writeRaster(Rn, 'Rn.2008.grd', overwrite = TRUE) | ||
|
||
### V. Final ETo calculation steps: | ||
Numerator <- 0.408*Delta*(Rn-G)+gamma*900/(MeanTemp+273)*U2*(Es-Ea) | ||
Denominator <- Delta+gamma*(1+0.34*U2) | ||
ETo <- Numerator/Denominator | ||
|
||
####### 13. Clean out negative ETo values: | ||
ETo[ETo < 0] <- 0.001 | ||
ETo <- mask(ETo, U2) | ||
writeRaster(ETo, filename = Filename, overwrite = TRUE) | ||
# YearAve.ETo <- calc(s, fun = mean, filename = paste0('Annual.average.ETo_2008.grd'), overwrite = TRUE) | ||
# YearTotal.ETo <- calc(s, fun = sum, filename = paste0('Year.total.ETo_2008.grd'), overwrite = TRUE) | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,76 @@ | ||
Daily.Crop.Curves <- | ||
function(Croplayer, StateNames, Stages, Kcb_tab, MaxHeight){ | ||
##### INPUTS: | ||
### Stages: a vector of length four, giving the whole number length in days of the following four seasons | ||
### 4 stages: initial, crop development, mid season, late season | ||
### Kcb_tab: a vector of length three, the Kcb for a given crop | ||
### StateNames: a vector of length of resultant list, same order as the nrows, of State IDs | ||
|
||
#### OUTPUTS: a data.frame with (x) the length of the growing season (in days), and the following variables: | ||
# DailyKcb$Season_day : # of day in the growing season | ||
# DailyKcb$Stage_ID : factor (1:4); stage of growing season | ||
# DailyKcb$daily_Kcb : daily Kcb values | ||
# DailyKcb$daily_Kcb : day_height - daily interpolated height of crop | ||
# ALSO: returns the exposed soil fraction (1-Fc), over the lenght of the growing season | ||
|
||
## Added Croplayer, need to access it to get rooting depths | ||
root.depth <- read.csv('crop.roots.csv') | ||
Crop <- Croplayer | ||
if (Croplayer == 'spring_barley' | Croplayer == 'fall_barley'){ | ||
Crop <- 'barley' | ||
} | ||
if (Croplayer == 'spring_oats' | Croplayer == 'fall_oats'){ | ||
Crop <- 'oats' | ||
} | ||
if (Croplayer == 'durum_wheat'){ | ||
Crop <- 'spring_wheat' | ||
} | ||
|
||
root.depth <- subset(root.depth, crop == Crop, select = c(min_rooting_depth, max_rooting_depth)) | ||
Season <- rowSums(Stages) | ||
|
||
Base <- lapply(Season, function(x) x = c(1:x)) | ||
Next <- Base; Kci <- Base; Height.Kci <- Base; DailyKcb <- Base; Growth_split <- Base; season.ID <- Base; day_height <- Base; MAX.Height <- Base; One.Minus.Fc <- Base # Initialize vectors | ||
Roots <- Base; Day_Roots <- Base | ||
|
||
for (i in 1:length(Next)){ | ||
|
||
B <- c(Stages[i,1], sum(c(Stages[i,1], Stages[i,2])), sum(c(Stages[i,1], Stages[i,2], Stages[i,3])), sum(Stages[i,])) | ||
|
||
Xs <- c(1, B[1], mean(c(B[1], B[1], B[2])), mean(c(B[1], B[2])), mean(c(B[2], B[2], B[3])), B[3], mean(c(B[3], B[4])), B[4]) | ||
|
||
Ys <- c(0.01, 0.02, Kcb_tab[[1]], Kcb_tab[[2]], sum(c(Kcb_tab[[1]], Kcb_tab[[2]], Kcb_tab[[3]])), Kcb_tab[[2]], mean(c(Kcb_tab[[2]], Kcb_tab[[3]])), Kcb_tab[[3]]) | ||
|
||
P <- data.frame(bezierCurve(Xs,Ys,500)) | ||
Q <- P[!duplicated(round(P$x)),] | ||
Q$x <- round(Q$x) | ||
|
||
Kci[[i]] <- round(Q$y, digits = 2) | ||
|
||
## Height calculation | ||
Growth_split[[i]] <- unlist(c(0.01, .015, rep(NA, times = ceiling(Stages[i,1]/2)-2), Kcb_tab[1], rep(NA, times = floor(Stages[i,1]/2)-1), | ||
rep(NA, times = Stages[i,2]-1), Kcb_tab[2], | ||
rep(Kcb_tab[2], times = Stages[i,3]), | ||
rep(NA, times = (Stages[i,4]-1)), Kcb_tab[3])) | ||
Height.Kci[[i]] <- spline(Base[[i]], Growth_split[[i]], xout = Base[[i]], method = 'natural', ties = mean)$y | ||
MAX.Height[[i]] <- max(Height.Kci[[i]]) | ||
day_height[[i]] <- round(Height.Kci[[i]]*MaxHeight/MAX.Height[[i]], 2) | ||
|
||
## Root growth: (currently goes from min rooting to max, and then halfway back in the post-season) | ||
Zs <- c(0.1, 0.15, mean(c(0.15, root.depth[[1]])), root.depth[[1]], root.depth[[2]], sum(c(mean(c(0.15, root.depth[[1]])), root.depth[[2]])), root.depth[[2]], root.depth[[1]]) | ||
R <- data.frame(bezierCurve(Xs,Zs,500)) | ||
S <- R[!duplicated(round(R$x)),] | ||
S$x <- round(S$x) | ||
|
||
Day_Roots[[i]] <- round(S$y, digits = 2) | ||
|
||
### Season.ID | ||
season.ID[[i]] <- as.factor(c(rep(1, times = Stages[i,1]), rep(2, times = Stages[i,2]), | ||
rep(3, times = Stages[i,3]), rep(4, times = Stages[i,4]))) | ||
DailyKcb[[i]] <- as.data.frame(cbind(Base[[i]], round(Kci[[i]], 2), round(Day_Roots[[i]], 2), day_height[[i]], as.factor(season.ID[[i]]))) | ||
names(DailyKcb[[i]]) <- c('Season_day', 'daily_Kcb', "daily_root.depth", "day_height", "season.ID") | ||
} | ||
names(DailyKcb) <- StateNames | ||
save(DailyKcb, file = paste0(Intermediates, 'Daily.Crop.Profile.', Croplayer, '.Rdata')) | ||
return(DailyKcb) | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,8 @@ | ||
Daily.Crop.Parameters <- | ||
function(Croplayer){ | ||
Calc.Basal.Crop.Coeff(Croplayer) | ||
KcMAX(Croplayer) | ||
KcMAX.fallow(Croplayer) | ||
Calc.Fc.Few(Croplayer) | ||
Fallow.Few.Calc(Croplayer) | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,35 @@ | ||
Estimate.Crop <- | ||
function(crop){ | ||
|
||
Subset.Growth.Season(crop) | ||
print(paste('Seasons subsetted and rescaled for', crop)) | ||
Daily.Crop.Parameters(crop) | ||
print(paste('Daily crop parameters estimated for', crop)) | ||
Final.Daily.ET.Calc(crop) | ||
print(paste('Daily crop water balance estimated for', crop)) | ||
|
||
Sum.Save.Daily.Evapotranspiration(crop, rainfed = TRUE) | ||
Sum.Save.Daily.Evapotranspiration(crop, rainfed = FALSE) | ||
print(paste('Daily ET rasters made for', crop)) | ||
|
||
Sum.Save.Water.Balances(crop, rainfed = FALSE, type = 'seasonal', BW.GW = FALSE) | ||
Sum.Save.Water.Balances(crop, rainfed = FALSE, type = 'annual', BW.GW = FALSE) | ||
Sum.Save.Water.Balances(crop, rainfed = TRUE, type = 'seasonal', BW.GW = FALSE) | ||
Sum.Save.Water.Balances(crop, rainfed = TRUE, type = 'annual', BW.GW = FALSE) | ||
Sum.Save.Water.Balances(crop, rainfed = FALSE, type = 'seasonal', BW.GW = TRUE) | ||
|
||
Generate.Land.Use(crop) | ||
print(paste('land use raster generated for', crop)) | ||
|
||
print(paste('Annual and seasonal water balance rasters saved for', crop)) | ||
|
||
SuperImpose.WB.on.LU(crop, rainfed = FALSE, type = 'seasonal', Growing.Season.GW.BW = FALSE) | ||
SuperImpose.WB.on.LU(crop, rainfed = FALSE, type = 'annual', Growing.Season.GW.BW = FALSE) | ||
SuperImpose.WB.on.LU(crop, rainfed = TRUE, type = 'seasonal', Growing.Season.GW.BW = FALSE) | ||
SuperImpose.WB.on.LU(crop, rainfed = TRUE, type = 'annual', Growing.Season.GW.BW = FALSE) | ||
SuperImpose.WB.on.LU(crop, rainfed = FALSE, type = 'seasonal', Growing.Season.GW.BW = TRUE) | ||
|
||
print(paste('Water balances superimposed on land use for', crop)) | ||
|
||
|
||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,34 @@ | ||
Fallow.Few.Calc <- | ||
function(Croplayer){ | ||
# Fc - the effective fraction of soil surface covered by vegetation [0 - 0.99], | ||
load(paste0(Intermediates, paste('Fallow.Season', Croplayer, 'MNRH_', 'Rdata', sep = '.'))); Template <- Fallow.Season; rm(Fallow.Season) | ||
Template <- lapply(Template, function(x) x[,(grep('layer', names(x)))]); Fc <- Template | ||
|
||
load(paste0(Intermediates, paste('KcMax.Fallow', Croplayer, 'Rdata', sep = '.'))) | ||
KcMax.fallow <- lapply(KcMax, function(x) x[,(grep('layer', names(x)))]) | ||
all.equal(sapply(Template, dim), sapply(KcMax.fallow, dim)) | ||
|
||
Off.season.vars <- c('winter_wheat', 'durum_wheat', 'fall_barley', 'fall_oats') | ||
|
||
if (Croplayer %in% Off.season.vars){ | ||
KcMin <- lapply(Template, function(x) c(rep(.15, times = (length(x))))) | ||
DayHeight <- lapply(Template, function(x) c(rep(.15, times = (length(x))))) | ||
Kcb <- lapply(Template, function(x) c(rep(1, times = (length(x))))) | ||
} | ||
if (!(Croplayer %in% Off.season.vars)){ | ||
KcMin <- lapply(Template, function(x) c(rep(.03, times = (length(x))))) | ||
DayHeight <- lapply(Template, function(x) c(rep(.05, times = (length(x))))) | ||
Kcb <- lapply(Template, function(x) c(rep(.07, times = (length(x))))) | ||
} | ||
|
||
# GR.KcMax <- lapply(Template, function(x) c(rep(.16, times = (length(x))))) | ||
|
||
for (i in 1:length(Fc)){ | ||
for (j in 1:length(DayHeight[[i]])){ | ||
Fc[[i]][,j] <- ((Kcb[[i]][j] - KcMin[[i]][j])/(KcMax.fallow[[i]][j] - KcMin[[i]][j]))^(1+0.5*DayHeight[[i]][j]) | ||
} | ||
} | ||
Fallow.Few <- Fc | ||
Fallow.Few <- lapply(Fc, function(x) 1-x[]) # This suffices for RAINFED crops | ||
save(Fallow.Few, file = paste0(Intermediates, paste('Fallow.Few', Croplayer, 'Rdata', sep = '.'))) | ||
} |
Oops, something went wrong.