This tutorial was originally developed as a lab activity in a graduate-level ecologically niche modelling course I taught at the University of Florida with my former PI, Rob Guralnick. The original gist can be found here: https://gist.github.com/hannahlowens/931fe265598a71b47ac0f597e3351d0d. I have updated it to use terra instead of raster, sp, and rgdal, and incorporated some lessons I have learned along the way to introduce spacial analysis in R and writing clearer code.

Setup

First thing’s first–we need to load the terra package and set the working directory.

# October 2023 Basic GIS Operations Before Niche Modeling
# R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
# Copyright (C) 2023 The R Foundation for Statistical Computing
# Platform: aarch64-apple-darwin20 (64-bit)

#Load the libraries you'll need.
library(terra) 
## terra 1.7.55
library(tidyterra)
## 
## Attaching package: 'tidyterra'
## The following object is masked from 'package:stats':
## 
##     filter
library(ggplot2)

# Set the working directory to the location of the data
setwd("~/Dropbox/ENMSeminar/Labs:Homeworks/Lab2/Lab2Data/")

Loading Data

Basics

Here are examples of how to read in vector data from .csv and .shp files, and raster data from .tif files.

#Read points and turn into shapefile
shermanSquirrels <- read.csv(file = "ShermanSquirrels.csv", header = T)
head(shermanSquirrels)
##         Species Longitude Latitude
## 1 Sciurus niger -83.36905 31.40090
## 2 Sciurus niger -82.60067 28.46997
## 3 Sciurus niger -81.46268 28.71139
## 4 Sciurus niger -81.18965 28.38080
## 5 Sciurus niger -81.18860 28.38107
## 6 Sciurus niger -81.18832 28.38104
shermanSquirrelsShp <- vect(shermanSquirrels, geom = c("Longitude", "Latitude"))
shermanSquirrelsShp
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 35, 1  (geometries, attributes)
##  extent      : -87.04154, -80.0361, 26.7061, 34.82  (xmin, xmax, ymin, ymax)
##  coord. ref. :  
##  names       :       Species
##  type        :         <chr>
##  values      : Sciurus niger
##                Sciurus niger
##                Sciurus niger
# Reading shapefiles
florida <- vect("FloridaBound1/FLORIDA.shp")
florida
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 1, 4  (geometries, attributes)
##  extent      : -1110383, 1298633, 77020.92, 2466420  (xmin, xmax, ymin, ymax)
##  source      : FLORIDA.shp
##  coord. ref. : NAD83(HARN) / Florida West (ftUS) 
##  names       : OBJECTID  DATESTAMP SHAPE_AREA SHAPE_LEN
##  type        :    <int>      <chr>      <num>     <num>
##  values      :        1 2000/05/16   1.57e+12 2.834e+07
# Read rasters
altitude <- rast("ETOPO_2022_v1_60s_N90W180_bed.tif")
altitude
## class       : SpatRaster 
## dimensions  : 10800, 21600, 1  (nrow, ncol, nlyr)
## resolution  : 0.01666667, 0.01666667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : ETOPO_2022_v1_60s_N90W180_bed.tif 
## name        : ETOPO_2022_v1_60s_N90W180_bed

Now let’s try plotting it.

# Let's try plotting
plot(altitude, main="Altitude", xlim = c(-88.5,-79), ylim = c(24, 32))
plot(shermanSquirrelsShp, pch = 16, cex = 1.5, add = T)
plot(florida, add = T)

Where did Florida go?! The problem becomes clear if you use crs() to look at the coordinate reference systems for each file.

# No coordinate reference system
crs(shermanSquirrelsShp)
## [1] ""
# Altitude is WGS 84
crs(altitude)
## [1] "COMPOUNDCRS[\"WGS 84 + EGM2008 height\",\n    GEOGCRS[\"WGS 84\",\n        DATUM[\"World Geodetic System 1984\",\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        CS[ellipsoidal,2],\n            AXIS[\"geodetic latitude (Lat)\",north,\n                ORDER[1],\n                ANGLEUNIT[\"degree\",0.0174532925199433]],\n            AXIS[\"geodetic longitude (Lon)\",east,\n                ORDER[2],\n                ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    VERTCRS[\"EGM2008 height\",\n        VDATUM[\"EGM2008 geoid\"],\n        CS[vertical,1],\n            AXIS[\"gravity-related height\",up,\n                LENGTHUNIT[\"metre\",1]],\n        ID[\"EPSG\",3855]]]"
# Florida polygon is NAD83
crs(florida)
## [1] "PROJCRS[\"NAD83(HARN) / Florida West (ftUS)\",\n    BASEGEOGCRS[\"NAD83(HARN)\",\n        DATUM[\"NAD83 (High Accuracy Reference Network)\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]],\n            ID[\"EPSG\",6152]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"Degree\",0.0174532925199433]]],\n    CONVERSION[\"unnamed\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",24.3333333333333,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-82,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.999941176470588,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",656166.666666667,\n            LENGTHUNIT[\"US survey foot\",0.304800609601219],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"US survey foot\",0.304800609601219],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"US survey foot\",0.304800609601219,\n                ID[\"EPSG\",9003]]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"US survey foot\",0.304800609601219,\n                ID[\"EPSG\",9003]]]]"

So we need to fix these inconsistencies.

Defining coordinate reference systems and re-projecting data

This is not the place to explain coordinate refernce systems and projections. Here’s a helpful explainer: https://pro.arcgis.com/en/pro-app/latest/help/mapping/properties/coordinate-systems-and-projections.htm

Let’s fix these files.

# Defining coordinate reference system
crs(shermanSquirrelsShp) <- crs(altitude)

# Reprojecting the shapefiles
shermanSquirrelsShp <- project(shermanSquirrelsShp, crs(florida))
altitude_prj <- project(altitude, crs(florida))
## 
|---------|---------|---------|---------|
=========================================

Now let’s try plotting it again. I used ext() to limit the extent of the altitude plot to that of the Florida shapefile. The buffer argument allows a little breathing room around the Florida shapefile, so we can see some of our squirrel observations that are outside of the state.

# Let's try plotting
plot(altitude_prj, main="Altitude", 
     xlim = ext(florida)[1:2], 
     ylim = ext(florida)[3:4], buffer = T)
plot(shermanSquirrelsShp, pch = 16, cex = 1.5, add = T)
plot(florida, add = T)

Great! Florida’s back! Now let’s check out some basic GIS operations that niche modelers often use.

An introduction to GIS operations in terra

Cropping and masking

Masking and cropping are almost identical to how raster set things up. Here’s an example. Raster cropping changes the extent of a raster to match the extent of the shapefile. Thus, it’s more computationally efficient to crop first, then mask if you want to limit raster data to a specific shape.

#Cropping rasters to extent of shapefile
ext(florida)
## SpatExtent : -1110382.95, 1298632.755, 77020.924, 2466419.852 (xmin, xmax, ymin, ymax)
ext(altitude_prj)
## SpatExtent : -54144746.0457589, 55451085.6168872, -97333092.2915943, 79668529.33174 (xmin, xmax, ymin, ymax)
altitudeCropped <- crop(altitude_prj, florida)
plot(altitudeCropped)

Remember that for rasters, masking changes cells outside a raster to another value (NA by default), although terra can also do inverse masking, which is great for those of us who work in marine systems.

#Masking rasters to shapefiles
altitudeCM <- mask(altitudeCropped, mask = florida)
plot(altitudeCM, main="Altitude", buffer = TRUE)
plot(florida, add = TRUE)
plot(shermanSquirrelsShp, add = TRUE)

Nice and tidy! Now let’s save that new raster file. There are no surprises here when comparing raster to terra.

#Save result as a .tif file
writeRaster(altitudeCM, filename = "AltitudeCroppedFlorida.tif", overwrite = TRUE)

Looking at the map, I see there are points outside our shapefile. How do you remove those?

Spatial analysis

First, we will identify which points are contained within the shapefile. The is.related() function gives you a vector of logical states (i.e. true or false) for the relationship you are interested in. In this case, for each point in the squirrel shapefile, we are asking if it falls within the florida shapefile.

floridaSquirrels <- is.related(shermanSquirrelsShp, florida, relation = "within")
plot(shermanSquirrelsShp, main = "Sherman Squirrels of Florida", 
     col = ifelse(floridaSquirrels, "red", "black"))
plot(florida, add = TRUE)

Ok, there are three points to remove using logical operations the same way you would any vector-based object.

flShermanSquirrelsShp <- shermanSquirrelsShp[floridaSquirrels]

It is also very straightforward to extract raster data to your point shapefile. As an example, let’s get the altitude at each observation point.

flShermanSquirrelsShp$Elevation <- extract(x = altitudeCM, y = flShermanSquirrelsShp, ID = FALSE)
head(flShermanSquirrelsShp)
##         Species Elevation
## 1 Sciurus niger  5.959705
## 2 Sciurus niger 12.509106
## 3 Sciurus niger 18.980927
## 4 Sciurus niger 18.980927
## 5 Sciurus niger 18.980927
## 6 Sciurus niger 48.024765

Now let’s create a new data.frame from the squirrel shapefile and save it as a .csv for futher use later.

flShermanSquirrels_Coords <- geom(flShermanSquirrelsShp)[,c("x", "y")] # Gets coordinates
flShermanSquirrels_Data <- data.frame(flShermanSquirrelsShp) # Gets shapefile data
flShermanSquirrels_Table <- cbind(flShermanSquirrels_Coords, flShermanSquirrels_Data)
colnames(flShermanSquirrels_Table)[1:2] <- c("Longitude", "Latitude")

head(flShermanSquirrels_Table)
##   Longitude Latitude       Species Elevation
## 1  463183.1  1504052 Sciurus niger  5.959705
## 2  828404.8  1591736 Sciurus niger 12.509106
## 3  916739.2  1472025 Sciurus niger 18.980927
## 4  917076.2  1472125 Sciurus niger 18.980927
## 5  917166.3  1472115 Sciurus niger 18.980927
## 6  672360.8  1347055 Sciurus niger 48.024765
write.csv(flShermanSquirrels_Table, "shermanSquirrelsFlorida.csv", row.names = FALSE)

More sophisticated mapping

For those of you who would like to make more sophisticated maps and are fans of the tidyverse, let’s make a fancier (but still pretty basic) map using the ggplot2 and tidyterra packages. Again, it is not the goal of this document to provide a full tutorial, but this will give you a taste of what’s possible.

# Reproject data to WGS 84
altitudeCM_upj <- project(altitudeCM, crs(altitude))
altitudeTidyMap <- trim(altitudeCM_upj)
florida_upj <- project(florida, altitude)
squirrel_upj <- project(flShermanSquirrelsShp, altitude)

# Define a "clear" color for no data values
my_col <- rgb(0, 0, 255, max = 255, alpha = 0, names = "clear")

squirrelMap <- ggplot() + 
  geom_spatraster(data = altitudeTidyMap, aes(fill = ETOPO_2022_v1_60s_N90W180_bed)) +
  scale_fill_viridis_c(name = "Altitude", option = "viridis", alpha = .7, na.value = my_col) +
  theme_minimal() + coord_sf(expand = TRUE, label_graticule = "SW") +
  ggtitle("Occurrences of Sherman's Fox Squirrel in Florida") +
  geom_sf(data = squirrel_upj) +
  geom_sf(data = florida_upj, bg = my_col, col = "black")
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
squirrelMap

And let’s save the map as a .png.

#Saving a map as a .png
png(filename="shermanSquirrelsFloridaMapUsingR.png")
squirrelMap
dev.off()
## quartz_off_screen 
##                 2

References

Hernangomez D (2023). tidyterra: tidyverse Methods and ggplot2 Helpers for terra Objects. R package version 0.4.0, https://doi.org/10.5281/zenodo.6572471, https://dieghernan.github.io/tidyterra/

Hijmans R (2023). terra: Spatial Data Analysis. R package version 1.7-55, https://CRAN.R-project.org/package=terra.

Wickham H (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York.