commit 046c8139e59ad975704c9a01852bb7dd7b42f725 Author: boB Rudis Date: Thu Jan 25 17:40:26 2018 -0500 initial commit diff --git a/.Rbuildignore b/.Rbuildignore new file mode 100644 index 0000000..42741f1 --- /dev/null +++ b/.Rbuildignore @@ -0,0 +1,11 @@ +^.*\.Rproj$ +^\.Rproj\.user$ +^\.travis\.yml$ +^README\.*Rmd$ +^README\.*html$ +^NOTES\.*Rmd$ +^NOTES\.*html$ +^\.codecov\.yml$ +^README_files$ +^doc$ +^terminator-animation\.gif$ \ No newline at end of file diff --git a/.codecov.yml b/.codecov.yml new file mode 100644 index 0000000..69cb760 --- /dev/null +++ b/.codecov.yml @@ -0,0 +1 @@ +comment: false diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..cce1f17 --- /dev/null +++ b/.gitignore @@ -0,0 +1,8 @@ +.DS_Store +.Rproj.user +.Rhistory +.RData +.Rproj +src/*.o +src/*.so +src/*.dll diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 0000000..76d9586 --- /dev/null +++ b/.travis.yml @@ -0,0 +1,31 @@ +language: r + +warnings_are_errors: true + +sudo: required + +cache: packages + +r: + - oldrel + - release + - devel + +apt_packages: + - libv8-dev + - xclip + +env: + global: + - CRAN: http://cran.rstudio.com + +after_success: + - Rscript -e 'covr::codecov()' + +notifications: + email: + - bob@rud.is + irc: + channels: + - "104.236.112.222#builds" + nick: travisci diff --git a/DESCRIPTION b/DESCRIPTION new file mode 100644 index 0000000..8170b55 --- /dev/null +++ b/DESCRIPTION @@ -0,0 +1,27 @@ +Package: terminator +Type: Package +Title: Compute Global Terminator (Day/Night) Bands +Version: 0.1.0 +Date: 2018-01-25 +Authors@R: c( + person("Bob", "Rudis", email = "bob@rud.is", role = c("aut", "cre"), + comment = c(ORCID = "0000-0001-5670-2640")), + person("Joe", "Gallagher", email = "joedgallagher@gmail.com", role = "aut") + ) +Encoding: UTF-8 +Maintainer: Bob Rudis +Description: A good description goes here otherwise CRAN checks fail. +URL: https://github.com/hrbrmstr/terminator +BugReports: https://github.com/hrbrmstr/terminator/issues +License: AGPL +Suggests: + testthat, + covr +Depends: + R (>= 3.2.0) +Imports: + purrr, + Rcpp, + utils +RoxygenNote: 6.0.1.9000 +LinkingTo: Rcpp diff --git a/NAMESPACE b/NAMESPACE new file mode 100644 index 0000000..46c6514 --- /dev/null +++ b/NAMESPACE @@ -0,0 +1,8 @@ +# Generated by roxygen2: do not edit by hand + +export(terminator) +export(terminator_lat_lon) +importFrom(Rcpp,sourceCpp) +importFrom(purrr,map_df) +importFrom(utils,globalVariables) +useDynLib(terminator) diff --git a/NEWS.md b/NEWS.md new file mode 100644 index 0000000..9b4679b --- /dev/null +++ b/NEWS.md @@ -0,0 +1,2 @@ +0.1.0 +* Initial release diff --git a/R/RcppExports.R b/R/RcppExports.R new file mode 100644 index 0000000..479163e --- /dev/null +++ b/R/RcppExports.R @@ -0,0 +1,19 @@ +# Generated by using Rcpp::compileAttributes() -> do not edit by hand +# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +#' Compute a single termiantor band +#' +#' Returns a dataframe of latitude and longitude for the line that separates illuminated +#' day and dark night for any given time +#' +#' @md +#' @param time time (numeric from `POSIXct`) for the computation (bands are time-dependent) +#' @param from,to,by latitude sequence setup +#' @return data frame +#' @references , +#' +#' @export +terminator <- function(time, from = -180, to = 180, by = 0.1) { + .Call('_terminator_terminator', PACKAGE = 'terminator', time, from, to, by) +} + diff --git a/R/terminator-package.R b/R/terminator-package.R new file mode 100644 index 0000000..8358e0d --- /dev/null +++ b/R/terminator-package.R @@ -0,0 +1,10 @@ +#' Compute Global Terminator (Day/Night) Bands +#' +#' @name terminator +#' @docType package +#' @author Bob Rudis (bob@@rud.is) +#' @importFrom purrr map_df +#' @importFrom utils globalVariables +#' @useDynLib terminator +#' @importFrom Rcpp sourceCpp +NULL \ No newline at end of file diff --git a/R/terminator_lat_lon.R b/R/terminator_lat_lon.R new file mode 100644 index 0000000..7a4a15a --- /dev/null +++ b/R/terminator_lat_lon.R @@ -0,0 +1,15 @@ +utils::globalVariables(c(".x")) + +#' Generate a full set of terminator frames +#' +#' @param day,from,to,by setup for the hours sequence +#' @export +terminator_lat_lon <- function(day = Sys.Date(), from=0, to=23, by=1) { + + purrr::map_df(seq(0, 23, 1), ~{ + out <- terminator(as.integer((as.POSIXct(day) + (60*60*.x))), -180, 190, 0.5) + out$frame <- .x + out + }) + +} \ No newline at end of file diff --git a/README.Rmd b/README.Rmd new file mode 100644 index 0000000..4846a94 --- /dev/null +++ b/README.Rmd @@ -0,0 +1,57 @@ +--- +output: rmarkdown::github_document +editor_options: + chunk_output_type: console +--- + +# terminator + +Compute Global Terminator (Day/Night) Bands + +## Description + +Compute global terminator (day/night) bands which can be overlayed as day and night regions on a ggplot2 world map. + +## What's Inside The Tin + +The following functions are implemented: + +- `terminator`: Compute a single termiantor band +- `terminator_lat_lon`: Generate a full set of terminator frames + +## Installation + +```{r eval=FALSE} +devtools::install_github("hrbrmstr/terminator") +``` + +```{r message=FALSE, warning=FALSE, error=FALSE, include=FALSE} +options(width=120) +``` + +## Usage + +```{r message=FALSE, warning=FALSE, error=FALSE, eval=TRUE} +library(terminator) +library(ggplot2) +library(gganimate) # devtools::install_github("dgrtwo/gganimate") + +# current verison +packageVersion("terminator") + +term_seq <- terminator_lat_lon() + +chart <- ggplot(term_seq, aes(frame = frame)) + + borders("world", colour = "gray90", fill = "gray85") + + geom_ribbon(aes(lat, ymax = lon), ymin = 90, alpha = 0.2) + + coord_equal(xlim = c(-180, 190), ylim = c(-58, 85), expand = 0) + + ggthemes::theme_map() + +gganimate( + chart, + interval = 0.1, ani.width=1000, ani.height=600, + filename = "terminator-animation.gif" +) +``` + +![](terminator-animation.gif) diff --git a/README.md b/README.md new file mode 100644 index 0000000..73d8de9 --- /dev/null +++ b/README.md @@ -0,0 +1,53 @@ + +# terminator + +Compute Global Terminator (Day/Night) Bands + +## Description + +Compute global terminator (day/night) bands which can be overlayed as +day and night regions on a ggplot2 world map. + +## What’s Inside The Tin + +The following functions are implemented: + + - `terminator`: Compute a single termiantor band + - `terminator_lat_lon`: Generate a full set of terminator frames + +## Installation + +``` r +devtools::install_github("hrbrmstr/terminator") +``` + +## Usage + +``` r +library(terminator) +library(ggplot2) +library(gganimate) # devtools::install_github("dgrtwo/gganimate") + +# current verison +packageVersion("terminator") +``` + + ## [1] '0.1.0' + +``` r +term_seq <- terminator_lat_lon() + +chart <- ggplot(term_seq, aes(frame = frame)) + + borders("world", colour = "gray90", fill = "gray85") + + geom_ribbon(aes(lat, ymax = lon), ymin = 90, alpha = 0.2) + + coord_equal(xlim = c(-180, 190), ylim = c(-58, 85), expand = 0) + + ggthemes::theme_map() + +gganimate( + chart, + interval = 0.1, ani.width=1000, ani.height=600, + filename = "terminator-animation.gif" +) +``` + +![](terminator-animation.gif) diff --git a/man/terminator.Rd b/man/terminator.Rd new file mode 100644 index 0000000..04a1996 --- /dev/null +++ b/man/terminator.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R, R/terminator-package.R +\docType{package} +\name{terminator} +\alias{terminator} +\alias{terminator-package} +\title{Compute a single termiantor band} +\usage{ +terminator(time, from = -180, to = 180, by = 0.1) +} +\arguments{ +\item{time}{time (numeric from \code{POSIXct}) for the computation (bands are time-dependent)} + +\item{from, to, by}{latitude sequence setup} +} +\value{ +data frame +} +\description{ +Returns a dataframe of latitude and longitude for the line that separates illuminated +day and dark night for any given time +} +\references{ +\url{https://github.com/joergdietrich/Leaflet.Terminator/blob/master/L.Terminator.js}, +\url{https://github.com/JoGall/terminator/blob/master/terminator.R} +} +\author{ +Bob Rudis (bob@rud.is) +} diff --git a/man/terminator_lat_lon.Rd b/man/terminator_lat_lon.Rd new file mode 100644 index 0000000..52db851 --- /dev/null +++ b/man/terminator_lat_lon.Rd @@ -0,0 +1,14 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/terminator_lat_lon.R +\name{terminator_lat_lon} +\alias{terminator_lat_lon} +\title{Generate a full set of terminator frames} +\usage{ +terminator_lat_lon(day = Sys.Date(), from = 0, to = 23, by = 1) +} +\arguments{ +\item{day, from, to, by}{setup for the hours sequence} +} +\description{ +Generate a full set of terminator frames +} diff --git a/src/.gitignore b/src/.gitignore new file mode 100644 index 0000000..22034c4 --- /dev/null +++ b/src/.gitignore @@ -0,0 +1,3 @@ +*.o +*.so +*.dll diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp new file mode 100644 index 0000000..255cddf --- /dev/null +++ b/src/RcppExports.cpp @@ -0,0 +1,31 @@ +// Generated by using Rcpp::compileAttributes() -> do not edit by hand +// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +#include + +using namespace Rcpp; + +// terminator +DataFrame terminator(int time, double from, double to, double by); +RcppExport SEXP _terminator_terminator(SEXP timeSEXP, SEXP fromSEXP, SEXP toSEXP, SEXP bySEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< int >::type time(timeSEXP); + Rcpp::traits::input_parameter< double >::type from(fromSEXP); + Rcpp::traits::input_parameter< double >::type to(toSEXP); + Rcpp::traits::input_parameter< double >::type by(bySEXP); + rcpp_result_gen = Rcpp::wrap(terminator(time, from, to, by)); + return rcpp_result_gen; +END_RCPP +} + +static const R_CallMethodDef CallEntries[] = { + {"_terminator_terminator", (DL_FUNC) &_terminator_terminator, 4}, + {NULL, NULL, 0} +}; + +RcppExport void R_init_terminator(DllInfo *dll) { + R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); +} diff --git a/src/terminator-main.cpp b/src/terminator-main.cpp new file mode 100644 index 0000000..72dd2b0 --- /dev/null +++ b/src/terminator-main.cpp @@ -0,0 +1,136 @@ +#define _USE_MATH_DEFINES + +#include +#include + +using namespace Rcpp; + +#define rad2deg(rad) ((rad * 180.0) / M_PI) +#define deg2rad(deg) ((deg * M_PI) / 180.0) +#define get_julian(time) (((double)time / 86400.0) + 2440587.5) +#define get_gmst(j_day) (fmod((18.697374558 + 24.06570982441908 * (j_day - 2451545.0)), 24.0)) + +inline std::vector < double > sun_ecliptic_position(double j_day) { + + // compute the position of the Sun in ecliptic coordinates days since start of J2000.0 + double n = j_day - 2451545.0; + + // mean longitude of the Sun + double L = fmod((280.460 + 0.9856474 * n), 360.0); + + // mean anomaly of the Sun + double g = fmod((357.528 + 0.9856003 * n), 360.0); + + // ecliptic longitude of Sun + double lambda = L + 1.915 * sin(deg2rad(g)) + 0.02 * sin(2 * deg2rad(g)); + + // distance from Sun in AU + double R = 1.00014 - 0.01671 * cos(deg2rad(g)) - 0.0014 * cos(2 * deg2rad(g)); + + std::vector< double > ret(2); + ret[0] = lambda; + ret[1] = R; + + return(ret); + +} + +inline double ecliptic_obliquity(double j_day) { + + double n = j_day - 2451545.0; + + // Julian centuries since J2000.0 + double T = n / 36525.0; + + // compute epsilon + return(23.43929111 - + T * (46.836769 / 3600.0 + - T * (0.0001831 / 3600.0 + + T * (0.00200340 / 3600.0 + - T * (0.576e-6 / 3600.0 + - T * 4.34e-8 / 3600.0))))); + +} + +std::vector< double > sun_equatorial_position(double lng, double obliq) { + + double alpha = rad2deg(atan(cos(deg2rad(obliq)) * tan(deg2rad(lng)))); + double delta = rad2deg(asin(sin(deg2rad(obliq)) * sin(deg2rad(lng)))); + + double lQuadrant = floor(lng / 90.0) * 90.0; + double raQuadrant = floor(alpha / 90.0) * 90.0; + + alpha = alpha + (lQuadrant - raQuadrant); + + std::vector< double > ret(2); + + ret[0] = alpha; + ret[1] = delta; + + return(ret); + +} + +inline double hour_angle(double lng, std::vector< double >sun_pos, double gst) { + return((gst + lng / 15.0) * 15.0 - sun_pos[0]); +} + +inline double longitude(double ha, std::vector< double >sun_pos) { + return(rad2deg(atan(-cos(deg2rad(ha)) / tan(deg2rad(sun_pos[1]))))); +} + +//' Compute a single termiantor band +//' +//' Returns a dataframe of latitude and longitude for the line that separates illuminated +//' day and dark night for any given time +//' +//' @md +//' @param time time (numeric from `POSIXct`) for the computation (bands are time-dependent) +//' @param from,to,by latitude sequence setup +//' @return data frame +//' @references , +//' +//' @export +// [[Rcpp::export]] +DataFrame terminator(int time, double from = -180, double to = 180, double by = 0.1) { + + // calculate latitude and longitude of terminator within specified range using time (in POSIXct format, e.g. `Sys.time()`) + double j_day = get_julian(time); + + double gst = get_gmst(j_day); + + std::vector< double > sunEclPos = sun_ecliptic_position(j_day); + + double eclObliq = ecliptic_obliquity(j_day); + + std::vector< double > sunEqPos = sun_equatorial_position(sunEclPos[0], eclObliq); + + std::vector< double > out_lat, out_lon; + + out_lat.reserve(4000); + out_lon.reserve(4000); + + int n=0; + + for (double i=from; i<=to; i+=by) { + n += 1; + out_lat.push_back(i); + out_lon.push_back( + longitude( + hour_angle(i, sunEqPos, gst), + sunEqPos + ) + ); + } + + out_lat.resize(n); + out_lon.resize(n); + + return( + DataFrame::create( + Named("lat") = out_lat, + Named("lon") = out_lon + ) + ); + +} \ No newline at end of file diff --git a/terminator-animation.gif b/terminator-animation.gif new file mode 100644 index 0000000..06fbb0b Binary files /dev/null and b/terminator-animation.gif differ diff --git a/terminator.Rproj b/terminator.Rproj new file mode 100644 index 0000000..446d9e1 --- /dev/null +++ b/terminator.Rproj @@ -0,0 +1,21 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX + +StripTrailingWhitespace: Yes + +BuildType: Package +PackageUseDevtools: Yes +PackageInstallArgs: --no-multiarch --with-keep.source +PackageBuildArgs: --resave-data +PackageRoxygenize: rd,collate,namespace diff --git a/tests/test-all.R b/tests/test-all.R new file mode 100644 index 0000000..a037cfa --- /dev/null +++ b/tests/test-all.R @@ -0,0 +1,2 @@ +library(testthat) +test_check("terminator") diff --git a/tests/testthat/test-terminator.R b/tests/testthat/test-terminator.R new file mode 100644 index 0000000..ab6f62f --- /dev/null +++ b/tests/testthat/test-terminator.R @@ -0,0 +1,6 @@ +context("basic functionality") +test_that("we can do something", { + + #expect_that(some_function(), is_a("data.frame")) + +})