Error calculating aspect from DEM in R using raster/terra

Question

I have a DEM file gathered using R. The raster and terra packages both have a function to calculate aspect, called terrain(). When using the terrain() function and the units="degrees" argument, the resulting raster should have a range of values from 0 to 360. However, sometimes when I run terrain, I get the values of 0 to 6.283185. Other times (on the same machine) I get a range of values from 0 to 360.

See the code below:

library(elevatr)
library(dplyr)
library(sf)
library(readxl)
library(raster)
library(terra)

tdir=tempdir()
countyurl = "https://www2.census.gov/geo/tiger/GENZ2018/shp/cb_2018_us_county_500k.zip"

if(file.exists(paste(tdir,"/cb_2018_us_county_500k.shp",sep=""))==F){
  download.file(countyurl,destfile = file.path(tdir,"Counties.zip"))
  unzip(file.path(tdir,"Counties.zip"),exdir=tdir)}
Herkimer = read_sf(paste(tdir,"/cb_2018_us_county_500k.shp",sep="")) %>%
  filter(STATEFP=="36") %>%
  filter(NAME == "Herkimer") %>%
  st_transform(.,crs=32618)

PRJ="+init=EPSG:4326"
DEM = get_elev_raster(locations = Herkimer,prj = PRJ,z = 10,neg_to_na = T) 
DEM = raster::mask(DEM,mask=Herkimer)
remove(PRJ)

terra_aspect = terra::terrain(DEM,opt='aspect',units='degrees')
terra_aspect

raster_aspect = raster::terrain(DEM,opt='aspect',units='degrees')
raster_aspect

With the resulting rasters:

terra_aspect
#> class      : RasterLayer 
#> dimensions : 3630, 1581, 5739030  (nrow, ncol, ncell)
#> resolution : 54.75867, 54.75867  (x, y)
#> extent     : 480758.6, 567332.1, 4711142, 4909916  (xmin, xmax, ymin, ymax)
#> crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
#> source     : memory
#> names      : aspect 
#> values     : 0, 6.283185  (min, max)

raster_aspect
#> class      : RasterLayer 
#> dimensions : 3630, 1581, 5739030  (nrow, ncol, ncell)
#> resolution : 54.75867, 54.75867  (x, y)
#> extent     : 480758.6, 567332.1, 4711142, 4909916  (xmin, xmax, ymin, ymax)
#> crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
#> source     : memory
#> names      : aspect 
#> values     : 0, 6.283185  (min, max)

Answer 1

You have a typo in your terrain() function calls so the output is in the default unit (radians). Just change the word 'units' to 'unit' so that:

terrain(DEM, opt='aspect', unit='degrees')

I also made a couple of minor changes so that your CRS is consistently WGS84/Pseudomercator and using the PROJ6 syntax, so here's the full script:

library(elevatr)
library(dplyr)
library(sf)
library(readxl)
library(raster)
library(terra)
library(sp)

tdir=tempdir()
countyurl = "https://www2.census.gov/geo/tiger/GENZ2018/shp/cb_2018_us_county_500k.zip"

if(file.exists(paste(tdir,"/cb_2018_us_county_500k.shp",sep=""))==F){
  download.file(countyurl,destfile = file.path(tdir,"Counties.zip"))
  unzip(file.path(tdir,"Counties.zip"),exdir=tdir)}
Herkimer = read_sf(paste(tdir,"/cb_2018_us_county_500k.shp",sep="")) %>%
  filter(STATEFP=="36") %>%
  filter(NAME == "Herkimer") %>%
  st_transform(.,crs=3857)

PRJ="EPSG:3857"
DEM = get_elev_raster(locations = Herkimer ,prj = PRJ, z = 10, neg_to_na = T) 
DEM = mask(DEM,mask=Herkimer)
remove(PRJ)

terra_aspect = terra::terrain(DEM,opt='aspect',unit='degrees')
terra_aspect

raster_aspect = raster::terrain(DEM,opt='aspect',unit='degrees')
raster_aspect

Results: