forked from eco4cast/neon4cast-ci
-
Notifications
You must be signed in to change notification settings - Fork 2
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #23 from eco4cast/stage2-patch
Stage2 patch
- Loading branch information
Showing
4 changed files
with
167 additions
and
18 deletions.
There are no files selected for viewing
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
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
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,136 @@ | ||
|
|
||
| #' | ||
| #' @param df dataframe of stage1 NEON GEFS forecasts for sites to forecast at | ||
| #' @param site_list a dataframe of the latitude and longitude for all site_ids in df | ||
| #' @param use_solar_geom logical for using solar geometry for daily SW calculation | ||
| #' @param psuedo logical for something... | ||
| to_hourly <- function(df, | ||
| site_list, | ||
| use_solar_geom = TRUE, | ||
| psuedo = FALSE){ | ||
|
|
||
| if(!psuedo){ | ||
| reference_datetime <- lubridate::as_datetime(df$reference_datetime)[1] | ||
| }else{ | ||
| reference_datetime <- NA | ||
| } | ||
|
|
||
| var_order <- names(df) | ||
|
|
||
| ensemble_maxtime <- df |> | ||
| dplyr::group_by(site_id, family, ensemble) |> | ||
| dplyr::summarise(max_time = max(datetime), .groups = "drop") | ||
|
|
||
| ensembles <- unique(df$ensemble) | ||
| datetime <- seq(min(df$datetime), max(df$datetime), by = "1 hour") | ||
| variables <- unique(df$variable) | ||
| sites <- unique(df$site_id) | ||
|
|
||
| full_time <- expand.grid(sites, ensembles, datetime, variables) |> | ||
| dplyr::rename(site_id = Var1, | ||
| ensemble = Var2, | ||
| datetime = Var3, | ||
| variable = Var4) |> | ||
| dplyr::mutate(datetime = lubridate::as_datetime(datetime)) |> | ||
| dplyr::arrange(site_id, ensemble, variable, datetime) |> | ||
| dplyr::left_join(ensemble_maxtime, by = c("site_id","ensemble")) |> | ||
| dplyr::filter(datetime <= max_time) |> | ||
| dplyr::select(-c("max_time")) |> | ||
| dplyr::distinct() | ||
|
|
||
| states <- df |> | ||
| dplyr::select(site_id, family, horizon, ensemble, datetime, variable, prediction) |> | ||
| dplyr::filter(!psuedo | (psuedo & horizon != "006") | (psuedo & datetime == max(df$datetime))) |> | ||
| dplyr::select(-horizon) |> | ||
| dplyr::group_by(site_id, family, ensemble, variable) |> | ||
| dplyr::right_join(full_time, by = c("site_id", "ensemble", "datetime", "family", "variable")) |> | ||
| dplyr::filter(variable %in% c("PRES", "RH", "TMP", "UGRD", "VGRD")) |> | ||
| dplyr::arrange(site_id, family, ensemble, datetime) |> | ||
| dplyr::mutate(prediction = imputeTS::na_interpolation(prediction, option = "linear")) |> | ||
| dplyr::mutate(prediction = ifelse(variable == "TMP", prediction + 273, prediction)) |> | ||
| dplyr::mutate(prediction = ifelse(variable == "RH", prediction/100, prediction)) |> | ||
| dplyr::ungroup() | ||
|
|
||
| fluxes <- df |> | ||
| dplyr::select(site_id, family, horizon, ensemble, datetime, variable, prediction) |> | ||
| dplyr::filter(horizon != "003") |> | ||
| dplyr::select(-horizon) |> | ||
| dplyr::group_by(site_id, family, ensemble, variable) |> | ||
| dplyr::right_join(full_time, by = c("site_id", "ensemble", "datetime", "family", "variable")) |> | ||
| dplyr::filter(variable %in% c("APCP","DSWRF","DLWRF")) |> | ||
| dplyr::arrange(site_id, family, ensemble, datetime) |> | ||
| tidyr::fill(prediction, .direction = "up") |> | ||
| dplyr::mutate(prediction = ifelse(variable == "APCP", prediction / (6 * 60 * 60), prediction), | ||
| variable = ifelse(variable == "APCP", "PRATE", variable)) |> | ||
| dplyr::ungroup() | ||
|
|
||
| if(use_solar_geom){ | ||
|
|
||
| fluxes <- fluxes |> | ||
| dplyr::left_join(site_list, by = "site_id") |> | ||
| dplyr::mutate(hour = lubridate::hour(datetime), | ||
| date = lubridate::as_date(datetime), | ||
| doy = lubridate::yday(datetime) + hour/24, | ||
| longitude = ifelse(longitude < 0, 360 + longitude, longitude), | ||
| rpot = downscale_solar_geom(doy, longitude, latitude)) |> # hourly sw flux calculated using solar geometry | ||
| dplyr::group_by(site_id, family, ensemble, date, variable) |> | ||
| dplyr::mutate(avg.rpot = mean(rpot, na.rm = TRUE), | ||
| avg.SW = mean(prediction, na.rm = TRUE))|> # daily sw mean from solar geometry | ||
| dplyr::ungroup() |> | ||
| dplyr::mutate(prediction = ifelse(variable == "DSWRF" & avg.rpot > 0.0, rpot * (avg.SW/avg.rpot),prediction)) |> | ||
| dplyr::select(any_of(var_order)) | ||
| } | ||
|
|
||
| hourly_df <- dplyr::bind_rows(states, fluxes) |> | ||
| dplyr::arrange(site_id, family, ensemble, datetime) |> | ||
| dplyr::mutate(variable = ifelse(variable == "TMP", "air_temperature", variable), | ||
| variable = ifelse(variable == "PRES", "air_pressure", variable), | ||
| variable = ifelse(variable == "RH", "relative_humidity", variable), | ||
| variable = ifelse(variable == "DLWRF", "surface_downwelling_longwave_flux_in_air", variable), | ||
| variable = ifelse(variable == "DSWRF", "surface_downwelling_shortwave_flux_in_air", variable), | ||
| variable = ifelse(variable == "PRATE", "precipitation_flux", variable), | ||
| variable = ifelse(variable == "VGRD", "eastward_wind", variable), | ||
| variable = ifelse(variable == "UGRD", "northward_wind", variable), | ||
| variable = ifelse(variable == "APCP", "precipitation_amount", variable), | ||
| reference_datetime = reference_datetime) |> | ||
| dplyr::select(any_of(var_order)) | ||
|
|
||
| return(hourly_df) | ||
|
|
||
| } | ||
|
|
||
| cos_solar_zenith_angle <- function(doy, lat, lon, dt, hr) { | ||
| et <- equation_of_time(doy) | ||
| merid <- floor(lon / 15) * 15 | ||
| merid[merid < 0] <- merid[merid < 0] + 15 | ||
| lc <- (lon - merid) * -4/60 ## longitude correction | ||
| tz <- merid / 360 * 24 ## time zone | ||
| midbin <- 0.5 * dt / 86400 * 24 ## shift calc to middle of bin | ||
| t0 <- 12 + lc - et - tz - midbin ## solar time | ||
| h <- pi/12 * (hr - t0) ## solar hour | ||
| dec <- -23.45 * pi / 180 * cos(2 * pi * (doy + 10) / 365) ## declination | ||
| cosz <- sin(lat * pi / 180) * sin(dec) + cos(lat * pi / 180) * cos(dec) * cos(h) | ||
| cosz[cosz < 0] <- 0 | ||
| return(cosz) | ||
| } | ||
|
|
||
| equation_of_time <- function(doy) { | ||
| stopifnot(doy <= 367) | ||
| 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 | ||
| return(et) | ||
| } | ||
|
|
||
| downscale_solar_geom <- function(doy, lon, lat) { | ||
|
|
||
| dt <- median(diff(doy)) * 86400 # average number of seconds in time interval | ||
| hr <- (doy - floor(doy)) * 24 # hour of day for each element of doy | ||
|
|
||
| ## calculate potential radiation | ||
| cosz <- cos_solar_zenith_angle(doy, lat, lon, dt, hr) | ||
| rpot <- 1366 * cosz | ||
| return(rpot) | ||
| } | ||
|
|
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