RSS Feed : http://blog.snap.uaf.edu

RSS Feed from http://blog.snap.uaf.edu

  • Custom images for Shiny dashboard valueBox icons

    The shinydashboard package provides functions like valueBox that conveniently display basic information like summary statistics. In addition to presenting a value and subtitle on a colored background, an icon may be included as well. However, the icon must come from either the Font Awesome or Glyphicon icon libraries and cannot be image files.

    I’ve provided a gist that shows how to achieve the use of custom icons with local image files stored in an app’s www/ directory. It involves overriding a couple functions in shiny and shinydashboard and adding a small bit of custom CSS. Ideally, functionality could be included in future versions of these two packages to allow this in a more robust and complete fashion. But for now, here is a way to do it yourself for value boxes.

    valueboxes1

    The gist above includes the app.R file to run the Custom Icons Shiny app demo and the override.R file which I have it source separately. The gist also includes an icons.R script to generate some statistics and probability themed icons from R. This is interesting and fun on its own. This is not needed to run the app, but the icons are of course needed. If you build the app locally you will have to run this script to generate the images and place them in your www/ folder. The live app demo contains them already.

    I’ve included both light and dark examples using icons that can be used to generally represent a distribution, mean, standard deviation, minimum, maximum, median or interquartile range. Note that app.R adds some custom CSS; it is not sufficient to override the definitions of icon and valueBox alone. I placed it inline for to reduce the number of required files, but it could alternatively be loaded from a file using includeCSS.

    This overriding functionality is only for valueBox widgets and the way in which a local image file is passed to icon is with a named list where the available names are src (required, image file name) and width (optional, defaults to '100%'). It is restrictive but demonstrates potential to use image files in place of icon library options without too much code refactoring.

    Read more »
  • mapmate 0.2.0

    mapmate is under development and blog posts can become outdated quickly. Up to date mapmate documentation and tutorial examples can be found at the official mapmate Github pages.

    mapmate has now been updating from version 0.1.0 to 0.2.0 on Github. The key change is the incorporation of new functions, help docs and code examples focused on network maps, which is a more complex map type not previously covered. The new tutorial content below provides a a couple basic code examples for making network maps with mapmate.

    While mapmate is aimed at still image sequence generation, allowing the user to exert full control over how still image sequences are used to produce animations subsequently (GUI video editor, ffmpeg, ImageMagick, etc) and not directly at animating from R, the examples here include the use of the animation package to help you quickly reproduce some basic animated gifs (as long as you have ImageMagick installed on your system). But the takeaway message is that mapmate now has better support for still image map sequence generation when using the network map type based on great circle arc path traversal.

    If you are new to mapmate there is also the introductory vignette. The content of this post is available on the mapmate Github pages as well.

    To install the package:

    devtools::install_github("leonawicz/mapmate")
    

    Networks

    The save_map function in the mapmate package offers the type="network" map type. This type of map displays networks or pathways defined by overlapping segments traversing along great circle arcs. This map type can be used to display arbitrary line segments as well if such data is provided, but the provided helper functions used here are aimed specifically at drawing great circles.

    Setup

    Obtain great circle arc endpoints

    The first example uses a flat map. Therefore it is important to break lines at the international dateline when preparing the data. In the second example, lines cross the dateline are not split because the segments will be plotted on the globe.

    First, begin with the network data set provided in mapmate is a simple data frame of lon/lat locations of various cities and corresponding weights related to population sizes. It must be expanded to a larger, more complex data frame that contains location pairs and, optionally, distances between locations in each pair. gc_endpoints assists with this by simulating some pairs. The resulting data frame contains endpoints of lines that are defined subsequently.

    library(mapmate)
    library(dplyr)
    
    set.seed(192)
    data(network)
    network
    
    distFun <- function(x) 1 - x/max(x)  # simple inverse distance weighting
    endpoints <- gc_endpoints(network, "lon", "lat")
    endpoints
    

    Obtain great circle arcs

    Next, we sample based on a combination of weights. This is an arbitrary and optional step. The example is given primarily to make the data set used in these examples smaller in size. More importantly, the gc_arcs helper function is used to further expand the endpoints data frame to one containing sequences of points describing great circle arcs instead of only their endpoints.

    The default number of points added between the endpoints of each great circle arc is n=50, but this can be changed and can also differ for each arc (e.g., based on distance between points). As noted, arcs are filled out planning for use with both flat maps and a globes. A group column is used to identify distinct arcs for plotting.

    # take a weighted sample, e.g., favoring larger averaged populations and
    # shorter distances
    endpoints <- mutate(endpoints, Dist_wts = distFun(Dist))
    endpoints <- sample_n(endpoints, 500, replace = TRUE, weight = (Pop_wts0 + Pop_wts1)/2 + 
        Dist_wts)
    
    # expand data frame from endpoints to arcs, each composed of a sequence of
    # points
    arcs_flat <- gc_arcs(endpoints, "lon0", "lat0", "lon1", "lat1", breakAtDateLine = TRUE)
    arcs_globe <- gc_arcs(endpoints, "lon0", "lat0", "lon1", "lat1")
    arcs_globe
    

    Obtain great circle arc path sequences

    Finally, gc_paths is used to further expand the great circle arcs data frame into one that contains sequences of ordered segments along each arc. The segments can vary in length between distinct arcs and can overlap one another within an arc. The group argument is required to identify distinct arcs in the input data frame. size is required to set an upper limit on the number of points constituting an arc segment; the actual length of the segments is chosen randomly and uniformly between 2 and size. Other arguments are optional.

    paths_flat <- gc_paths(arcs_flat, "group", size = 5)
    paths_globe <- gc_paths(arcs_globe, "group", size = 5)
    paths_globe
    

    Flat map network animation

    The direction of arc traversal can also be controlled by the direction argument if desired. This is useful for simulations if the input data are not yet randomized or if the directions simply need to be reversed. First some setup:

    n <- max(paths_flat$id)
    png.args <- list(width = 600, height = 300, bg = "black")
    clrs <- c("#1E90FF50", "#FFFFFF50", "#FFFFFF", "#1E90FF75")
    ylm <- range(paths_flat$lat)  # trimming empty southern map region
    

    Typically, I would leave the default arguments save.plot=TRUE and return.plot=FALSE, but here they are reversed for the purposes of the example. This returns a list of ggplot objects rather than saving png files. Note that I still included the png.args argument even though width and height will be discarded because save_map will take the background color specified by bg intended for png files and use it in the ggplot2 theme applied to the returned plots.

    The default background is transparent so I need to include this here since I have changed it to black. I have changed it to black in this example because I am making a reproducible gif for simplicity rather than doing any layering of separate image sequences in a video editor. Below, saveGIF from the animation package is used to make a simple gif from the plot sequence produced by save_map by looping over the list of returned plots.

    In summary, this is all a somewhat convoluted scenario to show you a short animated gif representing plots made by save_map. The intended use case for save_map is to simply export the sequence of png files and the user can do whatever they wish with those files subsequently. If literally all you want to do is make a short, simple animated gif of custom plots using a small amount of data, just use the animation package. You do not need mapmate for that. Also, if you are using animation in general, as with the code below, it is dependent on ImageMagick, which you will also have to install.

    gglist <- save_seq(paths_flat, id = "id", n.frames = n, ortho = FALSE, type = "network", 
        ylim = ylm, suffix = "2D", png.args = png.args, save.plot = FALSE, return.plot = TRUE)
    
    library(animation)
    # you may need to specify a different path on your Windows machine you may
    # also need to ensure convert.exe is part of your particular installation
    ani.options(convert = "C:/Program Files/ImageMagick-7.0.3-Q16/convert.exe")
    saveGIF(for (i in seq_along(gglist)) print(gglist[[i]]), "network2D.gif", interval = 1/20, 
        ani.width = 600, ani.height = 300)
    

    network2d

    This animation-dependent example and the animated gif are not meant to be distractions from the purpose of this package. Despite the example gif, any code here related to animation is beyond the scope of this tutorial. If you have trouble running it, see the animation documentation and just do the following example instead, which is the way mapmate is meant to be used:

    save_seq(paths_flat, id = "id", n.frames = n, ortho = FALSE, type = "network", 
        suffix = "2D", png.args = png.args)
    # Next, do whatever you want with the files, such as import them to a video
    # editing program
    

    Globe network animation

    Here is an example plotting network paths along great circle arcs on the globe. Remember that we use the other data set, which was generated with the default breakAtDateline=FALSE in gc_arcs. As a side note, if you redo the above example for flat maps using the unbroken data, paths_globe, you will see why the arc segments are handled differently when preparing data for flat maps vs. for globes.

    n <- max(paths_flat$id)
    png.args <- list(width = 600, height = 600, bg = "black")
    clrs <- c("#FFFFFF", "#FF450050", "#FF4500", "#FFFFFF50")
    
    gglist <- save_seq(paths_globe, id = "id", n.frames = n, col = clrs, type = "network", 
        pt.size = c(1, 1, 3, 2), suffix = "3D", png.args = png.args, save.plot = FALSE, 
        return.plot = TRUE)
    
    library(animation)
    ani.options(convert = "C:/Program Files/ImageMagick-7.0.3-Q16/convert.exe")
    saveGIF(for (i in seq_along(gglist)) print(gglist[[i]]), "network3D.gif", interval = 1/20, 
        ani.width = 600, ani.height = 600)
    

    network3d

    Again, normal usage is to just do the following and then use the saved still image sequence with full user control for whatever you like:

    save_seq(paths_globe, id = "id", n.frames = n, col = clrs, type = "network", 
        pt.size = c(1, 1, 3, 2), suffix = "3D", png.args = png.args)
    

    Here are a couple more advanced animations based on still image sequences of great circle arc network maps produced using mapmate code:

    mapmate 0.2.0 (Release date: 2016-11-15)

    • Added functions to assist with network maps: gc_endpoints, gc_arcs, and gc_paths.
    • Added help documentation and runnable examples of the network-related functions.
    • Added network data set to package.
    • Added unit tests for network-related functions.
    • Added tutorial/examples for network maps to the package Github pages.
    • Included simple animated gif examples in above page, piggybacked on animation package.
    • Bug fixes.
    Read more »
  • mapmate 0.1.0

    mapmate has now been updating from version 0.0.2 to 0.1.0 on Github. The biggest addition is a number of plotting options for making different kinds of maps. The new tutorial content below provides a number of code examples for making a variety of maps and also highlights current limitations associated with certain map types and settings.

    Since the code snippets below generate a lot of different plots, the plots are not contained directly in this blog post. See the full tutorial page to see everything. It will also be easier and cleaner to review the examples from there. There is also the accompanying introductory vignette. All can be accessed from the mapmate Github pages.

    bkgd

    To install the package:

    devtools::install_github("leonawicz/mapmate")
    

    Introduction and motivation

    The mapmate package is used for map- and globe-based data animation pre-production. Specifically, mapmate functions are used to generate and save to disk a series of map graphics that make up a still image sequence,
    which can then be used in video editing and rendering software of the user’s choice. This package does not make simple animations directly within R, which can be done with packages like animation. mapmate is more specific to maps, hence the name, and particularly suited to 3D globe plots of the Earth.

    The motivation for the package is to assist with heavy data processing and image production for large, complex data animation videos that can require a significant amount of computing resources to create. It is best suited to use on a Linux server with a significant number of CPU cores and a significant amount of RAM. It is all about video pre-production. Still image sequences created by R are meant to be rendered into video subsequently by dedicated video editing software. It makes the most sense to leverage the power of both programs for what they are best designed for: R for data analysis and plotting and something else for high quality video rendering.

    As mentioned, there are other solutions for producing animations completely from R with only implicit external software dependencies (e.g., using the rgl package or relying on ImageMagick convert). The only restriction is that you have to want to create an animation that is small and simple enough for other solutions to remain realistic options. At some point the size of data or the images or the total number of frame or some other limitation may creep in and stop you if your animation is just simply too big or complex. At that point, this package has the value of being dedicated to generating still image sequences to pass to another program. If you ever make videos which require the use of video editing software, then you will know at what point it makes sense to just make still image sequences that you can drop on your timeline to finish your video project rather than trying to force R to do the entire job itself.

    Nevertheless, while the strength of R for pre-production lies in its data analysis, processing and graphing capabilities, there are some limitations to how well save_map in the mapmate package will draw certain types of maps of certain types of data sets. The biggest limitation is in attempting to draw filled in polygons on a map using an orthographic projection; the 3D Earth globe view.

    Below is a series of usage examples of save_map applied to example data included in the package. The examples show how to make different maps based on the different arguments that can be passed to save_map. They also show how some of the options for map types included save_map are primarily included to reveal their current limitations in drawing maps. Points, paths and tiles tend to work well, as do contour lines for the most part, but polygons often display serious clipping issues.

    This package depends heavily on ggplot2 for plotting. While it is easy enough to write a function that returns an ordered set of lon/lat coordinates defining the circle that represents the edge the hemisphere in view for orthographically projected data, I have not yet found a way in this context to use that circle to essentially reclip the polygons. If I could do so, then I could potentially apply the reclipping inside save_map when ortho=TRUE prior to the reprojection done by ggplot2::coord_map. But not so easily, if the data are not actually projected to the orthographic projection until passed to and as done so by coord_map. The reliance on coord_map could cripple my ability to insert a reclipping step where I need to, but I haven’t delved into it yet since my main use of the package is not for polygons.

    Lastly, since these examples are simply intended to display what kind of map outputs are made given the various arguments that can be passed to save_map, the plot or image frame number/ID value never goes beyond 1 and save_map is not called iteratively. For code related to typical map sequence generation see the introductory vignette for examples: browseVignettes(package="mapmate") or view it here.

    Setup

    First some setup using the package data:

    library(mapmate)
    library(dplyr)
    library(RColorBrewer)
    pal <- rev(brewer.pal(11, "RdYlBu"))
    
    data(annualtemps)
    data(borders)
    data(bathymetry)
    
    id <- "frameID"
    temps <- mutate(annualtemps, frameID = Year - min(Year) + 1) %>% filter(frameID == 
        1)  # subset to first frame
    brdrs <- mutate(borders, frameID = 1)
    bath <- mutate(bathymetry, frameID = 1)
    

    Next for the plots. Note that the aptly named save_map is not typically used for returning plot objects but rather for saving them to disk, hence why every call to save_map in this set of examples must explicitly include the non-default save.plot and return.plot arguments.

    Flat world maps vs 3D Earth

    Use ortho=TRUE (default) or ortho=FALSE.

    save_map(temps, id = id, ortho = FALSE, col = "dodgerblue", type = "points", 
        save.plot = FALSE, return.plot = TRUE)
    save_map(temps, id = id, col = "#FF4500", type = "points", save.plot = FALSE, 
        return.plot = TRUE)
    

    Simple points

    The previous examples use point data. To confirm this works with points everywhere, depending on the resolution of the source data, here is another example using the bathymetry surface with type="points".

    save_map(bath, id = id, type = "points", save.plot = FALSE, return.plot = TRUE)
    

    Spatial/polygon boundary lines

    type="maplines" relies on ggplot2::geom_path internally, which is why it is suited to drawing ordered polygon border lines and why the original borders data frame contains the necessary additional variables. While the prepared data frames used in this set of examples are provided by the package, typically you would generate these yourself from your data which were originally in another form such as SpatialPolygonsDataFrame or RasterLayer objects. A common way to do this is by using ggplot2::fortify or, alternatively, duplicated functionality in the broom package. See those for examples or documentation.

    save_map(brdrs, id = id, type = "maplines", save.plot = FALSE, return.plot = TRUE)
    

    While we’re here, let’s also change the view using lon, lat, and rotation.axis.

    save_map(brdrs, id = id, lon = -70, lat = 40, rotation.axis = 0, type = "maplines", 
        save.plot = FALSE, return.plot = TRUE)
    

    Rasterized tiles

    In the previous examples there was no need to provide a third variable, a data variable. The plots are based on spatial locations only, whether points or paths. For gridded data, a z variable is needed and is specified with the z.name argument.

    save_map(bath, z.name = "z", id = id, col = pal, type = "maptiles", save.plot = FALSE, 
        return.plot = TRUE)
    

    Density maps and contour lines

    type="density" combines both 2D density and 3D contour mapping under a single heading for convenience and based on the default map appearance, but it is important to remember the distinction. The previous example using type="maptiles" is actually a special case and matches the plot output associated with the default arguments that are relevant to type="density".

    Filled contour-style density using tiles

    Default arguments include contour="none" and density.geom="tile".

    save_map(bath, z.name = "z", id = id, col = pal, type = "density", save.plot = FALSE, 
        return.plot = TRUE)
    

    Explicit contour lines can also be overlaid with the base map type.

    save_map(bath, z.name = "z", id = id, col = pal, type = "density", contour = "overlay", 
        save.plot = FALSE, return.plot = TRUE)
    

    Contour lines can also be drawn without including the base layer indicated by type.

    save_map(bath, z.name = "z", id = id, col = pal, type = "density", contour = "only", 
        save.plot = FALSE, return.plot = TRUE)
    

    Contour lines with points

    As a brief aside, these options also work for type="points". Remember that these contour lines are based on the spatial density of lon/lat coordinates, unlike the previous contour lines examples that make use of the data variable.

    save_map(temps, id = id, col = "red", type = "points", contour = "overlay", 
        save.plot = FALSE, return.plot = TRUE)
    save_map(temps, id = id, col = "blue", type = "points", contour = "only", save.plot = FALSE, 
        return.plot = TRUE)
    

    For type="points", the color col applies to both points and contour lines, but for type="density" col applies to the density map. Contour line overlays will remain as black lines. The exception is when contour="only", in which case col may still be a vector defining a color palette that will be applied to the contour lines in place of the underlying filled density map.

    Density maps of locations (no data variable)

    True density maps based on lon/lat locations rather than filled contours based on three variables can also be plotted by excluding z.name. The first call below includes z.name="z" and the second does not. As expected, these are completely different maps because they use fundamentally different data in different ways.

    save_map(temps, z.name = "z", id = id, col = pal, type = "density", contour = "overlay", 
        save.plot = FALSE, return.plot = TRUE)
    save_map(temps, id = id, col = pal, type = "density", contour = "overlay", save.plot = FALSE, 
        return.plot = TRUE)
    

    Of course, this would not be interesting to do when the data frame contains gridded data at all lon/lat locations on the world grid like in bathymetry.

    save_map(bath, id = id, lon = -70, lat = 50, col = pal, type = "density", contour = "overlay", 
        save.plot = FALSE, return.plot = TRUE)
    

    Tiles vs. polygons

    Relying internally on ggplot2::geom_polygon inside save_map can be much faster than ggplot::geom_tile to process the data and make maps. This is an important consideration for generating long sequences of many high-resolution maps. However, the trade off is that polygons may be drawn very poorly with much clipping, the degree of the problem depending on the given data. Of course, if tiles are drawn at an extremely coarse resolution based on the input data, they still won’t look good, and also won’t take as long to draw. It is what you decide to make of it with the data you choose to map. The below examples compare the use of tiles to polygons.

    save_map(temps, z.name = "z", id = id, col = pal, type = "density", contour = "overlay", 
        density.geom = "tile", save.plot = FALSE, return.plot = TRUE)
    save_map(temps, z.name = "z", id = id, col = pal, type = "density", contour = "overlay", 
        density.geom = "polygon", save.plot = FALSE, return.plot = TRUE)
    save_map(bath, z.name = "z", id = id, col = pal, type = "density", contour = "overlay", 
        density.geom = "tile", save.plot = FALSE, return.plot = TRUE)
    save_map(bath, z.name = "z", id = id, col = pal, type = "density", contour = "overlay", 
        density.geom = "polygon", save.plot = FALSE, return.plot = TRUE)
    

    Tiles may also interpolate where possible, but not extrapolate, which may not be desired behavior for a map; e.g., when drawing the continents and interpolating between them but still excluding the poles. It can lead to confusion about what is the underlying data. Yet polygons may still perform much worse, even in flat maps, as shown below.

    save_map(temps, id = id, col = pal, type = "density", contour = "overlay", density.geom = "tile", 
        save.plot = FALSE, return.plot = TRUE)
    save_map(temps, id = id, col = pal, type = "density", contour = "overlay", density.geom = "polygon", 
        save.plot = FALSE, return.plot = TRUE)
    save_map(temps, id = id, col = pal, type = "density", ortho = FALSE, contour = "overlay", 
        density.geom = "tile", save.plot = FALSE, return.plot = TRUE)
    save_map(temps, id = id, col = pal, type = "density", ortho = FALSE, contour = "overlay", 
        density.geom = "polygon", save.plot = FALSE, return.plot = TRUE)
    

    An example from scratch

    This is also apparent when using the simpler options type="polygons" and type="maptiles". The former relies on a SpatialPolygonsDataFrame object as the data source and the later on a RasterLayer. Here examples of each are compared and they are prepped from scratch using the appropriate class objects. Remember that these two map type options are essentially special cases of the density type option.

    # polygons
    library(rworldmap)
    library(rworldxtra)  # required for 'high' resolution map
    library(maptools)  # required for fortify to work
    # also recommend installing rgeos
    
    spdf <- joinCountryData2Map(countryExData, mapResolution = "high")
    [email protected]$id <- rownames([email protected])
    bio <- ggplot2::fortify(spdf, region = "id") %>% left_join(subset([email protected], 
        select = c(id, BIODIVERSITY)), by = "id") %>% mutate(frameID = 1) %>% rename(lon = long)
    
    # raster layer
    library(raster)
    proj4 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +to wgs84=0,0,0"
    z <- "BIODIVERSITY"
    # 1-degree resolution, still somewhat coarse
    r <- raster(extent(-180, 180, -90, 90), nrow = 180, ncol = 360, proj4)
    bio2 <- rasterize(spdf, r, field = z) %>% rasterToPoints %>% tbl_df() %>% setNames(c("lon", 
        "lat", z)) %>% mutate(frameID = 1)
    
    clrs <- c("royalblue", "purple", "orange", "yellow")
    save_map(bio, z.name = z, id = id, lon = -10, lat = 20, col = pal, type = "polygons", 
        save.plot = FALSE, return.plot = TRUE)
    save_map(bio2, z.name = z, id = id, lon = -10, lat = 20, col = pal, type = "maptiles", 
        save.plot = FALSE, return.plot = TRUE)
    

    In conclusion, there are a variety of methods available to draw different types of data on the globe from any arbitrary perspective. Flat maps are also an option but this is not as interesting so it is not focused on. The primary uses of save_map among what is shown here are for plotting points, polygon borders and tiles. The rest are interesting options and can work well enough for specific use cases but are included here mainly to highlight their current limitations.

    The one option not covered here (I will add this in a future update) is the type="network" option, which is used for plotting a time series of maps/globes showing moving great circle arcs. This is my personal favorite use case, but it is also the most complex so I am yet to fold code for that into this package and incorporate great circle arc network examples into the various help documentation, vignettes and tutorials. Something more to look forward to.

    mapmate 0.1.0 (Release date: 2016-11-01)

    • Add type argument options to save_map for points, contour lines, filled contour maps, density/intensity maps
    • Added ability to choose between tile- and polygon-based plotting with type="density".
    • Contour lines can be added on top of density maps, underneath points, or plotted alone as a substitute for either type="points" or type="density" maps.
    • Added a new web page document that provides plenty of examples with visuals on the usage of different save_map arguments
      and more thorough coverage of the limitations associated with using polygons.
    • Remade package website/Github pages to include home page, introductory vignette page, and page showcasing visual examples and current mapping limitations.
    • Added orthographic projection example Shiny app into package website page.
    Read more »
  • mapmate 0.0.2

    mapmate is an R package for map- and globe-based data animation pre-production. Specifically, mapmate functions are used to generate and save to disk a series of global map graphics that make up a still image sequence, which can then be used in video editing and rendering software of the user’s choice. This package does not make simple gif animations directly from R, which can be done with packages like animation. mapmate is more specific to maps, hence the name, and particularly suited to high-resolution png image sequences of 3D globe plots of the Earth.

    This is a development package. Available on Github. See the updated tutorial or vignette for plenty of example use cases and more detailed discussion of existing features and functionality.

    changes in mapmate 0.0.2 (Release date: 2016-10-26)

    • Updated functions, help documentation and examples.
    • Refactored introductory vignette.
    • Additional unit tests included.
    • Change to the behavior of get_lonlat_seq.
    • Additional restrictions imposed on acceptable inputs to functions.
    • Explicit, required id argument replaced the previously assumed presence of a frameID data frame column.
    • No more internal function conversion of id column name to frameID when originally named something else in any of the package functions. Non-standard evaluation is used to avoid dependence on a fixed name.
    • Added convenient wrapper function, save_seq, for maps or time series line plots processed in either series or parallel (Linux-only, via parallel::mclapply).
    • do_projection has been generalized to operate on data frames containing multiple unique plot frame ID values in the specified id column.
    • do_projection has been generalized to output the subsetted data frame with its original columns as before, or with keep=TRUE, the entire input data frame along with an additional boolean inview column.
    • Bug fixes

    The package also contains a simple Shiny app that displays how package function arguments for longitude, latitude, and rotation axis orientation work together to produce a 3D globe display of the Earth from a user-specified perspective.

    devtools::install_github("leonawicz/mapmate")
    library(mapmate)
    runOrtho()
    
    Read more »
  • Animate maps with mapmate: R package for map- and globe-based still image sequences

    [This post is out of date. A newer package version is available on Github and an updated tutorial here]

    Introduction to mapmate

    mapmate (map animate) is an R package for map animation. It is used to generate and save a sequence of plots to disk as a still image sequence intended for later use in data animation production.

    The mapmate package is used for map- and globe-based data animation pre-production. Specifically, mapmate functions are used to generate and save to disk a series of map graphics that make up a still image sequence, which can then be used in video editing and rendering software of the user’s choice. This package does not make simple animations directly within R, which can be done with packages like animation. mapmate is more specific to maps, hence the name, and particularly suited to 3D globe plots of the Earth.

    This introduction covers the following toy examples. Generate a sequence of still frames of:

    • map data for use in a flat map animation.
    • dynamic/temporally changing map data projected onto a static globe (3D Earth)
    • static map data projected onto rotating globe
    • dynamic map data projected onto rotating globe

    Functionality and fine-grain user control of inputs and outputs are limited in the current package version. Other features and functionality will be added in future package versions. While this is a very early developmental stage package, it can also be considered sufficiently complete for its purposes. The primary purpose I had for throwing this package together was to provide a relatively small and simple reproducible example of still image sequence generation similar to those I have used in data animations I’ve shared previously (see bottom of post for examples). The common request I receive is of course for R code.

    While I prefer a world in which all code is shared freely, there are cases where this can be challenging. Some of the R code I write that involves heavy duty data processing is written in the context of particular computing environments, ones which I know almost no one is working in. It’s not the most helpful to share non-reproducible code that no one can use, or code that is written for an environment in which I have access to hundreds of CPUs and terabytes of RAM, if someone wants to duplicate my work on their laptop. That’s why at times I have only been able to provide rough code templates to give people an idea of what I am coding. This package aims to provide some extremely basic but nonetheless explicit, reproducible examples.

    Ultimately, there is no free lunch. But I have done my best here to put together some code that will work in the usual computing environments, which is not parallelized, and which will not be unacceptably slow to process on the small example datasets provided in the package.

    Generate a still image sequence

    Let’s begin with some highly spatially aggregated map data in a data frame with columns lon, lat, and z.

    devtools::install_github("leonawicz/mapmate")
    library(mapmate)
    library(dplyr)
    library(purrr)
    data(annualtemps)
    annualtemps
    

    The save_map function below saves png files to disk as a still image sequence. Set your working directory or be explicit about your file destination. Make sure to obtain the range of data across not only space (a single map) but time (all maps). This is to ensure a constant color to value mapping with your chosen color palette across all maps in the set of png files. For simplicity, set the number of still frames to be equal to the number of unique time points in the data frame, in this one frame for each year.

    Since each point in this data frame represents a coarse grid cell, we use type="maptiles" in the call to save_map. save_map is used for tiles, polygon lines (type="maplines") and great circle arcs (type="network", not covered in this post). One of these three types must be specified and the data must be appropriate for the type.

    Tiles are the simplest, but while we are working with data frames, the data typically would originate from a rasterized object, hence the nice even spacing of the coordinates. Last of all, but important, the data frame must contain a frameID column. Here we will just add a frameID column based on the unique years in the data frame. For these static examples, arguments pertaining specifically to image sequences can be ignored for now (n.period, n.frames, z.range) and since we are not saving files to disk we can also ignore arguments like dir and suffix, but we will return to all these immediately after.

    Initial static 2D and 3D map examples

    library(RColorBrewer)
    pal <- rev(brewer.pal(11, "RdYlBu"))
    temps <- mutate(annualtemps, frameID = Year - min(Year) + 1)
    frame1 <- filter(temps, frameID == 1)  # subset to first frame
    
    save_map(frame1, ortho = FALSE, col = pal, type = "maptiles", save.plot = FALSE, 
        return.plot = TRUE)
    save_map(frame1, col = pal, type = "maptiles", save.plot = FALSE, return.plot = TRUE)
    

    unnamed-chunk-4-1
    unnamed-chunk-4-2

    Flat map sequence

    Dynamic data, static map

    save_map is not typically used in the manner above. The above examples are just to show static plot output, but there is no convenience in using mapmate for that when we can just use ggplot2 or other plotting packages directly. Below are some examples of how save_map is usually called as a convenient still image sequence generator.

    Notice the importance of passing the known full range of the data. The optional filename suffix argument is also used. For a flat map we still do not need to pay attention to the globe-specific arguments (lon, lat, rotation.axis) and those dealing with displaying data while globe rotation occurs across the image sequence can also be ignored (n.period).

    rng <- range(annualtemps$z, na.rm=TRUE)
    n <- length(unique(annualtemps$Year))
    suffix <- "annual_2D"
    temps <- split(temps, temps$frameID)
    walk(temps, ~save_map(.x, ortho=FALSE, col=pal, type="maptiles", suffix=suffix, z.range=rng))
    
    

    Globe sequence

    Dynamic data, static map

    For plotting data on the globe, set ortho=TRUE (default). In this example the globe remains in a fixed position and orientation. The default when passing a single, scalar longitude value to the lon argument is to use it as a starting point for rotation. Therefore, we need to pass our own vector of longitudes defining a custom longitude sequence to override this behavior. Since the globe is to remain in a fixed view, set lon=rep(-70, n). Note that a scalar lat argument does not behave as a starting point for rotation in the perpendicular direction so no explicitly repeating vector is needed.

    Also pay attention to n.period, which defines the number of frames or plot iterations required to complete one globe rotation (the degree increment). n.period can be any amount when providing a scalar lon value. The longitude sequence will propagate to a length of n.period. However, if providing a vector representing a custom path sequence of longitudes to lon, there is no assumption that the sequence is meant to be cyclical (do your own repetition if necessary).

    It could just be a single pass, custom path to change the globe view. In this case, n.period must be equal to the length of the lon and lat vectors or an error will be thrown. This is more a conceptual restriction than anything. With a custom path defined by lon and lat vectors, which may not be at all cyclical, it only makes sense to define the period as the length of the sequence itself.

    suffix <- "annual_3D_fixed"
    walk(temps, ~save_map(.x, lon = rep(-70, n), lat = 50, n.period = n, n.frames = n, 
        col = pal, type = "maptiles", suffix = suffix, z.range = rng))
    
    

    Static data, dynamic map

    Now let the globe rotate, drawing static data on the surface. In this example the nation borders are drawn repeatedly while allowing the Earth to spin through the still image sequence from frame to frame. walk from the purrr package is used to easily duplicate the borders data frame in a list while adding a frameID column that increments over the list.

    It's quite redundant and this is worse if the dataset is large, but this is how save_map currently accepts the data. Fortunately, it would make little difference if the data change from frame to frame, which is the main purpose of save_map. Datasets which do not vary from frame to frame generally serve the role of providing a background layer to other data being plotted, as with the nation borders here and as we will see later with a bathymetric surface.

    In the case of repeating data passed to save_map, it only makes sense when plotting on a globe and when using a perspective that has some kind of periodicity to it. This limits how many times the same data must be redrawn. For example, if the Earth completes one rotation in 90 frames, it does not matter that maps of time series data may be plotted over a course of thousands of frames. When layering those maps on top of a constant background layer on the rotating globe, only 90 plots of the background data need to be generated. Those 90 saved images can be looped in a video editor until they reach the end of the time series maps image sequence.

    The plot type has changed to maplines. z.range can be ignored because it only applies to map tiles (essentially square grid polygon fill coloring) which have values associated with surface pixels. It is irrelevant for polygon outlines as well as for network plots of line segment data. A single color can be used for the lines so the previously used palette for tiles has been replaced.

    data(borders)
    borders <- map(1:n, ~mutate(borders, frameID = .x))
    suffix <- "borders_3D_rotating"
    walk(borders, ~save_map(.x, lon = -70, lat = 50, n.period = 30, n.frames = n, 
        col = "orange", type = "maplines", suffix = suffix))
    
    

    Using one more example, return to the list of annual temperature anomalies data frames used previously. Those were used to show changing data on a fixed-perspective globe plot. Here we plot the first layer, repeatedly, as was just done with the borders dataset, allowing the Earth to rotate. Remember to update the frameID values. The key difference with this example of fixed data and a changing perspective is that providing z.range is crucial to maintaining constant color mapping.

    temps1 <- map(1:n, ~mutate(temps[[1]], frameID = .x))
    rng1 <- range(temps1[[1]]$z, na.rm = TRUE)
    suffix <- "year1_3D_rotating"
    walk(temps1, ~save_map(.x, lon = -70, lat = 50, n.period = 30, n.frames = n, 
        col = pal, type = "maptiles", suffix = suffix, z.range = rng1))
    
    

    Dynamic data, dynamic map

    Finally we have the case of changing data drawn on a map with a changing perspective. This example plots the full time series list of annual temperature anomalies data frames with globe rotation.

    suffix <- "annual_3D_rotating"
    walk(temps, ~save_map(.x, lon = -70, lat = 50, n.period = 30, n.frames = n, 
        col = pal, type = "maptiles", suffix = suffix, z.range = rng))
    
    

    Of course, there is no reason to draw static data on a static flat map or static (e.g., non-rotating) globe because that would yield a still image sequence of a constant image. The utility of the function is in generating a series of plots where either the data changes, the perspective changes, or both.

    Multiple image sequences

    Putting it all together, we can generate three still image sequences to subsequently be layered on top of one another in an animation. Below is an example using a bathymetry surface map as the base layer, the temperature anomalies that will be the middle layer, and the nation borders which will stack on top. This also is why the saved png images have a default transparent background. There is an expectation of eventual layering.

    data(bathymetry)
    bath <- map(1:n, ~mutate(bathymetry, frameID = .x))
    rng_bath <- range(bath[[1]]$z, na.rm = TRUE)
    pal_bath <- c("black", "steelblue4")
    
    walk(bath, ~save_map(.x, n.frames = n, col = pal_bath, type = "maptiles", suffix = "background", 
        z.range = rng_bath))
    walk(borders, ~save_map(.x, n.frames = n, col = "black", type = "maplines", 
        suffix = "foreground"))
    walk(temps, ~save_map(.x, n.frames = n, col = pal, type = "maptiles", suffix = "timeseries", 
        z.range = rng))
    
    

    Parallel processing

    For larger datasets and longer sequences, this is much faster the more it can be parallelized. mapmate is designed to add convenience for making relatively heavy duty animations. The emphasis is on images which will look sharp and lend themselves to smooth frame transitions. They may also do so while displaying a large amount of data. High-resolution images (even larger than 4K) can be used to allow zooming during animations without loss of quality. For all these reasons, processing large amounts of data and generating still image sequences can take a long time depending on how large, long, and complex the desired plots sequences are.

    The toy examples given earlier are not ideal for serious production. The code is provided for simple reproduction and tutorial purposes. The save_map function can be used more powerfully on a Linux server with a large number of CPU cores, and enough RAM to meet the needs of your data. It uses mclapply from the base R parallel package. It will not work on Windows. The choice was made for convenience.

    The example below is like the previous one, but using mclapply. It assumes you have a 32-CPU Linux server node. If you have multiple nodes, you could even go so far as to explore the Rmpi package to link across, say, 10 nodes to yield the power of 320 CPUs. But if you have access to that kind of computing power, you probably already know it and can explore it on your own. It is beyond the scope of this tutorial.

    mclapply(bath, save_map, n.frames = n, col = pal_bath, type = "maptiles", suffix = "background", 
        z.range = rng_bath, mc.cores = 32)
    mclapply(borders, save_map, n.frames = n, col = "orange", type = "maplines", 
        suffix = "foreground", mc.cores = 32)
    mclapply(temps, save_map, n.frames = n, col = pal, type = "maptiles", suffix = "timeseries", 
        z.range = rng, mc.cores = 32)
    
    

    Examples using mapmate

    Historical and projected global temperature anomalies

    Global UAF/SNAP Shiny Apps web traffic

    Flat map great circle animation example

    Read more »
  • R animation: global climate change

    I have posted a new R data animation video. It’s an example animation of modeled historical and projected global temperature change from 1850 – 2100. The data prep, analysis, full processing and generation of all sets of still frames for each layer in the video are done using R.

    Typically an ensemble of models would be used but this video is just to demonstrate a basic animation using one climate model, both with a monthly time series and a monthly 10-year moving average time series. If wondering about the y-axis range, the animation shows anomalies, or delta change, from the climate model’s historical baseline monthly average temperatures using a given climatology window.

    In a later video I will use annual and seasonal averages, which will display a smoother signal than monthly series.

    Read more »
  • Animated great circles on rotating 3D Earth: Example R code

    Making still frames for great circle animations in R

    Matthew Leonawicz

    Here I share R code I used to produce animated great circle arcs on top of a rotating 3D Earth. The code is not entirely reproducible but you should be able to use what is shared here to create your own video frames given your unique data and computing environment and resources.

    The WordPress blog is not the most elegant for displaying lots of code so go to the original full post.

    When I make great circle animations, at the core of the process is always an R function that transforms a series of coordinates describing points along a great circle arc into multiple series of great circle arc segments. The goal is simple: plot a series of line segments, saving each plot as a subsequent still frame, rather than plotting the original entire arc as a single plot. The input is generally a data table (much faster to work with than a data frame if you have a lot of data) with longitude and latitude columns where the coordinates in each row describe a subsequent point along one of my paths. I also use a third column to provide a unique group ID for each path to keep them distinct.

    Before getting to this process, here is some example code of how I formatted my data this way when using the geosphere package. Do not get bogged down in the details here. This is just an example for fuller context regarding my specific animation referenced above and I will not be focusing on it. Your data will be quite different. You will have to arrange it similarly, but obviously it will be in a different context.

    Data prep example

    As you can see below, all you need is a data frame with columns of longitude and latitude. I created a SpatialPoints object from these. In my case I wanted to connect great circle arcs between each of two specific locations and all other locations in my data set. I marked the row indices for these two locations.

    library(dplyr)
    library(geosphere)
    
    load("data.RData")  # a data frame, d, containing 'long' and 'lat' columns
    p <- SpatialPoints(cbind(d$long, d$lat), proj4string = CRS("+proj=longlat +datum=WGS84"))
    idx1 <- 69  # great circles from coords in all other rows to coords in this row
    idx2 <- 2648  # as above
    
    get_paths <- function(x, idx, ...) {
        gcInt <- function(x, x1, x2) {
            x <- gcIntermediate(x[x1, ], x[x2, ], ...)
            if (is.list(x)) {
                x <- x %>% purrr::map2(c(x1, x1 + 0.5), ~data.frame(.x, .y)) %>% 
                    bind_rows %>% setnames(c("long", "lat", "group"))
            } else x <- data.frame(x, x1) %>% setnames(c("long", "lat", "group"))
            x
        }
        purrr::map(setdiff(1:length(x), idx), ~gcInt(x, .x, idx)) %>% bind_rows
    }
    
    paths1 <- get_paths(p, idx1, addStartEnd = TRUE)
    paths2 <- get_paths(p, idx2, addStartEnd = TRUE)
    

    get_paths uses gcIntermediate from the geosphere package to obtain vectors of points defining great circle arcs between my two locations and all others. The resulting data frame has columns, long, lat, and group. I do not break the arcs at the dateline because I intend to draw them on a 3D globe, but if I were making a flat map I would do so. This is why the function handles the case of list output from gcIntermediate and adjusts the group ID by add 0.5 to one group.

    Already more than enough details. You can see what I am going for though. Given my data, I want to end up with columns of longitude and latitude defining great circle arcs and broken out by unique group IDs. I show this because it's highly likely you will want to use geosphere in a similar way even if you won't be connecting points in the same way I am here.

    Transforming great circle arcs into segments

    Okay. You have your data in the right format. Now you want to break up your groups of great circle arcs into a larger number of nested subgroups of great circle arc segments. Let's get some basic prep out of the way.
    Here are the packages I am using.

    The steup

    library(parallel)
    library(gridExtra)
    library(raster)
    library(data.table)
    library(dplyr)
    library(ggplot2)
    
    eb <- element_blank()
    theme_blank <- theme(axis.line = eb, axis.text.x = eb, axis.text.y = eb, axis.ticks = eb, 
        axis.title.x = eb, axis.title.y = eb, legend.position = "none", panel.background = eb, 
        panel.border = eb, panel.grid.major = eb, panel.grid.minor = eb, plot.background = element_rect(colour = "transparent", 
            fill = "transparent"))
    
    world <- map_data("world")
    

    The ggplot2 theme will allow for plotting without all the extraneous stuff like margins and axes and colors I don't want on my Earth. The world map is an aside and has nothing to do with the great circles. I include it here because it was part of my animation. It will be used from plotting the backdrop of nation boundaries on the globe. Having these underneath the great circle animation is helpful for geographic visual orientation.

    It's not clear yet, but the reason for including the raster package is because I also use a rasterized map layer of the earth's surface as the bottom layer underneath the nation boundaries. This is also an aside. It doesn't have anything to do with the great circle animation.

    Yes, I am using R's built-in parallel package. I happen to be do at least the parallelized operations in this project on a Linux server with 32 CPUs and 260 GB RAM. I'm sorry to say but I can't help you if you want to do this on a local Windows pc for example. If you rework the code for a much more restrictive environment, depending on your data, you may find yourself waiting forever for your output. I just do not recommend doing this type of thing outside of a beefy server environment, at least if time is something you value.

    Now, getting to the core of the process for real this time. Here is the actual function I used for the above animation to break great circle arcs into smaller segments.

    The main function

    df_segs <- function(d, seg.size, n.frames, replicates = 1, direction = "fixed") {
        n <- nrow(d)
        if (n < 3) 
            stop("Data not appropriate for this operation.")
        if (seg.size < 3) 
            stop("Segment size too small.")
        z <- round(runif(2, 2, seg.size))
        z[z > n] <- n
        n1 <- ceiling(diff(c((z[1] - z[2]), n))/z[1])
        if (n.frames - n1 < 100) 
            stop("Insufficient frames")
        offset <- sample(0:(n.frames - n1), replicates)
        
        f <- function(k, d, n, n1, z, offset) {
            ind2 <- z[1] * k
            ind1 <- max(ind2 - z[2], 1)
            if (ind2 > n) 
                ind2 <- n
            d <- slice(d, ind1:ind2)
            purrr::map(offset, ~mutate(d, group = ifelse(replicates == 1, group, 
                group + as.numeric(sprintf(".%d", k))), frameID = .x + k)) %>% bind_rows
        }
        
        if (direction == "reverse") 
            d <- mutate(d, long = rev(long), lat = rev(lat))
        if (direction == "random" && rnorm(1) < 0) 
            d <- mutate(d, long = rev(long), lat = rev(lat))
        d <- purrr::map(1:n1, ~f(.x, d, n, n1, z, offset)) %>% bind_rows %>% arrange(group, 
            frameID)
        d
    }
    

    Going through the code, you can see it requires a minimum of 100 desired frames (which would make for quite a short video). You can tweak the maximum number of points in a segment. The minimum is always two. Each segment will vary in length uniformly between the minimum and maximum.

    You can leave the data as is, with direction="fixed". This assumes the ordering of points/rows pertaining to each group is intended to show direction. Segments will be assembled in the same order. Alternatively, you can reverse or even randomize the order if it doesn’t matter for the given data.

    This is just one example function. You can make your own if this one does not generate the type of segments or provide the kind of random variation in segments you would like.

    Let’s go with 900 frames, which will result in about a 30 second video. Recall that in my case I had two different data sets. I combine them at the end of the above code snippet.

    n.frames <- 900
    set.seed(1)
    paths2 <- mutate(paths2, group = group + max(paths1$group))
    paths <- bind_rows(paths1, paths2) %>% split(.$group)
    

    Here I make the new table below. It has a subgroup ID column as well as a frame ID column. Group ID is now a decimal. To the left is the original great circle ID. To the right is the sequential segment ID. The latter is random and not all great circle arcs are broken up into the same number of segments. Some cover more rows of the table than others. They do not match up with the frame IDs and some recycling may be required if the number of desired frames is large enough.

    paths <- mclapply(paths, df_segs, seg.size = 5, n.frames = n.frames, replicates = 1, 
        direction = "random", mc.cores = 32) %>% bind_rows
    

    Alternatively, code like this will perform the same operation without using parallel processing. This may not be terribly problematic if the data set is not too large. I show this for comparison. But the need for parallel will be much greater in later steps.

    paths <- paths %>% split(.$group) %>% purrr::map(~df_segs(.x, 5, n.frames, replicates = 1, 
        direction = "random")) %>% bind_rows
    

    Aside: the background layertile

    Where I got the data

    I use a rasterized bathymetry surface based on a csv file I downloaded by using the marmap package. I won’t discuss that here. Please see the marmap package for details and examples if this is important. It is straightforward to use that package to download a local copy of a subregion of the data (or the whole map) at a specified resolution.

    Projecting background layer to 3D

    The project_to_hemisphere function below (adapted from mapproj here) projects points onto the 3D earth and identifies which ones are within the hemisphere field of view given the centroid focal point. Identifying the half which are out of view in any given frame allows me to toss them out for some added efficiency. This is not actually a slow process. The conversion of coordinates is not complex or demanding. The slow part is drawing the high resolution orthographic projection maps with ggplot2.

    d.bath <- read.csv("marmap_coord_-180;-90;180;90_res_10.csv") %>% data.table %>% 
        setnames(c("long", "lat", "z"))
    r <- raster(extent(-180, 180, -90, 90), res = 1/6)
    projection(r) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
    r <- setValues(r, d.bath$z)
    
    project_to_hemisphere <- function(lat, long, lat0, long0) {
        hold <- cbind(lat, long)
        x <- purrr::map(list(lat, lat0, long - long0), ~.x * pi/180)
        inview <- sin(x[[1]]) * sin(x[[2]]) + cos(x[[1]]) * cos(x[[2]]) * cos(x[[3]]) > 0
        data.table(long = hold[, 2], lat = hold[, 1], inview = inview)
    }
    

    Prepare for plotting

    I want a 120-frame period, which is about four seconds of video, in which time I will complete one rotation of the earth. I define a sequence of longitude values, but this I allow to repeat for the total number of video frames, n.frames, specified earlier. Note that I keep the latitude focus (41 degrees) constant in this animation, but there is no reason this cannot vary as well. In the save_maps function below, there is a hard-coded 23.4-degree fixed orientation, but this is something that can also be made variable.

    While the data we are interested in – the great circle arc segments – will plot sequentially over a series of 900 frames, both the nation boundaries and bathymetry surface map backdrops are constant so I only need to produce enough frames (120) for one rotation of each of those layers.

    I split the path segments table on the frame ID. I will generate still frames in parallel.

    In this animation I had originally chosen to aggregate my 10-minute resolution bathymetry surface from the marmap package because that resolution was a bit too fine and even in parallel it was taking an extremely long time for ggplot2 to draw the orthographic projection background maps. If I understand correctly, it is drawing a ton of tiny, filled polygons. It really grinds to a halt if you have a massive amount of them.

    Some setup

    n.period <- 120
    lon_seq <- rep(seq(0, 360, length.out = n.period + 1)[-(n.period + 1)], length = n.frames)
    lat_seq <- rep(41, length(lon_seq))
    paths <- paths %>% split(.$frameID)
    d.bath.agg <- r %>% aggregate(2) %>% rasterToPoints %>% data.table %>% setnames(c("long", 
        "lat", "z"))
    

    Next, in parallel across frames, I project my rasterized background map cells to the 3D surface and retain only the ones which are in view in a given frame as the earth spins. I use the full range of elevation data to retain a constant color palette mapping as cells at different extreme elevations move in and out of the field of view as the earth rotates. I store the range in z.range. I add a frame ID column. Finally I add a frame ID column for the nation boundaries table. The data is constant, but simply needs to be plotted 120 times in 3-degree shifts.

    d.tiles <- mclapply(1:n.period, function(i, x, lon, lat) {
        left_join(x, project_to_hemisphere(x$lat, x$long, lat[i], lon[i])) %>% filter(inview) %>% 
            dplyr::select(-inview) %>% mutate(frameID = i)
    }, x = d.bath.agg, lat = lat_seq, lon = lon_seq, mc.cores = 32)
    
    z.range <- purrr::map(d.tiles, ~range(.x$z, na.rm = TRUE)) %>% unlist %>% range
    d.world <- purrr::map(1:n.period, ~mutate(world, frameID = .x))
    

    Save plots

    I've showed the most important of two functions for making my great circle animation, df_segs, which I use to transformed a table of grouped rows of great circle arc coordinates into a much larger one with nested, grouped, sequential arc segments. The other critical function is save_maps.

    Here I have generalized it a bit by adding a type argument. It is primarily used for iterating over the data for each frame in a table produced by df_segs (the default type="network"). This saves a map image of sequential great circle arc segments for each frame.

    I added the options, type="maplines" and type="maptiles". The former is for the nation boundaries and the latter is for the rasterized map surface.

    I call the function in parallel three times to cover each type of output. I generate 120 frames (one earth rotation in three-degree increments) of both the nation boundaries and the surface tiles. The latter takes by far the longest amount of time to process, far longer than even the 900 frames of network traffic itself.

    The other main function

    save_maps <- function(x, lon_seq, lat_seq, col = NULL, type = "network", z.range = NULL) {
        if (is.null(col)) 
            col <- switch(type, network = c("#FFFFFF25", "#1E90FF25", "#FFFFFF", 
                "#1E90FF50"), maptiles = c("black", "steelblue4"), maplines = "white")
        i <- x$frameID[1]
        if (type == "network") 
            x.lead <- group_by(x, group) %>% slice(n())
        g <- ggplot(x, aes(long, lat))
        if (type == "maptiles") {
            if (is.null(z.range)) 
                z.range <- range(x$z, na.rm = TRUE)
            g <- ggplot(x, aes(long, lat, fill = z)) + geom_tile() + scale_fill_gradientn(colors = col, 
                limits = z.range)
        } else {
            g <- ggplot(x, aes(long, lat, group = group))
            if (type == "maplines") 
                g <- g + geom_path(colour = col)
            if (type == "network") 
                g <- g + geom_path(colour = col[2]) + geom_path(colour = col[1]) + 
                    geom_point(data = x.lead, colour = col[3], size = 0.6) + geom_point(data = x.lead, 
                    colour = col[4], size = 0.3)
        }
        g <- g + theme_blank + coord_map("ortho", orientation = c(lat_seq[i], lon_seq[i], 
            23.4))
        dir.create(outDir <- file.path("frames", type), recursive = TRUE, showWarnings = FALSE)
        png(sprintf(paste0(outDir, "/", type, "_%03d.png"), i), width = 4 * 1920, 
            height = 4 * 1080, res = 300, bg = "transparent")
        print(g)
        dev.off()
        NULL
    }
    
    mclapply(paths, save_maps, lon_seq, lat_seq, type = "network", mc.cores = 30)
    mclapply(d.world, save_maps, lon_seq, lat_seq, type = "maplines", mc.cores = 30)
    mclapply(d.tiles, save_maps, lon_seq, lat_seq, type = "maptiles", z.range = z.range, 
        mc.cores = 30)
    

    What about the video?

    For any project like this I simply drop each of the three sequences of sequentially numbered png files onto its own timeline track as a still image sequence in a standard video editor. I don't use R for this final step. While I could have plotted everything together on a single sequence of frames, it doesn't make sense to do so even when each sequence has equal length. This is because one sequence may take far longer to plot than another sequence. It is more efficient to make a different image sequence for separate layers and mix them in the editor afterward.

    You might also notice the very high pixel dimensions of the png outputs. I do this in cases where I plan to zoom during an animation. The larger the image, the more I can zoom in without degradation. You can imagine that this generates a lot of data when you have thousands of frames and each image is several megabytes.

    Conclusion: Don't try this at home

    Instead, try something like this on a server with plenty of CPUs and RAM. Otherwise try it with a small data set. Due to the particular computing environment I used and the resources needed to process the frames in an efficient and timely fashion, it is impossible to share this example animation code in a completely reproducible fashion. However, you should be able to use df_segs or some alteration of it with a properly formatted input table to generate the type of table you would want to pass subsequently to save_maps (or some alteration of that as well).

    That is the gist of this example and should be enough to get you started. Ultimately, you will have to adapt any code here to your unique data, its size and complexity, what exactly you want to do with it, and your available computing resources.

    Read more »
  • R animation: global SNAP Shiny app users

    For the blog readers, just a quick heads up that I have posted a new R data animation to YouTube. A complete post will follow, but for now here is the video. It displays the social network of SNAP Shiny app users over about the past year and a half using great circle trajectories on a rotating 3D Earth. It’s best in 1080p, but still somewhat degraded for streaming. I’ll post the raw source video later as well, which is crystal clear.

    I used geolocation data from Google Analytics. I routed all the traffic randomly through either Fairbanks, AK or Denver, CO to complete the network since those two places best describe where my apps come from. It’s great to see how many people connect to the apps and from where; too many people to plot as a single static graphic without some of the essence of the data being lost.

    As usual, the still frames of the rotating globe, its surface texture, country boundaries and great circle paths are all made in R. It’s basically a series of plots made with ggplot2.

    Read more »
  • R Shiny leaflet: using observers

    Using leaflet in R is fairly simple and there are plenty of online resources and examples available already. See the RStudio tutorial to get started if you haven’t already. This post presents a series of examples which build upon each other.

    ex_leaflet

    The code displays both point data and raster data on leaflet maps with Shiny app integration. The focus here is to show the utility of event observation on the server side for integrating the behavior of browser elements, namely, the leaflet map.

    In the course of these examples, I’ll also use a modal from the shinyBS package and the maps will be full size in the browser window with UI widgets laid on top using absolute panels. The full post contains interactive embedded Shiny widgets. See the full post here. It includes complete source code for these and other related apps using leaflet.

    The app does not do much in this version, but it provides a clear view of how to use observers in Shiny apps to ensure that various inputs and outputs remain mutually well behaved. Here I use observeEvent three distinct times to tie together the reactive behavior of instances of selectInput and leafletOutput. The entire ui.R at this point is:

    ui <- bootstrapPage(
      tags$style(type="text/css", "html, body {width:100%;height:100%}"),
      leafletOutput("Map", width="100%", height="100%"),
      absolutePanel(top=10, right=10,
        selectInput("location", "Community", c("", locs$loc), selected=""),
        conditionalPanel("input.location !== null && input.location !== ''",
          actionButton("button_plot_and_table", "View Plot/Table", class="btn-block"))
      )
    )
    

    In server.R note the three observers. They constitute the entire code block except for the initial, relatively simple call to renderLeaflet. Comments mark the role performed by each call to observeEvent, which:

    • Update the leafletOutput map when clicked
    • Update the selectInput when the map is clicked
    • Update the leafletOutput when the selectInput changes
    server <- function(input, output, session) {
      acm_defaults <- function(map, x, y) addCircleMarkers(map, x, y, radius=6, color="black", fillColor="orange", fillOpacity=1, opacity=1, weight=2, stroke=TRUE, layerId="Selected")
    
      output$Map <- renderLeaflet({
        leaflet() %>% setView(lon, lat, 4) %>% addTiles() %>%
          addCircleMarkers(data=locs, radius=6, color="black", stroke=FALSE, fillOpacity=0.5, group="locations", layerId = ~loc)
      })
    
      observeEvent(input$Map_marker_click, { # update the map markers and view on map clicks
        p <- input$Map_marker_click
        proxy <- leafletProxy("Map")
        if(p$id=="Selected"){
          proxy %>% removeMarker(layerId="Selected")
        } else {
          proxy %>% setView(lng=p$lng, lat=p$lat, input$Map_zoom) %>% acm_defaults(p$lng, p$lat)
        }
      })
    
      observeEvent(input$Map_marker_click, { # update the location selectInput on map clicks
        p <- input$Map_marker_click
        if(!is.null(p$id)){
          if(is.null(input$location) || input$location!=p$id) updateSelectInput(session, "location", selected=p$id)
        }
      })
    
      observeEvent(input$location, { # update the map markers and view on location selectInput changes
        p <- input$Map_marker_click
        p2 <- subset(locs, loc==input$location)
        proxy <- leafletProxy("Map")
        if(nrow(p2)==0){
          proxy %>% removeMarker(layerId="Selected")
        } else if(length(p$id) && input$location!=p$id){
          proxy %>% setView(lng=p2$lon, lat=p2$lat, input$Map_zoom) %>% acm_defaults(p2$lon, p2$lat)
        } else if(!length(p$id)){
          proxy %>% setView(lng=p2$lon, lat=p2$lat, input$Map_zoom) %>% acm_defaults(p2$lon, p2$lat)
        }
      })
    
    }
    
    shinyApp(ui, server)
    

    Care must be taken to ensure each observer does not update when it shouldn’t. One can easily trigger another back and forth in this context. Fortunately it is not too difficult to make the three observeEvent calls work well together. See the full post to use the app and continue with the interactive examples.

    Read more »
  • Upload shapefile to R Shiny app to extract leaflet map data

    In this post I share an R Shiny app which uses the leaftlet package for interactive maps. This app differs from prior apps I’ve made featuring leaflet maps. First, it displays rasterized map data rather than just point layers. Also, longitude and latitude sliders in the browser allow for cropping the map. Additionally, the user can upload a shapefile to crop and mask the rasterized data overlays in the leaflet window to the specific spatial data they wish to work with, and then extract and download that data.

    nwtapp1

    Explore the app here

    Features

    • Crop map extent using latitude and longitude sliders.
    • Upload a custom shapefile in the upload window to mask map layer.
    • Visually confirm your original shapefile and reprojected (if applicable) shapefile on the map look right.
    • The uploader compares your shapefile with the app data to ensure validity.
    • Shapefiles may contain points, lines, or polygons.
    • Combine crop and mask operations.
    • Use the color swatches to select from common color palettes or create your own.
    • Download the displayed map as a geotiff.
    • Use a modal to maintain clear separation between exploration in space at single time points (maps) and in time at single locations (time series charts).
    • Disable tooltips if you don’t need them.

    Another new feature in this app, and relatively new to the Shiny package, is the incorporation of a Shiny module. Shiny modules are handy ways to abstract, externalize and reuse Shiny app code. This is useful for keeping apps tidy when they grow large and complex, but also when the same code can be recycled across multiple apps.

    In this case I wrote a Shiny module that provides an interactive file uploader which appears in a modal window popup (shinyBS package) when a button is clicked. The file uploader is specifically geared toward accepting shapefiles. It displays previews of the original file and of the final shapefile as projected (if necessary) and overlain on top of a map. The module does some nominal validation of the shapefile to make sure the user has uploaded, for example, something which actually overlaps with the data in the app shown on the leaflet map. Once uploaded, the user can apply their shapefile as a mask to the app data.

    The color picker (shinyJS package) allows for creating your own color palette if not using a predefined palette.

    This app is restricted to climate data over Northwest Territories (NT), Canada. Some of the vetting of user uploaded shapefiles, like comparing spatial extents of a user-provided shapefile with this political boundary, would not have to occur if the app data had a worldwide extent. If you upload a shapefile of Australia, it will simply inform you that your shapefile does not overlap the available data. So if you want to play around with this feature, I have uploaded three example shapefiles for your convenience. Download and unzip the file. It contains one shapefile each of points, lines, and polygons which all intersect NT. The polygon also contains a hole. …And I at least know that these three shapefiles work! But please do add a Github issue if you find one that causes problems.

    Shapefiles of points, lines or polygons are all permitted and handled by the app appropriately. For now I only permit downloading rasters as geotiffs. I can see cases where a non-spatial format may be desired, such as a csv file containing columns of longitude, latitude and data at least in the case of points and perhaps lines, or when data refers to results of an aggregate function that is applied to an entire polygon or to named sub-polygons.

    Lastly, in addition to the modal popup window used to display the file uploader, a modal is also used for displaying time series plots pertaining to the currently selected leaflet map marker. Selecting a community on the map gives the user the option to open this modal and see location-specific data through time, which cannot be seen in the leaflet map itself.

    While the leaflet maps show single decadal time slices of climate information spatially, a number of select NT communities have annual-resolution climate data contained in the app which can be viewed distinct from the leaflet map side of things. Of course, if a user uploads an arbitrary shapefile of points, they can only extract and download from what is shown in the decadal map; there would be no similar data sets available for every possible map pixel, hence why user-uploaded points which are added to the leaflet map as map markers are not clickable like the points which come included in the app as a special data set.

    Full source code is available on Github.

    Read more »

Copyright Use-R.com 2012 - 2016 ©