converting sf to raster using predefined cell size and plot it

Question

nz_height and nz are present in spData package.

1. nz_height - sf class dataset - Top 101 highest points in New Zealand.
2. nz - sf class data - Polygons representing the 16 regions of New Zealand.

Below is the plot for visualization using tmap()

Can you please help me create

1. nz_raster (raster class) with cell size = 100km by 100km and each cell contains number of peaks in it
2. plot nz_raster.

# load packages
library(tidyverse)
library(sf)
library(raster)
library(spData)
library(tmap)


# vector plot of peaks
tm_shape(nz) + 
  tm_polygons(col = "white") + 
  tm_shape(nz_height) + 
  tm_symbols(col = "red") + 
  tm_scale_bar()

Answer 1

You could rasterize the points with the raster package or with its successor terra:

# Load packages
packs <- list("tidyverse", "raster", "sf", "terra")
lapply(packs, require, character.only = T)

# raster version
nz_height <- st_transform(nz_height, crs = "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs") %>% 
  as("Spatial")
nz_raster <- raster(resolution = 100000, crs = "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs", 
  xmn = 13600000, xmx = 15100000, ymn = -5300000, ymx = -4700000) %>% 
  rasterize(nz_height, ., field = "elevation", fun = "count", background = 0)

# terra version
nz_height <- st_transform(nz_height, crs = "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs") %>% 
  vect()
nz_raster <- rast(resolution = 100000, crs = "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs", 
  xmin = 13600000, xmax = 15100000, ymin = -5300000, ymax = -4700000) %>% 
  rasterize(nz_height, ., fun = length, background = 0)

The code uses a Mollweide equal area projection, i.e. the grid cells are equally large.

You can plot the raster object via the plot() function that is part of these packages. Other options are e.g. rasterVis::gplot(), ggplot2::geom_raster(), and tmap::tm_raster().