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)
#data$V1 are my quadrat numbers
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")
Best Answer
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) )
1) With your solution
And
2) with a SpatialPixelsDataFrame and the rasterFromXYZ function
Now you can compute a valid slope