Skip to content
Snippets Groups Projects
Commit c9f53b38 authored by og's avatar og
Browse files

lintr formatting

parent de10bb93
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
testcases/estuary/plot_output.R
+ 73
− 63
View file @ c9f53b38
......@@ -3,108 +3,118 @@
# SPDX-FileContributor Ovidio Garcia-Oliva <ovidio.garcia@hereon.de>
if (!require("ncdf4")) install.packages("ncdf4")
is.oxypom_ = T # True if OxyPOM is validated False if DiaMo is validated
################################################################################
## definition of local functions
Get.Year = function(...) format(as.Date(..., format="%d/%m/%Y"),"%Y")
as.POSIXct = function(...) base::as.POSIXct(...,,format="%Y-%m-%d %H:%M:%S")
get.year = function(...) format(as.Date(..., format = "%d/%m/%Y"), "%Y")
as.POSIXct = function(...) base::as.POSIXct(..., format = "%Y-%m-%d %H:%M:%S")
################################################################################
## loading the results of gotm
nc_data = nc_open('output.nc')
if("diamo_PAR"%in%names(nc_data$var)) is.oxypom_ = F # DiaMo is validated
nc_data = nc_open("output.nc")
temp = t(ncvar_get(nc_data, "temp"))
is_oxypom = TRUE # True if OxyPOM is validated False if DiaMo is validated
if("diamo_PAR"%in%names(nc_data$var)) is_oxypom = FALSE # DiaMo is validated
if(is.oxypom_){
oxy = t(ncvar_get(nc_data, "oxypom_DOxy"))
temp = t(ncvar_get(nc_data, "temp"))
if(is_oxypom){
oxy = t(ncvar_get(nc_data, "oxypom_DOxy"))
}else{
oxy = t(ncvar_get(nc_data, "diamo_OXY"))
oxy = t(ncvar_get(nc_data, "diamo_OXY"))
}
################################################################################
## loading observed values for temperature (temp) and dissolved oxygen (DO)
temp.obs.1 = read.delim("./data/FGG_Elbe_008!Wassertemperatur.txt", header=FALSE, comment.char="#")
temp.obs.2 = read.delim("./data/FGG_Elbe_009!Wassertemperatur.txt", header=FALSE, comment.char="#")
temp.obs.3 = read.delim("./data/LZ_AL!Wassertemperatur.txt", header=FALSE, comment.char="#")
temp_obs_1 = read.delim("./data/FGG_Elbe_008!Wassertemperatur.txt", header = FALSE, comment.char = "#")
temp_obs_2 = read.delim("./data/FGG_Elbe_009!Wassertemperatur.txt", header = FALSE, comment.char = "#")
temp_obs_3 = read.delim("./data/LZ_AL!Wassertemperatur.txt", header = FALSE, comment.char = "#")
temp.obs.1$V1 = as.POSIXct(temp.obs.1$V1)
temp.obs.2$V1 = as.POSIXct(temp.obs.2$V1)
temp.obs.3$V1 = as.POSIXct(temp.obs.3$V1)
temp_obs_1$V1 = as.POSIXct(temp_obs_1$V1)
temp_obs_2$V1 = as.POSIXct(temp_obs_2$V1)
temp_obs_3$V1 = as.POSIXct(temp_obs_3$V1)
temp.obs.1 = subset(temp.obs.1, Get.Year(V1)>=2004)
temp.obs.2 = subset(temp.obs.2, Get.Year(V1)>=2004)
temp.obs.3 = subset(temp.obs.3, Get.Year(V1)>=2004)
temp_obs_1 = subset(temp_obs_1, get.year(V1) >= 2004)
temp_obs_2 = subset(temp_obs_2, get.year(V1) >= 2004)
temp_obs_3 = subset(temp_obs_3, get.year(V1) >= 2004)
DO.obs.1 = read.delim("./data/FGG_Elbe_008!Sauerstoffgehalt_(Einzelmessung).txt", header=FALSE, comment.char="#")
DO.obs.2 = read.delim("./data/FGG_Elbe_009!Sauerstoffgehalt_(Einzelmessung).txt", header=FALSE, comment.char="#")
DO_obs_1 = read.delim("./data/FGG_Elbe_008!Sauerstoffgehalt_(Einzelmessung).txt", header = FALSE, comment.char = "#")
DO_obs_2 = read.delim("./data/FGG_Elbe_009!Sauerstoffgehalt_(Einzelmessung).txt", header = FALSE, comment.char = "#")
DO.obs.1$V1 = as.POSIXct(DO.obs.1$V1)
DO.obs.2$V1 = as.POSIXct(DO.obs.2$V1)
DO_obs_1$V1 = as.POSIXct(DO_obs_1$V1)
DO_obs_2$V1 = as.POSIXct(DO_obs_2$V1)
DO.obs.1 = subset(DO.obs.1, Get.Year(V1)>=2004)
DO.obs.2 = subset(DO.obs.2, Get.Year(V1)>=2004)
DO_obs_1 = subset(DO_obs_1, get.year(V1) >= 2004)
DO_obs_2 = subset(DO_obs_2, get.year(V1) >= 2004)
################################################################################
## plotting the data
x = 1:length(temp[,1])
x = as.Date(x,origin='2006-01-01')
x = as.Date(x, origin = "2006-01-01")
surface = length(temp[1,])-1 # layer number of the surface
surface = length(temp[1,]) - 1 # layer number of the surface
bottom = 1 # layer number of the bottom
col.sim='tomato'
col.sim.bottom = 'pink'
col.obs='black'
col.obs2='gray50'
col.obs3='gray50'
col_sim = "tomato"
col_sim_bottom = "pink"
col_obs = "black"
col_obs2 = "gray50"
col_obs3 = "gray50"
png(filename="./estuary_validation.png",
png(filename = "./estuary_validation.png",
width = 800, height = 400)
par(mfrow=c(2,1),mai=c(0.42,0.42,0.21,0.21),oma=2*c(1,1,0.5,0.5),las=1)
par(mfrow = c(2,1), mai = c(0.42,0.42,0.21,0.21), oma = 2 * c(1, 1, 0.5, 0.5), las = 1)
plot(x,rowMeans(temp[,bottom:surface]),type='n',
col=NA,ylim=c(-5,25),log='',lwd=0.5,
ylab='',
xlab='days'
plot(x, rowMeans(temp[,bottom:surface]),
type = "n",
col = NA,
ylim = c(-5, 25),
log = "",
lwd = 0.5,
ylab = "",
xlab = "days"
)
title('temperature at the surface degC',adj=0,line=0.1,cex=0.5,font.main=1)
points(as.Date(temp.obs.2$V1),temp.obs.2$V2,col=col.obs2,pch=0,cex=.5)
lines(as.Date(temp.obs.3$V1),temp.obs.3$V2,col=col.obs3,pch=20,cex=1)
points(as.Date(temp.obs.1$V1),temp.obs.1$V2,col=col.obs,pch=20,cex=1)
lines(x,temp[,bottom],lty=1,col=col.sim.bottom)
lines(x,temp[,surface],lty=1,col=col.sim,lwd=2)
title("temperature at the surface degC", adj = 0, line = 0.1, cex = 0.5, font.main = 1)
points(as.Date(temp_obs_2$V1), temp_obs_2$V2, col = col_obs2, pch = 0,cex = .5)
lines(as.Date(temp_obs_3$V1), temp_obs_3$V2, col = col_obs3, pch = 20,cex = 1)
points(as.Date(temp_obs_1$V1), temp_obs_1$V2, col = col_obs, pch = 20,cex = 1)
lines(x, temp[,bottom], lty = 1, col = col_sim_bottom)
lines(x, temp[,surface], lty = 1, col = col_sim, lwd=2)
cf = 1000/32 # conversion factor of mg L-1 to mmol-O2 L
plot(x,rowMeans(oxy[,1:2]),type='n',
col=NA,ylim=c(100,500),log='',lwd=0.5,
ylab='',
xlab='days',
plot(x, rowMeans(oxy[,1:2]),
type = "n",
col = NA,ylim = c(100, 500),
log = "",
lwd = 0.5,
ylab = "",
xlab = "days",
)
title('dissolved oxygen concentration mmol-O2 L-1',adj=0,line=0.1,cex=0.5,font.main=1)
points(as.Date(DO.obs.2$V1),(1000/32)*DO.obs.2$V2,col=col.obs2,pch=0,cex=.5)
points(as.Date(DO.obs.1$V1),cf*DO.obs.1$V2,col=col.obs,pch=20,cex=1)
lines(x,oxy[,bottom],lty=1,col=col.sim.bottom)
lines(x,oxy[,surface],lty=1,col=col.sim,lwd=2)
legend('bottomleft',
cex=1,
legend=c('E9-obs','E8-obs','AL-obs','surface-sim', 'bottom-sim'),
bty = 'n',
title("dissolved oxygen concentration mmol-O2 L-1", adj = 0, line = 0.1, cex = 0.5, font.main = 1)
points(as.Date(DO_obs_2$V1), (1000 / 32) * DO_obs_2$V2, col = col_obs2, pch = 0, cex = .5)
points(as.Date(DO_obs_1$V1), cf * DO_obs_1$V2, col = col_obs, pch = 20, cex = 1)
lines(x, oxy[,bottom], lty = 1, col = col_sim_bottom)
lines(x, oxy[,surface], lty = 1, col = col_sim, lwd = 2)
legend("bottomleft",
cex = 1,
legend = c("E9-obs", "E8-obs", "AL-obs", "surface-sim", "bottom-sim"),
bty = "n",
horiz = T,
border=NA,
pt.cex=c(.5,1,NA,NA,NA),
col=c(col.obs2,col.obs,col.obs3,col.sim,col.sim.bottom),
pch=c(0,20,0,0,0),
pt.lwd=c(1,1,1,1,1),
lwd=c(0,0,1,2,1)
border = NA,
pt.cex = c(.5, 1, NA, NA, NA),
col=c(col_obs2, col_obs, col_obs3, col_sim, col_sim_bottom),
pch = c(0, 20, 0, 0,0),
pt.lwd=c(1, 1, 1, 1, 1),
lwd=c(0, 0, 1, 2, 1)
)
dev.off()
\ No newline at end of file
testcases/estuary/setup_data.R
+ 20
− 20
View file @ c9f53b38
......@@ -2,21 +2,21 @@
# SPDX-License-Identifier: CC0-1.0
# SPDX-FileContributor Ovidio Garcia-Oliva <ovidio.garcia@hereon.de>
setwd('./data/')
setwd("./data/")
################################################################################
## definition of local functions
Get.Year = function(...) format(as.Date(..., format="%d/%m/%Y"),"%Y")
as.POSIXct = function(...) base::as.POSIXct(...,,format="%Y-%m-%d %H:%M:%S")
get.year = function(...) format(as.Date(..., format = "%d/%m/%Y"), "%Y")
as.POSIXct = function(...) base::as.POSIXct(..., format = "%Y-%m-%d %H:%M:%S")
################################################################################
## loading observed values for temperature (temp) wind velocity (wind) and wind
## direction (dirw) to create forcing file
temp = read.delim("Cuxhaven_DWD!Lufttemperatur.txt", header=FALSE, comment.char="#")
wind = read.delim("Cuxhaven_DWD!Windgeschwindigkeit.txt", header=FALSE, comment.char="#")
dirw = read.delim("Cuxhaven_DWD!Windrichtung.txt", header=FALSE, comment.char="#")
temp = read.delim("Cuxhaven_DWD!Lufttemperatur.txt", header = FALSE, comment.char = "#")
wind = read.delim("Cuxhaven_DWD!Windgeschwindigkeit.txt", header = FALSE, comment.char = "#")
dirw = read.delim("Cuxhaven_DWD!Windrichtung.txt", header = FALSE, comment.char = "#")
temp$V1 = as.POSIXct(temp$V1)
wind$V1 = as.POSIXct(wind$V1)
......@@ -24,34 +24,34 @@ dirw$V1 = as.POSIXct(dirw$V1)
## transforming temperatures at 9m above ground to 2 meter values
## assuming linear profile
h.station = 9
h.model = 2
T.gradient = -0.0065
h_station = 9
h_model = 2
T_gradient = -0.0065
bias = 5
temp$temp = temp$V2+(h.model-h.station)*T.gradient + bias
temp$temp = temp$V2 + (h_model - h_station) * T_gradient + bias
## transforming wind speed at 9m above ground to 10 meter values
## assuming logarithmic profile
h.station = 9
h.model = 10
h_station = 9
h_model = 10
w.exponent = 0.14
wind$V3 = wind$V2*(h.model/h.station)**w.exponent
wind$V3 = wind$V2 * (h_model / h_station)^w.exponent
## merging wind files
wind = merge(wind,dirw,by='V1',suffixes = c('.vel','.dir'))
wind = merge(wind, dirw, by = "V1", suffixes = c(".vel", ".dir"))
## calculating wind components
wind$u10 = wind$V3*cos(pi*wind$V2.dir/180)
wind$v10 = wind$V3*sin(pi*wind$V2.dir/180)
wind$u10 = wind$V3 * cos(pi*wind$V2.dir / 180)
wind$v10 = wind$V3 * sin(pi*wind$V2.dir / 180)
## creating the meteofile
meteofile = merge(temp,wind)
meteofile = merge(temp, wind)
meteofile = subset(meteofile,
Get.Year(V1)>=2004 & !is.na(temp) & !is.na(v10) & !is.na(u10),
select=c(V1,temp,u10,v10)
get.year(V1) >= 2004 & !is.na(temp) & !is.na(v10) & !is.na(u10),
select=c(V1, temp, u10, v10)
)
meteofile$V1 = meteofile$V1 + 1 # adding 1 sec to avoid a bug with gotm config
colnames(meteofile) = paste0("#", colnames(meteofile))
write.csv(meteofile,'meteofile.csv',row.names=F,quote=F)
write.csv(meteofile, "meteofile.csv", row.names=FALSE ,quote=FALSE)
testcases/light/model_comparison.R
+ 63
− 54
View file @ c9f53b38
......@@ -4,79 +4,88 @@
if (!require("ncdf4")) install.packages("ncdf4")
################################################################################
## definition of local functions
Get.Year = function(...) format(as.Date(..., format="%d/%m/%Y"),"%Y")
################################################################################
## loading the results of gotm
system('ln -f fabm.new.yaml fabm.yaml')
system('./gotm')
nc_data = nc_open('output.nc')
bpar = t(ncvar_get(nc_data, "light_par"))
bswr = t(ncvar_get(nc_data, "light_swr"))
system("ln -f fabm.new.yaml fabm.yaml")
system("./gotm")
nc_data = nc_open("output.nc")
bpar = t(ncvar_get(nc_data, "light_par"))
bswr = t(ncvar_get(nc_data, "light_swr"))
bphy = t(ncvar_get(nc_data, "oxypom_ALG1"))
system('ln -f fabm.ref.yaml fabm.yaml')
system('./gotm')
nc_data = nc_open('output.nc')
rpar = t(ncvar_get(nc_data, "light_par"))
rswr = t(ncvar_get(nc_data, "light_swr"))
system("ln -f fabm.ref.yaml fabm.yaml")
system("./gotm")
nc_data = nc_open("output.nc")
rpar = t(ncvar_get(nc_data, "light_par"))
rswr = t(ncvar_get(nc_data, "light_swr"))
rphy = t(ncvar_get(nc_data, "oxypom_ALG1"))
################################################################################
## plotting the data
col.sim='tomato'
col.ref='black'
col_sim = "tomato"
col_ref = "black"
N = length(rpar[1,])
png(filename="./light_validation.png",
png(filename = "./light_validation.png",
width = 600, height = 600)
par(mfrow=c(3,2),mai=2*c(0.42,0.42,0.21,0.21),oma=2*c(1,1,0.5,0.5),las=1)
hist(100*(bpar-rpar)/(rpar+bpar),
main='difference in par relatiave to reference (%)',
xlim=c(-100,100),
xlab='',
freq = F,
border=NA,
n=50)
hist(100*(bphy-rphy)/(rphy+bphy),
xlim=c(-100,100),
xlab='',
main='difference in ALG1 relatiave to reference (%)',
freq=F,
border=NA,
n=50)
plot(rpar,bpar,col='lightgray',pch=20,
xlab='reference ALG1',
ylab='dobgc_light ALG1')
abline(a=0,b=1)
par(mfrow = c(3, 2), mai = 2 * c(0.42, 0.42, 0.21, 0.21), oma = 2 * c(1, 1, 0.5, 0.5), las = 1)
hist(100 * (bpar - rpar) / (rpar + bpar),
main = "difference in par relatiave to reference (%)",
xlim = c(-100, 100),
xlab = "",
freq = FALSE,
border = NA,
n = 50
)
hist(100 * (bphy - rphy) / (rphy + bphy),
xlim = c(-100, 100),
xlab = "",
main = "difference in ALG1 relatiave to reference (%)",
freq = FALSE,
border = NA,
n = 50
)
plot(rpar, bpar,
col = "lightgray",
pch = 20,
xlab = "reference ALG1",
ylab = "dobgc_light ALG1"
)
abline(a = 0, b = 1)
plot(rphy,bphy,col='lightgray',pch=20,
xlab='reference ALG1',
ylab='dobgc_light ALG1')
abline(a=0,b=1)
plot(rowMeans(rpar[,1:(N/2+1)]),type='l',col=col.ref,
main='par in the surface',
xlab='',
ylab='W m-2')
lines(rowMeans(bpar[,1:(N/2+1)]),type='l',col=col.sim)
plot(rphy, bphy,
col = "lightgray",
pch = 20,
xlab = "reference ALG1",
ylab = "dobgc_light ALG1"
)
abline(a = 0, b = 1)
plot(rowMeans(rpar[,1:(N / 2 + 1)]),
type = "l",
col = col_ref,
main = "par in the surface",
xlab = "",
ylab = "W m-2"
)
lines(rowMeans(bpar[,1:(N / 2 + 1)]), type = "l", col = col_sim)
plot(rowMeans(rphy[,1:(N/2+1)]),type='l',col=col.ref,
main='ALG1 in the surface',
xlab='',
ylab='mmol-C m-3')
lines(rowMeans(bphy[,1:(N/2+1)]),type='l',col=col.sim)
plot(rowMeans(rphy[,1:(N / 2 + 1)]),
type = "l",
col = col_ref,
main = "ALG1 in the surface",
xlab = "",
ylab = "mmol-C m-3"
)
lines(rowMeans(bphy[,1:(N / 2 + 1)]),type = "l",col = col_sim)
dev.off()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

A HIFIS Service | Privacy | Imprint | Support | Documentation | Changelog