Compare NCEP CRUNCEP
David LeBauer
4/27/2015
Get Climate data from CIAT in Columbia
- BETYdb record: https://www.betydb.org/sites/1232
- lat: 3.0127854
- lon: -76.5006448
Extract data subset:
qsub -I
cd /home/groups/ebimodeling/met/cruncep/
module load nco
ncks -O -d lon,-76.75,-76.25 -d lat,2.75,3.25 out/all.nc ciat.nc
ncrename -d lat,latitude -d lon,longitude ciat.nc
#install_github("pecanproject/pecan/modules/data.atmosphere")
library(PEcAn.data.atmosphere)
#install_github("pecanproject/pecan/models/biocro")
library(PEcAn.BIOCRO)##
## Attaching package: 'PEcAn.BIOCRO'
##
## The following objects are masked _by_ '.GlobalEnv':
##
## cf2biocro, met2model.BIOCROlibrary(data.table)
library(ncdf4)
library(lubridate)##
## Attaching package: 'lubridate'
##
## The following objects are masked from 'package:data.table':
##
## hour, mday, month, quarter, wday, week, yday, yearlibrary(udunits2)
ciat.nc <- nc_open("raw_data/ciat.nc")
ciat.cfmet <- load.cfmet(met.nc = ciat.nc, lat = 3.01, lon = -76.5, start.date = "1983-01-01", end.date = "1983-12-31")## Warning: All formats failed to parse. No formats found.## Warning: All formats failed to parse. No formats found.ciat.hourly <- cfmet.downscale.time(ciat.cfmet)
ciat.biocromet <- cf2biocro(ciat.hourly, longitude = -76.5, zulu2solarnoon = TRUE)
save(ciat.biocromet, file = "./refined_data/ciat.met.Rdata")Compare NCEP and CRUNCEP met
library(data.table)
library(ggplot2)
theme_set(theme_bw())
library(BioCro)## Loading required package: lattice
## Loading required package: coda
## Loading required package: knitr
## Loading required package: roxygen2
##
## Attaching package: 'BioCro'
##
## The following object is masked from 'package:PEcAn.data.atmosphere':
##
## lightMEload("./refined_data/ICTP.hourly.climate.Rd")
load("refined_data/ciat.met.Rdata")
head(ICTP.hourly.climate)## year doy hour SolarR Temp RH WS precip
## [1,] 1983 1 0 0 22.80000 0.9066987 0.3398889 0
## [2,] 1983 1 1 0 21.64020 0.9146447 0.3398889 0
## [3,] 1983 1 2 0 20.75026 0.9250000 0.3398889 0
## [4,] 1983 1 3 0 20.19082 0.9370590 0.3398889 0
## [5,] 1983 1 4 0 20.00000 0.9500000 0.3398889 0
## [6,] 1983 1 5 0 20.19082 0.9629410 0.3398889 0head(ciat.biocromet)## year doy hour SolarR Temp RH WS precip
## 1 1982 365 19 0.0000000 22.83648 0.7663728 0 0.02937205
## 2 1982 365 20 0.0000000 22.73468 0.7657333 0 0.03307075
## 3 1982 365 21 0.0000000 22.63289 0.7650594 0 0.03676944
## 4 1982 365 22 0.0000000 22.53109 0.7643507 0 0.04046813
## 5 1982 365 23 0.1930545 22.37245 0.7634611 0 0.04326784
## 6 1982 365 24 0.7472740 22.14559 0.7622690 0 0.04498876ICTP <- data.table(ICTP.hourly.climate)
ICTP <- ICTP[, lapply(.SD, as.numeric)]
met <- rbindlist(list(cbind(source = "NCEP", ICTP),
cbind(source = "CRUNCEP", ciat.biocromet)),
use.names = TRUE, fill = FALSE)
met <- cbind(met[,list(day = doy+hour/24)], met)ggplot() + geom_line(data = met, aes(day, SolarR, color = source)) + xlim(c(180,210))## Warning: Removed 16054 rows containing missing values (geom_path).
ggplot() + geom_line(data = met, aes(day, RH, color = source)) + xlim(c(180,210))## Warning: Removed 16054 rows containing missing values (geom_path).
ggplot() + geom_line(data = met, aes(day, Temp, color = source)) + xlim(c(180,210))## Warning: Removed 16054 rows containing missing values (geom_path).
ggplot() + geom_line(data = met, aes(day, precip, color = source)) + xlim(c(180,210))## Warning: Removed 16054 rows containing missing values (geom_path).
ggplot() + geom_line(data = met, aes(day, WS, color = source)) + xlim(c(180,210))## Warning: Removed 16054 rows containing missing values (geom_path).
Effect of climate data source on yield and partitioning
NCEP
Parms (from cassava.Rmd)
doy_break <-c(2, 4, 5, 6, 8, 12) * 30 # assumign 30 days in a month
ncep <- as.data.frame(met[source == "NCEP"])
ncep$source <- NULL
ncep$day <- NULL
TT <- cumsum(ncep$Tavg)
thermal_break <- c(TT[60], TT[120], TT[150], TT[180], TT[240], TT[360])
cassavaphenology <- phenoParms(
tp1 = thermal_break[1], tp2 = thermal_break[2],
tp3 = thermal_break[3], tp4 = thermal_break[4],
tp5 = thermal_break[5], tp6 = thermal_break[6],
kStem1 = 0.495, kLeaf1 = 0.495, kRoot1 = 0.01, kRhizome1 = 0.0,
kStem2 = 0.30, kLeaf2 = 0.40, kRoot2 = 0.015, kRhizome2 = 0.285,
kStem3 = 0.20, kLeaf3 = 0.20, kRoot3 = 0.03, kRhizome3 = 0.57,
kStem4 = 0.25, kLeaf4 = 0.11, kRoot4= 0.04, kRhizome4 = 0.60,
kStem5 = 0.25, kLeaf5 = 0.05, kRoot5 = 0.03, kRhizome5 = 0.67,
kStem6 = 0.05, kLeaf6 = 0.05, kRoot6 = 0.045, kRhizome6 = 0.855)
cassava_senescence <- willowseneParms(senLeaf = 4000,
senStem = 10000,
senRoot = 10000,
senRhizome = 10000,
leafdeathrate = 0.2)
cassava_initial <- list(iRhizome = 1.0, iStem = 0.0)
cassava_photo <- list(Rd = 1.1)
ncep.cassava <- willowGro(WetDat = ncep,
day1 = 5, dayn=360, lat = 3.5,
willowphenoControl = cassavaphenology,
seneControl=cassava_senescence,
iPlantControl = cassava_initial,
photoControl = cassava_photo)
plot(ncep.cassava)CRUNCEP
cruncep <- as.data.frame(met[source == "CRUNCEP"])
cruncep$day <- cruncep$source <- NULL
TT <- cumsum(cruncep$Tavg)
thermal_break <- c(TT[60], TT[120], TT[150], TT[180], TT[240], TT[360])
cassavaphenology <- phenoParms(
tp1 = thermal_break[1], tp2 = thermal_break[2],
tp3 = thermal_break[3], tp4 = thermal_break[4],
tp5 = thermal_break[5], tp6 = thermal_break[6],
kStem1 = 0.495, kLeaf1 = 0.495, kRoot1 = 0.01, kRhizome1 = 0.0,
kStem2 = 0.30, kLeaf2 = 0.40, kRoot2 = 0.015, kRhizome2 = 0.285,
kStem3 = 0.20, kLeaf3 = 0.20, kRoot3 = 0.03, kRhizome3 = 0.57,
kStem4 = 0.25, kLeaf4 = 0.11, kRoot4= 0.04, kRhizome4 = 0.60,
kStem5 = 0.25, kLeaf5 = 0.05, kRoot5 = 0.03, kRhizome5 = 0.67,
kStem6 = 0.05, kLeaf6 = 0.05, kRoot6 = 0.045, kRhizome6 = 0.855)
cruncep.cassava <- willowGro(WetDat = cruncep,
day1 = 5, dayn=360, lat = 3.5,
willowphenoControl = cassavaphenology,
seneControl=cassava_senescence,
iPlantControl = cassava_initial,
photoControl = cassava_photo)
plot(cruncep.cassava)
plot(cruncep.cassava$Stem, type = 'l')
lines(ncep.cassava$Stem, col = 'red')## Warning: Removed 16054 rows containing missing values (geom_point).## Warning: Removed 16054 rows containing missing values (geom_path).## Warning: Removed 16054 rows containing missing values (geom_point).## Warning: Removed 16054 rows containing missing values (geom_path).## Warning: Removed 16054 rows containing missing values (geom_point).## Warning: Removed 16054 rows containing missing values (geom_path).## Warning: Removed 16054 rows containing missing values (geom_point).## Warning: Removed 16054 rows containing missing values (geom_path).## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.## Warning: Removed 8039 rows containing missing values (stat_smooth).## Warning: Removed 8015 rows containing missing values (stat_smooth).## Warning: Removed 16054 rows containing missing values (geom_point).
Compare with site data
ciat.nc <- nc_open("raw_data/ciat.nc")
ciat.cfmet <- load.cfmet(met.nc = ciat.nc, lat = 3.01, lon = -76.5, start.date = "1991-05-01", end.date = "1993-4-31")## Warning: All formats failed to parse. No formats found.## Warning: All formats failed to parse. No formats found.biocro.met <- cf2biocro(ciat.cfmet, longitude = -76.5, zulu2solarnoon = TRUE)
ciat.cfmet <- cbind(ciat.cfmet, SolarR = biocro.met$SolarR)
cruncep_monthly <- ciat.cfmet[,list(rainfall = round(sum(precipitation_flux) * 60 * 60 * 3),
solar = round(mean(SolarR)),#ud.convert(round(sum(surface_downwelling_shortwave_flux_in_air)), "W*3h/m2", "MJ/m2"),
Tmax = round(max(air_temperature) -273.15, 1),
Tmin = round(min(air_temperature) -273, 1),
Tmean = round(mean(air_temperature) -273, 1)),
by = 'year,month']
santander <- fread("raw_data/santander_monthly_met.csv")
santander$pan <- NULL
met <- rbindlist(list(cbind(source = 'station', santander), cbind(source = 'cruncep', cruncep_monthly)), use.names = FALSE, fill = FALSE)
met[,`:=`(date = lubridate::ymd(paste(year, month, '01')))]## source year month rainfall solar tmax tmin tmean date
## 1: station 1991 5 147 493 29.2 19.5 24.1 1991-05-01
## 2: station 1991 6 52 521 29.9 19.2 24.5 1991-06-01
## 3: station 1991 7 38 530 29.9 18.1 24.1 1991-07-01
## 4: station 1991 8 80 528 30.4 16.4 24.3 1991-08-01
## 5: station 1991 9 152 511 30.8 18.0 24.6 1991-09-01
## 6: station 1991 10 161 544 29.9 17.8 23.7 1991-10-01
## 7: station 1991 11 190 417 29.0 18.6 23.4 1991-11-01
## 8: station 1991 12 116 513 29.8 18.5 24.4 1991-12-01
## 9: station 1992 1 66 537 30.0 18.2 24.6 1992-01-01
## 10: station 1992 2 181 449 30.3 18.8 24.4 1992-02-01
## 11: station 1992 3 68 498 31.3 19.2 25.2 1992-03-01
## 12: station 1992 4 194 459 30.0 19.0 24.6 1992-04-01
## 13: station 1992 5 174 446 29.6 19.3 24.4 1992-05-01
## 14: station 1992 6 41 479 30.2 18.4 25.1 1992-06-01
## 15: station 1992 7 56 522 29.8 17.0 24.5 1992-07-01
## 16: station 1992 8 24 465 31.0 17.4 25.2 1992-08-01
## 17: station 1992 9 98 501 30.0 18.0 24.4 1992-09-01
## 18: station 1992 10 77 496 29.7 18.1 24.3 1992-10-01
## 19: station 1992 11 134 430 29.1 18.1 23.4 1992-11-01
## 20: station 1992 12 129 478 29.5 18.9 24.1 1992-12-01
## 21: station 1993 1 131 533 29.6 18.1 24.1 1993-01-01
## 22: station 1993 2 86 441 29.8 18.1 24.0 1993-02-01
## 23: station 1993 3 121 455 29.8 18.4 23.6 1993-03-01
## 24: station 1993 4 252 460 29.4 18.8 23.7 1993-04-01
## 25: cruncep 1991 5 115 423 25.8 20.6 22.7 1991-05-01
## 26: cruncep 1991 6 25 407 26.3 21.0 23.1 1991-06-01
## 27: cruncep 1991 7 41 431 25.4 20.2 22.6 1991-07-01
## 28: cruncep 1991 8 63 446 25.7 19.7 22.5 1991-08-01
## 29: cruncep 1991 9 92 453 26.2 20.9 23.0 1991-09-01
## 30: cruncep 1991 10 65 443 26.6 20.9 23.1 1991-10-01
## 31: cruncep 1991 11 110 422 25.2 20.4 22.5 1991-11-01
## 32: cruncep 1991 12 64 384 26.6 20.8 23.7 1991-12-01
## 33: cruncep 1992 1 25 424 26.9 20.4 23.4 1992-01-01
## 34: cruncep 1992 2 44 461 27.5 20.8 24.0 1992-02-01
## 35: cruncep 1992 3 39 470 27.7 21.9 24.5 1992-03-01
## 36: cruncep 1992 4 68 446 26.9 21.8 24.0 1992-04-01
## 37: cruncep 1992 5 52 429 25.6 20.9 22.9 1992-05-01
## 38: cruncep 1992 6 21 421 26.8 21.7 24.1 1992-06-01
## 39: cruncep 1992 7 84 419 26.3 20.5 23.2 1992-07-01
## 40: cruncep 1992 8 25 424 27.2 21.1 23.6 1992-08-01
## 41: cruncep 1992 9 168 453 26.4 21.1 23.1 1992-09-01
## 42: cruncep 1992 10 34 442 25.5 20.8 22.7 1992-10-01
## 43: cruncep 1992 11 123 425 25.5 20.8 22.9 1992-11-01
## 44: cruncep 1992 12 94 402 26.0 20.7 23.2 1992-12-01
## 45: cruncep 1993 1 72 437 26.4 20.2 23.1 1993-01-01
## 46: cruncep 1993 2 105 462 26.4 20.7 23.3 1993-02-01
## 47: cruncep 1993 3 83 467 25.6 19.3 22.3 1993-03-01
## 48: cruncep 1993 4 72 447 25.7 20.4 22.8 1993-04-01
## source year month rainfall solar tmax tmin tmean datelibrary(ggthemes)
theme_set(theme_linedraw())
ggplot(data = met,aes(date, rainfall, color = source)) + geom_line() + geom_point(size = 3) + ylim(c(0,200)) + ggtitle("Monthly Rainfall (mm)")## Warning: Removed 1 rows containing missing values (geom_path).## Warning: Removed 1 rows containing missing values (geom_point).
ggplot(data = met,aes(date, solar, color = source)) + geom_line() + geom_point(size = 3) + ylim(c(0, 600)) + ggtitle("Monthly Solar Radiation (W/m2)")
ggplot(data = met) + geom_line(aes(date, tmax, color = source), linetype = 2) +
geom_point(aes(date, tmax, color = source), alpha = 0.5, size = 3 )+
geom_line(aes(date, y = tmean, color = source)) +
geom_point(aes(date, y = tmean, color = source), size = 3) +
geom_line(aes(date, y = tmin, color = source), linetype = 3) +
geom_point(aes(date, y = tmin, color = source), alpha = 0.5, size = 3) +
xlab("Temp (C)") + ggtitle("Monthly mean, min, max temperature")
library(plotrix)
model <- met[source == 'cruncep']
ref <- met[source == 'station']
taylor.diagram(ref = ref$rainfall, model = model$rainfall, main = "Montly Precip (mm/month)", ref.sd = TRUE)
taylor.diagram(ref = ref$tmean, model = model$tmean, col = "blue", ref.sd = TRUE,
main = "Monthly Temp:\nmean(blue), min(light blue)m max(dark blue))")
taylor.diagram(ref = ref$tmax, model = model$tmax, add = TRUE, col = "darkblue")
taylor.diagram(ref = ref$tmin, model = model$tmin, add = TRUE, col = "lightblue")
ans <- rbindlist(list(
ref[,list(tmin = mean(tmin), tmax = mean(tmax), tmean = mean(tmean), precip = mean(precip), solar = mean(solar))],
model[,list(tmin = mean(tmin), tmax = mean(tmax), tmean = mean(tmean), precip = mean(precip), solar = mean(solar))],
data.table(tmin = mean(ref$tmin - model$tmin), tmax = mean(ref$tmax - model$tmax), tmean = mean(ref$tmean - model$tmean), precip = mean(ref$rainfall - model$rainfall), solar = mean(ref$solar - model$solar))))
ans2 <- round(t(as.matrix(ans)), 1)
colnames(ans2) <- c("station", "cruncep", "difference")
kable(data.frame(ans2))| station | cruncep | difference | |
|---|---|---|---|
| tmin | 18.3 | 20.7 | -2.4 |
| tmax | 29.9 | 26.3 | 3.7 |
| tmean | 24.3 | 23.2 | 1.1 |
| precip | 34.9 | 34.9 | 45.2 |
| solar | 487.8 | 434.9 | 52.8 |