Map Axes and Kernel Density Estimates (KDEs)

Abstract

The purpose of this document is to instruct you on how to add kernel density estimates to your plot axes. In this case, axes will be applied to modeled geographic data.

1 Getting started

First you’ll need to ensure that you’ve loaded the appropriate packages for the task. As I will be manipulating light pollution data, my ‘geographic’ packages are as follows:

if(!require(pacman)) install.packages("pacman")

Loading required package: pacman

p_load(raster, terra, sf)

Tip

For brevity, I have used the p_load(...) function from the pacman package to load:

• raster
• terra
• sf

Axes, are best added using the ggplot2 and its extensions, ggspatial and ggExtra.

p_load(ggplot2, ggspatial, ggExtra)

2 Preparing your data

First, I will load my data. Here, I have saved an old workspace that contains a shapefile loaded as ‘Bristol_data’. This file contains a field that relates to light pollution values in Bristol.

You can simply replace this step with a process you prefer to use to load your shapefile data.

load(url("https://github.com/Richtea84/StudentHelp/raw/refs/heads/main/Bristol_data.RData"))

2.1 Optional step:

I want to create a choropleth that features distributions for “High” and “Low” light pollution values in my data.

Bristol_data$Hi_Low <- ifelse(Bristol_data$LP < median(Bristol_data$LP), "Low", "High")

We can preview the result using ggplot2 (Figure 1).

ggplot(data = Bristol_data)+
  geom_sf(aes(color = Hi_Low))+
  ggspatial::annotation_north_arrow()+
  theme_minimal()

Figure 1: Previewing the data that will form the KDE axes plots. Note: point geometries must be active - polygons can be added later.

3 Adding KDEs to the axes

The first step in the process is to save the initial ggplot as an object. Here, it is called map_obj:

ggplot(data = Bristol_data)+
  geom_point(aes(x = lon, y = lat, color = Hi_Low))+
  theme_minimal() -> map_obj

We now use the ggMarginal(...) function to on the object and specify basic colour and fill aesthetics (Figure 2).

ggMarginal(map_obj, groupColour = T, groupFill = T)

Figure 2: The point geometries with accompanying histograms. Here, a generic plot has been generated that hasn’t been snapped to a geographic projection yet.

3.1 Incorporating geographic projections

To project your data as a map, you need to invoke the geom_sf(...) function in your ggplot object. In this example, we will attempt to show polygon geometries instead of the points - the data we’ll be using is called geog_data.

To prepare this data, first we ensure the geometries are set to longitude and latitude (WGS: 4326). A data join can then be accomplished (object called geog_data_2:

# Syncing the projection:
geog_data <- st_transform(geog_data, 4326)

# Joining the polygonal data to the points geometries
geog_data_2 <- merge(x = geog_data, as.data.frame(st_drop_geometry(Bristol_data)), by = "lso11cd")

Tip

The merge(...) function doesn’t work for two sf objects. Therefore, you need to convert points data into a data frame that excludes geometries.

3.2 The map

Finally, we can produce the map (Figure 3).

ggplot(data = Bristol_data)+
  geom_sf(data = geog_data)+
  geom_point(aes(x = lon, y = lat, color = Hi_Low), alpha = 0,show.legend = FALSE)+
  labs(fill = "Light Pollution")+
  theme_minimal()+
  theme(legend.position = "bottom") -> map_obj

map_obj+
  geom_sf(data = geog_data_2, aes(fill = Hi_Low), alpha = 0.7)+
  ggspatial::annotation_north_arrow()-> map_obj2

ggMarginal(map_obj2, groupColour = T, groupFill = T)

Figure 3: The final rendered map featuring polygon geometeries. Note that the KDEs are still driven by the centroid geometries.