This tutorial will guide you in a step-by-step process to mask and crop a raster from shapefile in R.

Introduction

Spatial analysis strongly depends defining a study area. Mask rasters, i.e. rasters with boolean values, simplify operations for retrieving a specific area. In this tutorial we will create mask and crop a raster from a political boundary shapefile. This will be used to redefine our “study area” (output from the Create a raster with built-up area from GHSL data tutorial).

1. Load required libraries

In this tutorial we will use the following libraries: raster: Title Geographic Data Analysis and Modeling Version 2.6-7, November 12, 2017 rgdal: Bindings for the ‘Geospatial’ Data Abstraction, Version 1.3-4, August 3, 2018 leaflet: Create Interactive Web Maps with the JavaScript ‘Leaflet’, Version 2.0.2, August 27, 2018

Run the following lines to load the libraries.

# install.packages('raster')    # Run only once
# install.packages('rgdal')     # Run only once
# install.packages('leaflet')   # Run only once
library(raster)
library(rgdal)
library(leaflet)

2. Prepare data

If you want to replicate this tutorial, you can download the data from our GitHub repository: Mask and crop a raster from_shapefile. Notice that you will need the contents of the _input file. As you will notice inputs include a TIF file (raster) with the built-up are in Denpasar Metro Area and a SHP file (shapefile) with the administrative boundaries of Denpasar.

First of all, we will import and explore the data.

built_up  <- raster('_input/built_up_denpasar.tif')
denpasar  <- readOGR('_input/denpasar_metro_area/denpasar_metro_area.shp')
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/ricardo/GitHub/ST_preprocessing/mask_raster_from_shp/_input/denpasar_metro_area/denpasar_metro_area.shp", layer: "denpasar_metro_area"
## with 128 features
## It has 11 fields
## Warning in readOGR("_input/denpasar_metro_area/denpasar_metro_area.shp"):
## Z-dimension discarded
# If needed, transform the polygon to match CRS between layers. 
denpasar <- spTransform(x = denpasar, CRSobj = '+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')
plot(denpasar)

plot(built_up)

3. Mask the raster

Once that the raster and the shapefile are in the same CRS, creating a mask is quite straightforward.

masked <- mask(x = built_up, mask = denpasar)
plot(masked)

4. Crop the raster

As you can see from the figure below, the raster now has plenty of white space (cells with NA values). If that is your case, you might need to crop the raster. In this tutorial, we will use the shapefile extent to crop the raster.

cropped <- crop(x = masked, y = extent(denpasar))
plot(cropped)

5. Plot

Lastly, we will plot the final result using leaflet.

pal <- colorNumeric(c("#556270", "#4ECDC4", "#C7F464", "#FF6B6B", "#C44D58"),
                    values(cropped),
                    na.color = "transparent")

leaflet() %>%
  addProviderTiles(providers$OpenStreetMap.BlackAndWhite) %>%
  addRasterImage(cropped, colors = pal, opacity = 0.5) %>%
  addLegend(pal = pal, values = values(cropped),
    title = "Built-up area")

6. Save your raster file

Save your raster file in your preferred directory.

writeRaster(cropped, filename = "_output/denpasar_builtup.tif", overwrite = TRUE)

That’s it!