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)