From c9f53b38e1b9440b550ad2b467b71c82aba2b58a Mon Sep 17 00:00:00 2001 From: og <ovidio.garcia> Date: Fri, 21 Mar 2025 07:11:52 +0100 Subject: [PATCH] lintr formatting --- testcases/estuary/plot_output.R | 136 ++++++++++++++++------------- testcases/estuary/setup_data.R | 40 ++++----- testcases/light/model_comparison.R | 117 +++++++++++++------------ 3 files changed, 156 insertions(+), 137 deletions(-) diff --git a/testcases/estuary/plot_output.R b/testcases/estuary/plot_output.R index 91f6507..3202bdf 100755 --- a/testcases/estuary/plot_output.R +++ b/testcases/estuary/plot_output.R @@ -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 diff --git a/testcases/estuary/setup_data.R b/testcases/estuary/setup_data.R index ac23119..2ff02cb 100644 --- a/testcases/estuary/setup_data.R +++ b/testcases/estuary/setup_data.R @@ -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) diff --git a/testcases/light/model_comparison.R b/testcases/light/model_comparison.R index 3fb485b..e7620a3 100755 --- a/testcases/light/model_comparison.R +++ b/testcases/light/model_comparison.R @@ -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() -- GitLab