This tutorial will guide you in a step-by-step process to create create a distance raster from spatial points. This tutorial considers that you are already familiar with the following concepts. If you are not, please visit the links.

Introduction

Urban studies commonly depend on assessing how city services are distributed in the city. However, technical teams commonly struggle to generate, maintain or retrieve detailed spatial information. The following R code explains how to create a raster (GeoTif) in which each pixel includes the minimum distance (crow-flight) to the closest hospital. In this tutorial we will use data from hospitals in Mexico City and the output from the Mask and crop raster from shapefile 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 sp: Classes and Methods for Spatial Data, Version 1.3-1, June 5, 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('sp')        # Run only once
# install.packages('leaflet')   # Run only once
library(raster)
library(rgdal)
library(sp)
library(leaflet)

2. Prepare data

If you want to replicate this tutorial, please download the data from our GitHub repository: distance_raster. 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 CSV file, with the location of all public hospitals in Denpasar.

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

built_up  <- raster("_input/denpasar_builtup.tif")
plot(built_up)

hospitals <- read.csv("_input/hospitals.csv")
plot(hospitals[c("lng", "lat")])

Please note from both plots that hospitals and built-up have different CRS.

3. Estimate the distance

The distance estimation requires only a couple of lines. However, you might need to transform your data to match CRS between layers. In our example we will transform all hospital points to the built-up raster’s CRS.

ghsl_crs <- '+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs'

# Read the hospitals data 
points <- SpatialPointsDataFrame(coords = hospitals[c("lng", "lat")], data = hospitals,
                                 proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))

# Transform the CRS
points <- spTransform(x = points, CRSobj = ghsl_crs)

Once that we have the same CRS in both sources, we can estimate the distances with the following code. Please note that the distanceFromPoints function will estimate the shortest distance (in meters) to the nearest hospital.

# Distance from points to each raster cell
distances <- distanceFromPoints(object = built_up, xy = points)
plot(distances)

# Lastly, we will "mask" the raster to keep only those cells inside CDMX
distance_raster <- mask(x = distances, mask = built_up)
plot(distance_raster)

4. Plot

Lastly, we will plot the final result using leaflet.

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

leaflet() %>%
  addProviderTiles(providers$OpenStreetMap.BlackAndWhite) %>%
  addRasterImage(distance_raster, colors = pal, opacity = 0.5) %>%
  addLegend(pal = pal, values = values(distance_raster),
    title = "Distance in meters")

5. Save your raster file

Save your raster file in your preferred directory.

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

That’s it!