magick animation example

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 <https://github.com/joergdietrich/Leaflet.Terminator/blob/master/L.Terminator.js>,
 #'     <https://github.com/JoGall/terminator/blob/master/terminator.R>

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}

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()
 

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 <https://github.com/joergdietrich/Leaflet.Terminator/blob/master/L.Terminator.js>,
 //'     <https://github.com/JoGall/terminator/blob/master/terminator.R>
@@ -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
       )