Creating Multiple Density Maps with the Same Extent

Introduction

In this article, we will explore how to create multiple density maps from points using the smooth_map function from the tmaptools package. The goal is to have all rasters have the same extent, given by a shapefile. We will cover the necessary steps, including data preparation, reprojection, and resampling.

Prerequisites

Before starting, ensure you have the required packages installed:

tmaptools
rgdal
sf
raster

You can install these packages using R’s package manager:

install.packages("tmaptools")
install.packages("rgdal")
install.packages("sf")
install.packages("raster")

Data Preparation

Let’s start with a reproducible example of a dataframe (df) containing points with their respective longitude, latitude, and var:

library(tmaptools)
library(rgdal)
library(sf)
library(raster)

# Create dataframe
df <- setNames(data.frame(matrix(ncol = 3, nrow = 7)), c("longitude", "latitude", "var"))

df$longitude <- c(-53.30002, -50.47749, -55.85561, -57.88447, -55.83864, -58.84610, -57.49215)
df$latitude <- c(-13.037530,  -13.023480, -13.416200, -13.659120, -12.758670, -13.114460, -14.622520)
df$var <- c("A", "A", "A", "A", "B", "B", "B")

# Convert var to factor
df$var <- as.factor(df$var)

# Read shapefile of Mato Grosso, Brazil
study_area <- readOGR("shape_MT.shp", layer="shape_MT") # Will load the shapefile 

create empty stack
s <- list()

# Loop over factor levels
for (i in levels(df$var)) { 
  a <- subset(df, df$var == i)
  my.sf.point <- st_as_sf(x = a, coords = c("longitude", "latitude"), crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
  density_map <- smooth_map(my.sf.point, cover = study_area) 
  r <- density_map$raster
  
  # Print extent of the raster
  print(extent(r))
  
  s[[i]] <- r
}

s[["A"]] <- resample(s[["A"]], s[["B"]], method='bilinear')

stack(s)

Reprojection and Resampling

To achieve the desired extent for all density maps, we need to reproject them to a common coordinate reference system (CRS). In this example, we’ll use the ESRI:102033 South_America_Albers_Equal_Area_Conic CRS.

# Define CRS
crs.laea <- CRS("+proj=aea +lat_1=-5 +lat_2=-42 +lat_0=-32 +lon_0=-60 +x_0=0 +y_0=0 +ellps=aust_SA +units=m +no_defs")

# Reproject my.sf.point to ESRI:102033
my.sf.point <- st_transform(my.sf.point, 102033)

Next, we’ll resample the rasters to match the extent of the study area:

# Define extent of the study area
extent(r) <- c(-175916.2, 1048374, 1610251, 2845008) # extent of the study area

# Plot raster
plot(r)

# Resample rasters to match study area extent
s[["A"]] <- resample(s[["A"]], s[["B"]], method='bilinear')

stack(s)

Conclusion

Creating multiple density maps with the same extent involves reprojecting and resampling them. By following these steps, you can ensure that your rasters have a consistent spatial reference system, making it easier to compare and analyze them.

In this article, we explored how to create density maps from points using smooth_map and reprojection techniques to achieve common CRS for all maps. The code provided is designed to be easy to follow and can be adapted to fit your specific requirements.

Last modified on 2023-11-21