2

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")
  • @Roxane Please use the edit link below your question to add any extra info in response to comments. The answer area is reserved for answers, not for extra information. – Midavalo Feb 27 '17 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? – Kazuhito Feb 27 '17 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 ? – Roxane Feb 27 '17 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? – Roxane Feb 27 '17 at 17:03
  • 1
    @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). – Roxane Feb 28 '17 at 16:53
2

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

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)

enter image description here

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) 

enter image description here

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(). – Roxane Mar 1 '17 at 16:00

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.