Terrain() function for computing slope and aspect from elevation data always returns NA

I'm facing a problem with the terrain() function in the raster package (in R) since yesterday. I'm working with elevation data and coordinates in Lambert 72 (espg registry : "+init=epsg:31300"). Here is the script I'm using, first to create the raster, and then to calculate (or trying to do so) the slope and aspect from it:

``````data = read.table("Coordinates_piquets.txt", h = T, sep = "\t")
unitname(data[,3]) <- c("metre", "metres")
X = data[,1]
Y = data[,2]
Z = data[,3]
``````

The `"coordinates_piquets.txt"` file is made up of 3 rows and 184 lines. First row: latitude; second row: longitude; third row: elevation.

``````library(raster)
data = matrix(c(X,Y,Z),  ncol=3,  byrow=FALSE)
e = extent(data[,1:2])
r=raster(e, ncol=3, nrow=184, crs = CRS("+init=epsg:31300"))
print(r)

class       : RasterLayer
dimensions  : 184, 3, 552  (nrow, ncol, ncell)
resolution  : 65.93271, 2.17544  (x, y)
extent      : 86799.25, 86997.04, 171858, 172258.3  (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:31300 +proj=lcc +lat_1=49.83333333333334 +lat_2=51.16666666666666 +lat_0=90 +lon_0=4.356939722222222 +x_0=150000.01256 +y_0=5400088.4378 +ellps=intl +towgs84=-106.8686,52.2978,-103.7239,0.3366,-0.457,1.8422,-1.2747 +units=m +no_defs

x = rasterize(data[,1:2], r, data[,3], fun=mean)
summary(x)
layer
Min.    160.3210
1st Qu. 174.3855
Median  179.6935
3rd Qu. 185.6319
Max.    199.7430
NA's    400.0000

slope_asp = terrain(x, opt=c('slope', 'aspect'), unit='degrees', neighbors=8)
plot(slope_asp)
summary(slope_asp)

[,1] [,2]
Min.      NA   NA
1st Qu.   NA   NA
Median    NA   NA
3rd Qu.   NA   NA
Max.      NA   NA
NA's     552  552
``````

Does anyone have an idea of why I'm getting NA values? I don't get any error or warning message and the plots I'm trying to draw from the terrain() function results are empty. I'm very new to GIS so the solution could be straightforward.

str(slope_asp) returns:

``````Formal class 'RasterBrick' [package "raster"] with 12 slots
..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
.. .. ..@ name        : chr ""
.. .. ..@ datanotation: chr "FLT4S"
.. .. ..@ byteorder   : chr "little"
.. .. ..@ nodatavalue : num -Inf
.. .. ..@ NAchanged   : logi FALSE
.. .. ..@ nbands      : int 1
.. .. ..@ bandorder   : chr "BIL"
.. .. ..@ offset      : int 0
.. .. ..@ toptobottom : logi TRUE
.. .. ..@ blockrows   : int 0
.. .. ..@ blockcols   : int 0
.. .. ..@ driver      : chr ""
.. .. ..@ open        : logi FALSE
..@ data    :Formal class '.MultipleRasterData' [package "raster"] with 14 slots
.. .. ..@ values    : num [1:552, 1:2] NA NA NA NA NaN NA NA NaN NA NA ...
.. .. ..@ offset    : num 0
.. .. ..@ gain      : num 1
.. .. ..@ inmemory  : logi TRUE
.. .. ..@ fromdisk  : logi FALSE
.. .. ..@ nlayers   : int 2
.. .. ..@ dropped   : NULL
.. .. ..@ isfactor  : logi [1:2] FALSE FALSE
.. .. ..@ attributes: list()
.. .. ..@ haveminmax: logi TRUE
.. .. ..@ min       : num [1:2] Inf Inf
.. .. ..@ max       : num [1:2] -Inf -Inf
.. .. ..@ unit      : chr ""
.. .. ..@ names     : chr [1:2] "slope" "aspect"
..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
.. .. ..@ type      : chr(0)
.. .. ..@ values    : logi(0)
.. .. ..@ color     : logi(0)
.. .. ..@ names     : logi(0)
.. .. ..@ colortable: logi(0)
..@ title   : chr(0)
..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
.. .. ..@ xmin: num 86799
.. .. ..@ xmax: num 86997
.. .. ..@ ymin: num 171858
.. .. ..@ ymax: num 172258
..@ rotated : logi FALSE
..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
.. .. ..@ geotrans: num(0)
.. .. ..@ transfun:function ()
..@ ncols   : int 3
..@ nrows   : int 184
..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
.. .. ..@ projargs: chr "+init=epsg:31300 +proj=lcc +lat_1=49.83333333333334 +lat_2=51.16666666666666 +lat_0=90 +lon_0=4.356939722222222 +x_0=150000.012"| __truncated__
..@ history : list()
..@ z       : list()
``````

slope_asp\$slope returns:

``````class       : RasterLayer
dimensions  : 184, 3, 552  (nrow, ncol, ncell)
resolution  : 65.93271, 2.17544  (x, y)
extent      : 86799.25, 86997.04, 171858, 172258.3  (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:31300 +proj=lcc +lat_1=49.83333333333334 +lat_2=51.16666666666666 +lat_0=90 +lon_0=4.356939722222222 +x_0=150000.01256 +y_0=5400088.4378 +ellps=intl +towgs84=-106.8686,52.2978,-103.7239,0.3366,-0.457,1.8422,-1.2747 +units=m +no_defs
data source : in memory
names       : slope
values      : NA, NA  (min, max)
``````

The story is the same for aspect...

The problem is now solved when working with the entire dataset (I had to define the nrow and ncol arguments). I'm now trying to create a loop to obtain slope and aspect values for each one of my 143 quadrats. Should I write ncol=1 and nrow=1 then? There seems to be a problem again...

``````data = read.table("Coordinates_elevation.txt", sep = "\t")

results = as.data.frame(matrix(nrow = 143, ncol = 3))
row.names(results) = c(1:143)
colnames(results) <- c('Elevation', 'Slope', 'Aspect')

for (i in 1:143){
as.factor(data\$V1)
as.factor(i)
subdata=data[data\$V1==i,]
X = subdata[,2]
Y = subdata[,3]
Z = subdata[,4]

subdata1 = matrix(c(X,Y,Z),  ncol=3,  byrow=FALSE)
e = extent(subdata1[,1:2])
r=raster(e, nrow=1, ncol=1, crs = CRS("+init=epsg:31300"))

x = rasterize(subdata1[,1:2], r, subdata1[,3], fun=mean)
slope_asp=terrain(x, opt=c('slope', 'aspect'), unit='degrees', neighbors=8)

#when trying with a single i, I'm again getting NA values for slope and aspect...

slope = slope_asp\$slope
results[i,2] = slope
aspect=slope_asp\$aspect
results[i,3] = aspect
elevation = subdata1[,3]
results[i,1]=elevation
}

#when running the entire loop, I'm getting "Error in `[<-.data.frame`(`*tmp*`, i, 2, value = <S4 object of class "RasterLayer">) : incompatible types (from S4 to logical) in subassignment type fix".

write.table(results, "results_slopeaspect.txt", sep = "\t")
``````
– Midavalo
Commented Feb 27, 2017 at 16:25
• Looking at your raster dimension, I wonder if you really wanted to create (3 x 184) raster grids over the (200m x 400m) extent. Each cell is very wide (x=65.9m) and short (y=2.18m). Could you check your raster() inputs? Commented Feb 27, 2017 at 16:44
• @Kazuhito In fact I have much more "coordinates-elevation" data (more than 9000) but it was really computationnaly intensive so I'm working with only 184 data (borders of my quadrats of 12.5x12.5 meters) while elaborating a correct script. Maybe the problem is that my plot isn't of rectangular shape ? Commented Feb 27, 2017 at 17:02
• @Kazuhito There is maybe a way to merge my script with a command making R understand the shape of the terrain with the expand.grid() function? Commented Feb 27, 2017 at 17:03
• @JeffreyEvans The unit argument in the terrain() function? Its description says it's a choice between 'degrees', 'radians' or 'tangent' (for slope and aspect calculation). Commented Feb 28, 2017 at 16:53

Simply use a SpatialPixelsDataFrame and the rasterFromXYZ function of the raster package to create the raster ( Creating a DEM from regularly / irregularly spaced points (R and Python) )

``````data = read.table("test.txt", h = T, sep = ",") # example with regularly spaced points
library(raster)
X = data\$x
Y = data\$y
Z = data\$z
data = matrix(c(X,Y,Z),  ncol=3,  byrow=FALSE)
e = extent(data[,1:2])
r=raster(e, ncol=3, nrow=25, crs = CRS("+init=epsg:31370"))
x = rasterize(data[,1:2], r, data[,3], fun=mean)
plot(x)
``````

And

``````slope_asp = terrain(x, opt=c('slope', 'aspect'), unit='degrees', neighbors=8)
summary(slope_asp)
[,1] [,2]
Min.      NA   NA
1st Qu.   NA   NA
Median    NA   NA
3rd Qu.   NA   NA
Max.      NA   NA
NA's      75   75
``````

2) with a SpatialPixelsDataFrame and the rasterFromXYZ function

``````df = data.frame(X,Y,Z)
theraster = rasterFromXYZ(df)
crs(theraster) = "+init=epsg:31370"
plot(theraster)
``````

Now you can compute a valid slope

``````slope_asp = terrain(theraster, opt=c('slope', 'aspect'), unit='degrees', neighbors=8)
summary(slope_asp)
[,1] [,2]
Min.      25  225
.....
``````
• Thank you! I tried with that function at first but it didn't worked because "x cell sizes are not regular" (message error). That's why I tried with rasterize() and terrain(). Commented Mar 1, 2017 at 16:00