Understanding Sf and Geospatial Mapping in R for Accurate Arctic Maps with Circular Masks

Understanding Sf and Geospatial Mapping in R

=====================================================

As a technical blogger, it’s essential to delve into the world of sf, a powerful geospatial package for R. In this article, we’ll explore the basics of sf and apply its capabilities to create an Arctic map with a circular mask.

Introduction to Sf

sf (Simple Features) is a lightweight package that provides a flexible and efficient way to work with geometric data in R. It’s built on top of the Geospatial Analysis Library (Geographe), which offers robust support for geospatial operations, including spatial joins, buffering, and projections.

Installing Sf

Before we begin our journey into sf, you’ll need to install it using the package manager of your choice. For most users, this can be done via the CRAN website.

# Install sf
install.packages("sf")

Geospatial Data in R

Geospatial data is a type of spatial data that represents geographic locations and features. In R, we often work with these datasets using packages like sf and ggplot2.

Reprojection of Sf Objects

Reprojection is the process of transforming a geospatial dataset from one coordinate system to another. This is essential when working with different projection systems or when combining data from various sources.

Creating an Arctic Map with a Circular Mask

To create our arctic map, we’ll start by importing necessary libraries and loading the rnaturalearth::ne_coastline package for demonstration purposes.

library(sf)
library(tidyverse)

# Load example coastline dataset for illustration
coastlines <- rnaturalearth::ne_coastline(scale = 50, returnclass = "sf") 
st_crs(coastlines) <- "EPSG:4326"

Next, we’ll define our arctic bounding box using data frames. This will serve as the base for our circular mask.

# Create Arctic and non-Arctic bounding boxes
arcticBox <- data.frame(id = "A",
                        lon = c(-180, 180, 180, -180, -180),
                        lat = c(90, 90, 50, 50, 90))

otherBox <- data.frame(id = "B",
                       lon = c(-180, 180, 180, -180, -180),
                       lat = c(50, 50, -90, -90, 50))

boxes <- bind_rows(arcticBox, otherBox)

We’ll then convert this data into sf objects and combine them using st_combine.

# Convert bounding boxes to sf objects
boxesSF <- st_as_sf(boxes, coords = c("lon", "lat"), crs = "EPSG:4326") %>% 
  group_by(id) %>% summarise(geometry = st_combine(geometry)) %>% 
  st_cast("POLYGON")

# Display combined bounding box data
print(boxesSF)

With our bounding boxes in place, we’ll now create the Arctic map. We’ll use ggplot2 to visualize this data.

# Create Arctic map using ggplot2
p1 <- ggplot() + 
  geom_sf(data = coastlines, color = "blue", size = 0.25) +
  geom_sf(data = boxesSF, mapping = aes(fill = id), color = NA) +
  scale_fill_manual(values = alpha(colour = c("blue", "grey"), 
                                 alpha = c(0,1))) +
  coord_sf(ylim = c(0,90)) +
  theme(legend.position = "none")

Reprojecting the Map

To reproject our map to the North Pole Azimuthal Equidistant projection, we’ll add a new layer to our ggplot2 plot.

# Add North Pole Azimuthal Equidistant projection
p1 + coord_sf(crs = sf::st_crs("ESRI:102016"))

Removing Masked Areas

However, something strange is happening when reprojecting. The mask appears to be missing!

To resolve this issue, we’ll need to revisit our bounding boxes and ensure that they’re correctly defined.

# Update arcticBox for correct geometry
arcticBox <- data.frame(id = "A",
                        lon = c(-180, 180),
                        lat = c(90, -90))

boxesSF <- st_as_sf(arcticBox) %>% 
  group_by(id) %>% summarise(geometry = st_combine(geometry)) %>% 
  st_cast("POLYGON")

By making these adjustments to our bounding boxes and reprojecting the map, we should now see a circular mask in place.

Conclusion

Reproducing an arctic map with a circular mask requires careful consideration of reprojection, spatial joins, and masking polygons. By utilizing sf and ggplot2, we can effectively create such maps while maintaining accurate geographic data representation.

In conclusion, when working with geospatial data in R, it’s essential to understand the capabilities of packages like sf and how they’re used in conjunction with other libraries, including ggplot2.

Last modified on 2024-05-02