1

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)
3
  • You are using two different CRS with different units, EPSG:4326 (degrees) and EPSG: 32618 (meters). Try changing the definition of PRJ to PRJ="+init=EPSG:32618" and see if that solves your problem.
    – Kartograaf
    Mar 23 at 21:18
  • 1
    Good suggestion, but still no avail. The code is a reprex, were you able to recreate the error? Mar 23 at 21:22
  • Got it working, looks like it was just a typo!
    – Kartograaf
    Mar 23 at 22:03

1 Answer 1

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:

enter image description here

3
  • Amazing. I guess the function still runs for some reason with what should be an unknown argument? Mar 23 at 23:31
  • It appears that this is precisely the case. I reported the bug here so the powers that be can patch it up. github.com/rspatial/raster/issues/259
    – Kartograaf
    Mar 23 at 23:37
  • I appreciate that! I’m glad to know it wasn’t totally my stupidity. I can’t remember the last time I posted to Stack Exchange with the error being due to a simple typo. Mar 23 at 23:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.