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