Lucinda Roberts | March, 2024

Prepared in partial fulfillment of the requirements of GEOG 590: R for Earth System Science

1 Project Summary

1.1 R in Cartography

As a cartographer, I am often asked to break down complicated data stories into easily digestible information. Sometimes, this means I use visual communication to spatialize or depict qualitative stories. This is where the majority of my traditional cartography has focused. However, being an effective cartographer and spatial storyteller involves being competent at handling, cleaning, and distilling large datasets into digestible stories.

There are many toolboxes in R which allow for advanced cartographic visualizations. The Rayshader package, for example, is an opensource pachage that uses raytracing and hillshading algorithms to create stunning 2D and 3D maps. In the video linked below, Kati Perry, a cartographer at the Washington Post and UO Alumn, talks through how she used R to programmatically create a series of maps for the Associated Press showing Arctic ice loss over time.

NACIS video

1.2 Trade Project

My motivation behind this work is a project I have been developing for the last 3 years in collaboration with the InfoGraphics Lab and 2 economics professors– Dr. Woan Foong Wong and Dean of the Lundquist School of Business Bruce Blonigen. Drs. Wong and Blonigen have an NSF grant to research how the invention of the steam engine changed the relative importance of global cities to international trade. So, I have been working with Joanna Merson in the InfoGraphics Lab to create a heuristic model of global historic sailing routes in the early 1850s. Previous researchers have created similar models, but none were geographers, so they failed to accurately accomodate issues such as the earth’s projection or spatial limitations to travel.

1.3 This project

This project was motivated by my attempt to learn R for the first time and guided by my intention of getting more familiar with an update we made to our underlying dataset for the Trade project. The products, then, are less focused on how to perform earth system science using R, but rather trying new techniques for how to communicate earth system science using packages in R. I do not have this communication down quite yet, so these images are not intended to represent scientifically informative graphics, but rather a demonstration of my capabilities as a designer using this programming language.

The goal was to produce an R Markdown file which could be integrated integrated into my portfolio which highlights, not just my skills in R, but how I approach learning new applications. As such, it was important to me that the file not only be reflective of my learning progress, but that it also matches the overall aesthetic of my portfolio. Therefore, I focused on creating a custom cascading style sheet (CSS) file, using a custom font from Google fonts, and making sure the design and layout were easy to read and scan.

2 Methodology

In the sections below, I walk through how I came up with my three iterations on developing a data visualization of wind speed and direction in R.

2.1 Installing Libraries

Below are the libraries I used to install, manipulate, and plot my maps.

library(ncdf4) # For handling Net CDF files
library(sf) # Simple features in R
library(stars)
library(tmap) # Quickly drawing data, including simple datasets
library(rnaturalearth) # Quickly loading data from Natural Earth
library(metR) # Has a cool package for visualizing vectors
library(dplyr) 
library(classInt)
library(lattice)
library(ggplot2) # for plotting data
library(RColorBrewer) # for incorporating cartographically tried-and-true colors

2.2 Opening NetCDF File

The data I used was the ERA5 wind dataset from Copernicus. This is a dataset of atmospheric conditions at a very small temporal and spatial scale, which is incredibly useful for the purposes of our model.

# set path and filename
ncpath <- ".\\Data\\"
ncname <- "Wind_5Dec"  
ncfname <- paste(ncpath, ncname, ".nc", sep="")
dname <- "Wind_Dec"


# open a netCDF file
ncin <- nc_open(ncfname)
print(ncin)
## File .\Data\Wind_5Dec.nc (NC_FORMAT_64BIT):
## 
##      2 variables (excluding dimension variables):
##         short u10[longitude,latitude,time]   
##             scale_factor: 0.000725724893463742
##             add_offset: 1.09752226938921
##             _FillValue: -32767
##             missing_value: -32767
##             units: m s**-1
##             long_name: 10 metre U wind component
##         short v10[longitude,latitude,time]   
##             scale_factor: 0.000684122905143568
##             add_offset: -1.47976558440179
##             _FillValue: -32767
##             missing_value: -32767
##             units: m s**-1
##             long_name: 10 metre V wind component
## 
##      3 dimensions:
##         longitude  Size:1440 
##             units: degrees_east
##             long_name: longitude
##         latitude  Size:721 
##             units: degrees_north
##             long_name: latitude
##         time  Size:3 
##             units: hours since 1900-01-01 00:00:00.0
##             long_name: time
##             calendar: gregorian
## 
##     2 global attributes:
##         Conventions: CF-1.6
##         history: 2024-03-13 10:04:42 GMT by grib_to_netcdf-2.25.1: /opt/ecmwf/mars-client/bin/grib_to_netcdf.bin -S param -o /cache/data9/adaptor.mars.internal-1710324282.5448966-27416-6-d11edc9d-5e9e-4dc0-b32f-9dc14aac69fd.nc /cache/tmp/d11edc9d-5e9e-4dc0-b32f-9dc14aac69fd-adaptor.mars.internal-1710324273.6142778-27416-6-tmp.grib

2.3 Processing UV Data

In a vector, ‘U’ refers to the direction of a vector and ‘V’ refers to its magnitude.This section demonstrates how I pre-process UV data from the netCDF files.

# Read the longitude, latitude, and wind components
longitude <- ncvar_get(ncin, "longitude") -180
latitude <- ncvar_get(ncin, "latitude")
u10 <- ncvar_get(ncin, "u10")
v10 <- ncvar_get(ncin, "v10")


# Create a data frame
wind_data <- data.frame(
  lon = rep(longitude, each = length(latitude)),
  lat = rep(latitude, times = length(longitude)),
  u = as.vector(u10),
  v = as.vector(v10)
)

# Sample data to reduce rendering time
sample_size <- 1500  # Adjust this number to control the density of wind vectors
sample_indices <- sample(nrow(wind_data), sample_size)
sampled_wind_data <- wind_data[sample_indices, ]

2.4 Simple plot

This was my first attempt to plot out the wind data using GGplot. Here, I used the “arrow” geometry segment to represent wind vectors. This was mostly a proof-of-concept for knowing I had properly loaded the netCDF file and selected the U/V components.

# Plot wind vectors
ggplot(sampled_wind_data, aes(lon, lat, xend = lon + u, yend = lat + v)) +
  geom_segment(arrow = arrow(length=unit(.2, 'cm'))) +
  xlab("Longitude") +
  ylab("Latitude") +
  theme_bw() +
  coord_fixed() +
  ggtitle("Wind Vectors")

2.5 Adding in coastlines

Next, I really wanted to add coastlines to the back of my map to contextualize the scale. I was able to add coastlines using the tmap library. However, data protted through that library does not have a “rotate” function and I could not represent the wind vectors the way I wanted to.

One consistent challenge with this project was the speed of loading times for rendering the images. To ameliorate this for the class final, I chose to use a sample of the wind data to minimize the number of points the program had to render.

# Plot wind vectors and coastline together
data("World") #Load world data for tmap

wind_sf <- st_as_sf(wind_data, coords = c("lon", "lat"), crs = 4369)
## Warning in CPL_crs_from_input(x): GDAL Message 1: CRS EPSG:4369 is deprecated.
## Its non-deprecated replacement EPSG:4965 will be used instead. To use the
## original CRS, set the OSR_USE_NON_DEPRECATED configuration option to NO.
# Sample data to reduce rendering time
sample_size <- 1000  # Adjust this number to control the density of wind vectors
sample_indices <- sample(nrow(wind_data), sample_size)
sampled_wind_sf <- wind_sf[sample_indices, ]

tm_shape(World) +
   tm_fill() +
#   tm_borders() +
tm_shape(sampled_wind_sf) +
  tm_symbols(size = 0.1) +
  tm_layout(title = "Wind points with Coastline")

2.6 Swirling Seas

Since working with wind data, one of my biggest sources of cartographic inspiration has been NASA’s swirling seas data visualization, which shows the movement of winds across the oceans. I found a YouTuber named MilosMakesMaps who has a tutorial on how to recreate this visualization in R: https://www.youtube.com/watch?v=d8VZMjhBFwE&t=947s. His tutorial focused on Europe, so I adjusted his workflow to accommodate global winds.

## File .\Data\Wind_5Dec1940.nc (NC_FORMAT_64BIT):
## 
##      2 variables (excluding dimension variables):
##         short u100[longitude,latitude,time]   
##             scale_factor: 0.000804639933316039
##             add_offset: 4.97864736265053
##             _FillValue: -32767
##             missing_value: -32767
##             units: m s**-1
##             long_name: 100 metre U wind component
##         short v100[longitude,latitude,time]   
##             scale_factor: 0.000802896417644545
##             add_offset: -4.77635664840413
##             _FillValue: -32767
##             missing_value: -32767
##             units: m s**-1
##             long_name: 100 metre V wind component
## 
##      3 dimensions:
##         longitude  Size:1440 
##             units: degrees_east
##             long_name: longitude
##         latitude  Size:721 
##             units: degrees_north
##             long_name: latitude
##         time  Size:1 
##             units: hours since 1900-01-01 00:00:00.0
##             long_name: time
##             calendar: gregorian
## 
##     2 global attributes:
##         Conventions: CF-1.6
##         history: 2024-02-01 20:53:31 GMT by grib_to_netcdf-2.25.1: /opt/ecmwf/mars-client/bin/grib_to_netcdf.bin -S param -o /cache/data1/adaptor.mars.internal-1706820811.296947-18417-3-5305a79e-5fe6-4c2c-8ddc-784cc99c3265.nc /cache/tmp/5305a79e-5fe6-4c2c-8ddc-784cc99c3265-adaptor.mars.internal-1706820810.697953-18417-3-tmp.grib
coastline <- ne_coastline(scale = "small", returnclass = "sf")

ggplot(data = wind_data) +
  geom_sf(
        data = coastline,
        fill = "grey50",
        color = "grey",
        size = .25,
        alpha = .99
    ) +
  #coord_sf(crs = "+proj=robin") +
  metR::geom_streamline(
    aes(
      x = wind_data$lon,
      y = wind_data$lat,
      dx = wind_data$u,
      dy = wind_data$v,
      color = sqrt(after_stat(dx)^2 + after_stat(dy)^2),
      alpha = sqrt(after_stat(dx)^2 + after_stat(dy)^2),
      linewidth = sqrt(after_stat(dx)^2 + after_stat(dy)^2),
    ),
    L = 15,
    res = 5,
    n = 60,
    arrow = NULL,
    lineend = "round",
    inherit.aes = F
  ) +
    scale_color_gradientn(
      name = "Speed",
      colours = hcl.colors(
        12, "BluYl"
      ),
      breaks = waiver()
    ) + 
  scale_alpha(
    range = c(.1, .8)
  ) +
  scale_linewidth(
    range = c(.1, .5)
  ) +
  guides(
    alpha = "none",
    linewidth = "none",
    color = guide_colorbar(
      direction = "vertical",
      title.position = "left",
      label.position = "right",
      title.hjust = .5,
      label.hjust = 0,
      nrow = 1
    )
  ) +
  theme_void() +
  theme(
    legend.title = element_text(
      size = 11,
      color = "white"
    ),
    legend.text = element_text(
      size = 9,
      color = "white"
    ),
    plot.title = element_text(
      size = 17,
      color = "white",
      hjust = .1,
      vjust = -1
    ),
    plot.background = element_rect(
      fill = "black",
      color = NA
    ),
  ) +
  labs(
    title = "Winds in Dec 1940",
    subtitle = "Source: Climate change service, ERA5 hourly data December, 1940",
    x = "",
    y = ""
  )
## Warning: Use of `wind_data$lon` is discouraged.
## ℹ Use `lon` instead.
## Warning: Use of `wind_data$lat` is discouraged.
## ℹ Use `lat` instead.
## Warning: Use of `wind_data$u` is discouraged.
## ℹ Use `u` instead.
## Warning: Use of `wind_data$v` is discouraged.
## ℹ Use `v` instead.

3 Discussion and Next Steps

Overall, I am impressed by the capabilities of R to handle vast amounts of data and the R developer community for coming up with so many advanced tools for data visualization. Overall, I felt as though my ability to refine the visualizations was limited by my unfamiliarity with the netCDF format and the complex geometries of wind data as a whole. I am proud of what I was able to achieve in this visualization, but I do think I could use more Moving forward, I want to pick a dataset that I am more familiar with while I am learning the language and then I want to revisit and refine this project.

In the future, I want to further develop the aethetic placement of the title and legend. Eventually, I would like to recreate the “swirling seas” visualization using rayshader so that I can then create animations of the moving wind data. Overall, learning to develop these visualizations in R can help develop more reproducible workflows for long-term projects and is a useful step towards making my science more open.