Tidy up leap year implementation, and other date-related cleanup #1623

ashiklom · 2017-09-01T14:28:58Z

Description

New function PEcAn.utils::days_in_year(year) calculates the days in a year given a year (or vector of years). I used this function to replace many occurrences of ifelse(leap_year(year), 366, 355) many places in the code. I also simplified many statements related to leap year if-else clauses, where it was simple enough to do so.

Also, new function PEcAn.data.atmsphere::eccentricity_obliquity(doy) does the appropriate math. This was refactored out of several places that had this fairly complex equation duplicated.

Also, all of the files I touched have had trailing whitespace removed (courtesy of RStudio), which extends the review diff, but is good coding practice.

Motivation and Context

See issue #801.

Types of changes

Bug fix (non-breaking change which fixes an issue)
New feature (non-breaking change which adds functionality)
Breaking change (fix or feature that would cause existing functionality to change)

Checklist:

My change requires a change to the documentation.
I have updated the CHANGELOG.md.
I have updated the documentation accordingly.
I have read the CONTRIBUTING document.
I have added tests to cover my changes.
All new and existing tests passed.

New function `PEcAn.utils::days_in_year(year)` calculates the days in a year given a year (or vector of years). I used this function to replace many occurrences of `ifelse(leap_year(year), 366, 355)` many places in the code. I also simplified many statements related to leap year if-else clauses, where it was simple enough to do so. Also, new function `PEcAn.data.atmsphere::eccentricity_obliquity(doy)` does the appropriate math. This was refactored out of several places that had this fairly complex equation duplicated. Also, all of the files I touched have had trailing whitespace removed (courtesy of RStudio), which extends the review diff, but is good coding practice.

For completeness, even though it doesn't actually do anything yet.

mdietze · 2017-09-01T14:37:35Z

models/ed/R/met2model.ED2.R

-    f                <- pi/180 * (279.5 + 0.9856 * doy)
-    et               <- (-104.7 * sin(f) + 596.2 * sin(2 * f) + 4.3 * sin(4 * f) - 429.3 * 
-                           cos(f) - 2 * cos(2 * f) + 19.3 * cos(3 * f)) / 3600  # equation of time -> eccentricity and obliquity
+    et <- PEcAn.data.atmosphere::eccentricity_obliquity(doy)


I think it's great to encapsulate this, but you could probably go even further and encapsulate the whole solar geometry calculation of solar zenith angle as well (i.e. everything up to cosz below)

Yeah, I was thinking that too. I went through it quickly and pulled this out because it was really obvious and I saw it in a few places. But yeah, I'll go ahead and do that, so wait to merge this until I do.

… date-cleanup

mdietze · 2017-09-01T16:31:13Z

@ashiklom you've got a ROxygen fail on models/jules/man/detect.timestep.Rd

dlebauer · 2017-09-01T17:14:20Z

models/dalec/R/met2model.DALEC.R

-           dt <- (366 * 24 * 60 * 60) / length(sec), # leap year 
-           dt <- (365 * 24 * 60 * 60) / length(sec)) # non-leap year
+    diy <- PEcAn.utils::days_in_year(year)
+    dt <- diy * 24 * 60 * 60 / length(sec)


wouldn't the convention be to use

udunits::ud.convert(diy, 'days', 'seconds') / length(sec)

dlebauer · 2017-09-01T17:16:25Z

models/ed/R/met2model.ED2.R

-        dtmp <- rep(1:366, each = 86400 / dt)
-      }
+      diy <- PEcAn.utils::days_in_year(y)
+      ytmp <- rep(y, diy * 86400 / dt)


another case where seconds / day is hardcoded. Better to use ud.convert?

dlebauer · 2017-09-01T17:17:19Z

models/ed/R/met2model.ED2.R

@@ -146,7 +145,7 @@ met2model.ED2 <- function(in.path, in.prefix, outfolder, start_date, end_date, l
    if (useCO2) {
      CO2 <- c(rep(CO2[1], toff), CO2)[1:slen]
    }
-    
+
    ## build time variables (year, month, day of year)
    skip <- FALSE
    nyr  <- floor(length(sec) / 86400 / 365 * dt)


looks like a kludge. Could ud.convert work here too?

dlebauer · 2017-09-01T17:18:07Z

models/ed/R/met2model.ED2.R

-
+
+    diy <- PEcAn.utils::days_in_year(year)
+    dt <- diy * 24 * 60 * 60 / length(sec)


use ud.convert? This pattern seems to be repeated a lot. Maybe we need another function seconds_in_year?

dlebauer · 2017-09-01T17:19:39Z

models/ed/R/met2model.ED2.R

-    cosz             <- sin(lat * pi/180) * sin(dec) + cos(lat * pi/180) * cos(dec) * cos(h)
-    cosz[cosz < 0]   <- 0
-
+    cosz <- PEcAn.data.atmosphere::solar_angle(doy, lat, lon, dt)


does this function return solar_angle, or cos(zenith angle)? I would expect that solar_angle to return zenith_angle, and perhaps function should be named cos_solar_zenith_angle???

dlebauer · 2017-09-01T17:22:40Z

models/jules/R/write.config.JULES.R

@@ -628,7 +628,8 @@ detect.timestep <- function(met.dir,met.regexp,start_date){
      tlen <- grep("time =", met.header)
      if (length(tlen) > 0) {
        tlen <- as.numeric(gsub(pattern = "[^[:digit:]]", "", met.header[tlen]))
-        dt <- 86400 / round(tlen/(365 + lubridate::leap_year(as.Date(start_date))))
+        diy <- PEcAn.utils::days_in_year(lubridate::year(as.Date(start_date)))
+        dt <- 86400 / round(tlen/(diy))


another repeated could this be wrapped in a function called dt?

dlebauer · 2017-09-01T17:23:59Z

models/preles/R/runPRELES.jobsh.R

      VPD <- VPD * 0.01  # convert to Pascal
-      
+


these transformations should also use ud.convert

dlebauer · 2017-09-01T17:24:35Z

models/sipnet/R/met2model.SIPNET.R

-                   366 * 24 * 60 * 60 / length(sec), # leap year
-                   365 * 24 * 60 * 60 / length(sec)) # non-leap year
+
+      dt <- diy * 24 * 60 * 60 / length(sec)


dlebauer · 2017-09-01T17:25:39Z

modules/assim.batch/R/get.da.data.R

@@ -43,33 +43,33 @@ calculate.nee.L <- function(yeardoytime, model.i.nee, observed.flux, be, bu) {
 get.da.data <- function(out.dir, ameriflux.dir, years, be, bu, ensemble.size = 199) {
  ensemble.size <- 500


is this hard-coded value of ensemble size a bug?

dlebauer · 2017-09-01T17:29:02Z

modules/data.atmosphere/R/download.MsTMIP_NARR.R

-    
-    ntime <- ifelse(year %% 4 == 0, 2923, 2919) # should use lubridate::leap_year?
+
+    ntime <- ifelse(lubridate::leap_year(year), 2923, 2919)   ## ANS: Where do these numbers come from?


these numbers are the number of 3 hour intervals in a year -1 (365*8=2920)

perhaps can be replaced by days_in_year(year) * 8 - 1?

dlebauer · 2017-09-01T17:31:36Z

modules/data.atmosphere/R/solar_angle.R

+#' @param doy Day of year
+#' @param lat Latitude
+#' @param lon Longitude
+#' @param dt Timestep


could be better documented.

What is the return value?

cite a canonical source of the equations used

add description

Should the name be cos_solar_zenith_angle to be more explicit?

dlebauer · 2017-09-01T17:31:56Z

modules/data.atmosphere/R/solar_angle.R

+#' @author Alexey Shiklomanov
+#' @param doy Day of year
+#' @export
+eccentricity_obliquity <- function(doy) {


dlebauer · 2017-09-01T17:33:14Z

modules/data.atmosphere/inst/scripts/extract2driver.R

-			eccentricity * 
-			(cos(declination) * 
-				cos(deg2rad(LAT)) * 
+	radiation = 1367.0 *


magic number !!!also uses = instead of <-

This is an old script (not a package function definition) that hasn't been touched in a while, so I'm just gonna leave it. There's a lot more that's wrong with it, which I don't feel like spending time fixing.

dlebauer

Looks like a good start

I made a few comments - it seems that there is an opportunity to clean up some more magic numbers as seconds/ day is repeatedly calculated or given as a magic number. dt <- diy * 24 * 60 * 60 / length(sec)

I found a few bits of code that you didn't explicitly touch but are worth looking at.

Also not sure what the current thinking on roll-your-own vs. use existing packages, but there is already the insol package that calculates solar zenith angle and potential radiation

@dlebauer

Per @dlebauer recommendations.

Use it in all the `met2model` code instead of hard-coded math.

Also, Rstudio removed trailing whitespace from a bunch of files.

Also, some more whitespace fixes by RStudio.

ashiklom · 2017-09-05T15:21:11Z

@dlebauer I addressed many of your comments, but to address all of the cases where we have had-coded date-time conversions would take a while. I definitely agree that all of our met processing and met2model code could use some serious cleanup and refactoring, but I don't have the time to sift through every line to figure out what the conversions are (especially since I wrote none of this original code, so it's often hard for me to divine exactly what's going on). I think this PR makes decent progress on a lot of the required cleanup, and it's worth merging now and continuing the process in later pull requests. If you or @mdietze disagree, you can request additional changes to the PR, but it'll be a little while before I get to them since I have a lot of other stuff to catch up on.

dlebauer · 2017-09-05T15:47:35Z

models/preles/R/runPRELES.jobsh.R


      ## Format/convert inputs
      ppfd   <- tapply(PPFD, doy, mean, na.rm = TRUE)  # Find the mean for the day
      tair   <- udunits2::ud.convert(tapply(Tair, doy, mean, na.rm = TRUE), "kelvin", "celsius")  # Convert Kelvin to Celcius
      vpd    <- udunits2::ud.convert(tapply(VPD, doy, mean, na.rm = TRUE), "Pa", "kPa")  # pascal to kila pascal
      precip <- tapply(Precip, doy, sum, na.rm = TRUE)  # Sum to daily precipitation
      co2    <- tapply(CO2, doy, mean)  # need daily average, so sum up day
-      co2    <- co2 / 1e+06  # convert to ppm
+      co2    <- co2 / 1e+06  # convert to ppm. ANS: Convert from what? Mole-fraction to ppm is multiplying by 1e6, not dividing


Is this an error that is currently causing invalid calculations? If so, who can resolve it?

@tonygardella wrote this coupler, so I think he can check and resolve this.

I'll fix it!

ashiklom added 3 commits September 1, 2017 07:32

Met: Replace year %% 4 with leap_year

48a40b3

Update Changelog.

9983509

ashiklom requested a review from dlebauer September 1, 2017 14:30

ashiklom self-assigned this Sep 1, 2017

ashiklom added the Type: Enhancement label Sep 1, 2017

Add leap year fix to met2model.CLM45.

6516e96

For completeness, even though it doesn't actually do anything yet.

mdietze approved these changes Sep 1, 2017

View reviewed changes

mdietze and others added 4 commits September 1, 2017 10:51

Merge branch 'develop' into date-cleanup

8822bd1

Met: Add PEcAn.data.atmosphere::solar_angle function

a150568

Merge branch 'date-cleanup' of https://github.com/ashiklom/pecan into…

a58269b

… date-cleanup

Update changelog

9625c84

dlebauer reviewed Sep 1, 2017

View reviewed changes

models/preles/R/runPRELES.jobsh.R

VPD <- VPD * 0.01 # convert to Pascal

Copy link

Member

dlebauer Sep 1, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

these transformations should also use ud.convert

dlebauer reviewed Sep 1, 2017

View reviewed changes

Fix documentation

b767d88

dlebauer reviewed Sep 1, 2017

View reviewed changes

dlebauer requested changes Sep 1, 2017

View reviewed changes

ashiklom added 2 commits September 5, 2017 09:41

Improve documentation for solar angle math

a04a3d1

Per @dlebauer recommendations.

Add utils::seconds_in_year function

b5ce058

Use it in all the `met2model` code instead of hard-coded math.

ashiklom added 3 commits September 5, 2017 10:21

Replace +/- 273.15 with udunits T conversions

cfaeda0

Also, Rstudio removed trailing whitespace from a bunch of files.

Change manual unit conversions to udunits2

7f32455

Also, some more whitespace fixes by RStudio.

assim.batch: Fix ignored ensemble size argument

fa071c7

ashiklom added 2 commits September 5, 2017 11:24

ED: Cleanup met2model.ED2 library calls

bd1ce91

data.atmosphere: More cleanup of met NAMESPACE

f6ff4a2

dlebauer approved these changes Sep 5, 2017

View reviewed changes

dlebauer reviewed Sep 5, 2017

View reviewed changes

dlebauer approved these changes Sep 5, 2017

View reviewed changes

ashiklom added 2 commits September 5, 2017 12:06

Merge with upstream

d719dd2

Update changelog

749a72e

mdietze merged commit 009e165 into PecanProject:develop Sep 5, 2017

ashiklom mentioned this pull request Sep 7, 2017

incorrect identification of leap years; multiple implementations. #801

Closed

3 tasks

ashiklom deleted the date-cleanup branch October 3, 2017 21:34

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Tidy up leap year implementation, and other date-related cleanup #1623

Tidy up leap year implementation, and other date-related cleanup #1623

ashiklom commented Sep 1, 2017 •

edited

mdietze Sep 1, 2017

ashiklom Sep 1, 2017

mdietze commented Sep 1, 2017

dlebauer Sep 1, 2017

dlebauer Sep 1, 2017

dlebauer Sep 1, 2017

dlebauer Sep 1, 2017

dlebauer Sep 1, 2017

dlebauer Sep 1, 2017

dlebauer Sep 1, 2017

dlebauer Sep 1, 2017

dlebauer Sep 1, 2017

dlebauer Sep 1, 2017

dlebauer Sep 1, 2017

dlebauer Sep 1, 2017

dlebauer Sep 1, 2017

dlebauer Sep 1, 2017

ashiklom Sep 5, 2017

dlebauer left a comment

ashiklom commented Sep 5, 2017

dlebauer Sep 5, 2017

mdietze Sep 5, 2017

tonygardella Sep 5, 2017



		diy <- PEcAn.utils::days_in_year(year)
		dt <- diy * 24 * 60 * 60 / length(sec)

		@@ -43,33 +43,33 @@ calculate.nee.L <- function(yeardoytime, model.i.nee, observed.flux, be, bu) {
		get.da.data <- function(out.dir, ameriflux.dir, years, be, bu, ensemble.size = 199) {
		ensemble.size <- 500


		ntime <- ifelse(year %% 4 == 0, 2923, 2919) # should use lubridate::leap_year?

		ntime <- ifelse(lubridate::leap_year(year), 2923, 2919) ## ANS: Where do these numbers come from?

Tidy up leap year implementation, and other date-related cleanup #1623

Tidy up leap year implementation, and other date-related cleanup #1623

Conversation

ashiklom commented Sep 1, 2017 • edited

Description

Motivation and Context

Types of changes

Checklist:

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mdietze commented Sep 1, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dlebauer left a comment

Choose a reason for hiding this comment

ashiklom commented Sep 5, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ashiklom commented Sep 1, 2017 •

edited