diff --git a/R/RcppExports.R b/R/RcppExports.R index 479163e..30f8c60 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -3,12 +3,12 @@ #' Compute a single termiantor band #' -#' Returns a dataframe of latitude and longitude for the line that separates illuminated +#' Returns a dataframe of longitude and latitude 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 +#' @param from,to,by longitude sequence setup #' @return data frame #' @references , #' diff --git a/README.Rmd b/README.Rmd index 86d3dab..913a77f 100644 --- a/README.Rmd +++ b/README.Rmd @@ -60,6 +60,66 @@ mb ## rcpp 493.219 642.67 752.3545 747.8535 828.3085 1044.269 100 ``` +Rather than use the built-in map data, we'll use an `rnaturalearth` so we don't need to mess with extended longitudes. + +```{r wintri, fig.width=10} +library(ggalt) +library(rnaturalearth) # ropensci/rnaturalearth + +world_map <- fortify(countries110) + +ggplot() + + geom_cartogram( + data=world_map, map=world_map, aes(x=long, y=lat, map_id=id), + fill=NA, color="#2b2b2b" + ) -> gg + +for (i in 0:23) { + + gg <- gg + geom_line( + data=terminator(as.integer((as.POSIXct(Sys.Date()) + (60*60*i))), -180, 180, 0.1), + aes(lon, lat), color="blue" + ) +} + +gg + coord_proj("+proj=wintri") + ggthemes::theme_map() +``` + +Animation with `magick`: + +```{r magick, eval=FALSE} +library(magick) +library(rnaturalearth) # ropensci/rnaturalearth +library(tidyverse) + +world_map <- fortify(countries110, region="name") +world_map <- filter(world_map, id != "Antarctica") + +x <- image_graph(width=1000*2, height=500*2, res=144) +pb <- progress_estimated(24) +for (i in 0:23) { + pb$tick()$print() + ggplot() + + geom_cartogram( + data=world_map, map=world_map, aes(x=long, y=lat, map_id=id), + fill=NA, color="#2b2b2b", size=0.125 + ) + + geom_ribbon( + data=terminator(as.integer((as.POSIXct(Sys.Date()) + (60*60*(i)))), -180, 180, 0.1), + aes(lon, ymin=lat, ymax=90), fill="lightslategray", alpha=1/2 + ) + + scale_x_continuous(limits=c(-180, 180)) + + coord_quickmap() + + ggthemes::theme_map() -> gg + print(gg) +} +dev.off() +x <- image_animate(x) +image_write(x, "magick-preview.gif") +``` + +![](magick-preview.gif) + Using Joe's animation example: ```{r anim, eval=FALSE} diff --git a/README.md b/README.md index e06f2fb..f885c01 100644 --- a/README.md +++ b/README.md @@ -59,6 +59,73 @@ mb ## rcpp 493.219 642.67 752.3545 747.8535 828.3085 1044.269 100 ``` +Rather than use the built-in map data, we’ll use an `rnaturalearth` so +we don’t need to mess with extended longitudes. + +``` r +library(ggalt) +library(rnaturalearth) # ropensci/rnaturalearth + +world_map <- fortify(countries110) +``` + + ## Regions defined for each Polygons + +``` r +ggplot() + + geom_cartogram( + data=world_map, map=world_map, aes(x=long, y=lat, map_id=id), + fill=NA, color="#2b2b2b" + ) -> gg + +for (i in 0:23) { + + gg <- gg + geom_line( + data=terminator(as.integer((as.POSIXct(Sys.Date()) + (60*60*i))), -180, 180, 0.1), + aes(lon, lat), color="blue" + ) +} + +gg + coord_proj("+proj=wintri") + ggthemes::theme_map() +``` + +![](README_files/figure-gfm/wintri-1.png) + +Animation with `magick`: + +``` r +library(magick) +library(rnaturalearth) # ropensci/rnaturalearth +library(tidyverse) + +world_map <- fortify(countries110, region="name") +world_map <- filter(world_map, id != "Antarctica") + +x <- image_graph(width=1000*2, height=500*2, res=144) +pb <- progress_estimated(24) +for (i in 0:23) { + pb$tick()$print() + ggplot() + + geom_cartogram( + data=world_map, map=world_map, aes(x=long, y=lat, map_id=id), + fill=NA, color="#2b2b2b", size=0.125 + ) + + geom_ribbon( + data=terminator(as.integer((as.POSIXct(Sys.Date()) + (60*60*(i)))), -180, 180, 0.1), + aes(lon, ymin=lat, ymax=90), fill="lightslategray", alpha=1/2 + ) + + scale_x_continuous(limits=c(-180, 180)) + + coord_quickmap() + + ggthemes::theme_map() -> gg + print(gg) +} +dev.off() +x <- image_animate(x) +image_write(x, "magick-preview.gif") +``` + +![](magick-preview.gif) + Using Joe’s animation example: ``` r @@ -66,7 +133,7 @@ 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) + + geom_ribbon(aes(lon, ymax = lat), ymin = 90, alpha = 0.2) + coord_equal(xlim = c(-180, 190), ylim = c(-58, 85), expand = 0) + ggthemes::theme_map() diff --git a/README_files/figure-gfm/unnamed-chunk-4-1.png b/README_files/figure-gfm/unnamed-chunk-4-1.png new file mode 100644 index 0000000..033397a Binary files /dev/null and b/README_files/figure-gfm/unnamed-chunk-4-1.png differ diff --git a/README_files/figure-gfm/wintri-1.png b/README_files/figure-gfm/wintri-1.png new file mode 100644 index 0000000..f46a65b Binary files /dev/null and b/README_files/figure-gfm/wintri-1.png differ diff --git a/magick-preview.gif b/magick-preview.gif new file mode 100644 index 0000000..cede944 Binary files /dev/null and b/magick-preview.gif differ diff --git a/src/terminator-main.cpp b/src/terminator-main.cpp index 7fe8b27..e0f4181 100644 --- a/src/terminator-main.cpp +++ b/src/terminator-main.cpp @@ -77,22 +77,22 @@ inline std::vector< double > sun_equatorial_position(double lng, double obliq) { } -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 hour_angle(double lat, std::vector< double >sun_pos, double gst) { + return((gst + lat / 15.0) * 15.0 - sun_pos[0]); } -inline double longitude(double ha, std::vector< double >sun_pos) { +inline double latitude(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 +//' Returns a dataframe of longitude and latitude 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 +//' @param from,to,by longitude sequence setup //' @return data frame //' @references , //' @@ -123,7 +123,7 @@ DataFrame terminator(int time, double from = -180, double to = 180, double by = n += 1; out_lon.push_back(i); out_lat.push_back( - longitude( + latitude( hour_angle(i, sunEqPos, gst), sunEqPos )