mercredi 11 juillet 2018

How to translate GRASS code for southness and westness into R code?

This is my first time posting so let me know how to improve the question if need be.

I have spent a long time searching online for help and asked lab mates but I'm up against a wall. See below for code to calculate continuous westness and southness variables from elevation raster. This code was written for GRASS, which I am unfamiliar with. I would like to keep all coding in R (I'm completing a large project, many coding components will later be shared with other R users). The original code writer is very difficult to get a hold of to ask questions which is why I am here.

Original code sent to me is:

system('r.in.gdal.exe -oe input=C:/DEM.tif output=dem')     #0 or 1
execGRASS('r.slope.aspect', parameters=list(elevation='dem', slope='slope1', 
          aspect='aspect_grass'))   # 0=East, 90=North, 180=West, 270=South, 
system('r.mapcalc.exe "slope= round(slope1)" ')
system('r.mapcalc.exe "aspect = round((450 - aspect_grass) % 360)" ')
system('r.mapcalc.exe  "aspect_s1 = if(slope <= 1, 0, -1 * 
       sin(aspect_grass))"')  # uses degrees
system('r.mapcalc.exe  "aspect_w1 = if(slope <= 1, 0, -1 * 
       cos(aspect_grass))"')
system('r.mapcalc.exe "aspect_s = round(100 * aspect_s1)" ')
system('r.mapcalc.exe "aspect_w = round(100 * aspect_w1)" ')

*** Note - I am confused by the fourth command - I think aspect = round((450 - aspect_grass) % 360) might be unnecessary for this operation?

This is my version so far in R for southness but the 'if' statement when calculating southness and westness has stumped me (I believe it's C++ format above?). Slope and aspect rasters have been created by me for the same study area, calculated from the same DEM raster using the raster::terrain function.

slope <- raster('H:/slope.img')
aspect <- raster('H:/aspect.img')

slope_round <- round(slope)
aspect_s1 <- if(slope_round <=1) {
(0)
} else {
(-1 * sin(aspect))
}
aspect_s <- round(100 * aspect_s1)

The if statement returns this error:

Error in if (slope_round <= 1) { : argument is not interpretable as logical

I am not yet comfortable with R 'if' statements and am not sure I'm interpreting the C++ code correctly. Would very much appreciate some assistance with this translation!

Aucun commentaire:

Enregistrer un commentaire