0

Background I created some digital elevation models and subsequently made a 2 layer raster for each containing the calculated aspect and slope using raster::terrain. Now I want to calculate the area of each which satisfies these two criteria:

  • 135 < Aspect degrees < 215
  • Slope degrees > 2

Here is a link to a simple example of my dataset https://www.dropbox.com/sh/5mfp421lg9jwkkh/AACG-_pTh56ie6zYvf6RRmRHa?dl=0

What I tried I load both layers simultaneously as a stack. I first tried to do a simpler selection just to get a hang of the notation (being a beginner in R) on a sample dataset (I know its outside my own parameters set out above, but I just want to get it right first)

InvestigatedArea <- stack(paste0(getwd(),"/2150913_AspectSlope1.grd"))
result1 <- InvestigatedArea$slope[InvestigatedArea$slope[] ==  "2.336640"]
result2 <- InvestigatedArea$aspect[InvestigatedArea$aspect[] ==  "83.33418"]

Both these result in NULL answers, but I know my sample data has these values in them. Visible below in slots 81 in both instances.

> InvestigatedArea$slope[]
  [1]        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA
 [11]        NA        NA        NA        NA        NA        NA        NA        NA 15.014699 27.780954
 [21] 28.800718 21.376579 14.421870  9.720741  9.799111 16.171215 22.605927 25.304127 23.419540 18.468254
 [31] 14.223759 11.670895  9.239475        NA        NA 24.906059 32.080647 27.778091 17.200335 11.943377
 [41] 12.618942 17.589796 24.648991 27.494242 26.727062 22.572401 15.775166 10.440066  8.643677  6.733919
 [51]        NA        NA 30.189089 30.811302 19.766258  8.661880  9.401229 16.495119 24.468115 28.694721
 [61] 25.355656 20.456161 14.836404  9.168124  5.176007  4.732281  4.400262        NA        NA 26.651682
 [71] 24.032007 12.866996  5.050277 10.976402 20.835245 29.398857 29.420214 19.465899 10.608719  4.032556
 [81]  2.336640  2.177207  2.326176  3.879496        NA        NA        NA        NA        NA        NA
 [91]        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA
[101]        NA        NA

> InvestigatedArea$aspect[]
  [1]        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA
 [11]        NA        NA        NA        NA        NA        NA        NA        NA 307.59409 300.68588
 [21] 298.03110 296.04544 290.26505 293.71802 315.41357 322.72733 321.72003 318.79593 318.05820 319.33633
 [31] 316.50836 313.69717 320.50980        NA        NA 298.20007 303.81354 310.67441 317.35971 305.58899
 [41] 300.54056 309.42044 316.02017 322.52969 325.32404 329.24329 337.43942 330.18988 310.18207 309.40103
 [51]        NA        NA 298.71249 295.38638 302.32114 321.94833 303.58517 296.89182 301.38943 307.88290
 [61] 318.19205 322.84161 329.89163 358.89279 357.97546 296.30090 277.03336        NA        NA 299.59601
 [71] 284.44580 273.19974 277.92230 291.86438 292.40649 294.21609 300.87869 316.21036 322.60995 311.60651
 [81]  83.33418  89.59325 248.77850 236.84702        NA        NA        NA        NA        NA        NA
 [91]        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA
[101]        NA        NA

To be honest, I'm not even sure I am pointing to the right position where the actual points are stored. If I look at the structure of the rasterstack in the RStudio Environment browser, I cant actually identify the place where the data values are stored. I do know there is data, because I can plot it: Simple plot() of "InvestigatedArea"

So I cant even begin to calculate specific areas, because I can't even get past this first hurdle. I know it's probably just a stupid notation mistake. What must I change in my notation to get a list of positions that contain the data that satisfies my conditions?

  • This InvestigatedArea$slope[] == "2.336640" is trying to compare character values, not numeric ones. Its like testing if sqrt(2) == "1.414214" - that's FALSE. Try something like InvestigatedArea$aspect[] < 215 as the expression. Exact comparison of decimal numbers is ill-advised, and is R FAQ 7.31 cran.r-project.org/doc/FAQ/… – Spacedman Jun 7 '18 at 11:58
  • @Spacedman you are absolutely right (an amateur mistake) and if I went straight to my original requirements, instead of trying to test a simple case first, I would have avoided the problem altogether. :-( Now, I at least get a result, but its a vector that contains all the values from the original that satisfy the condition. What I actually want is the position number where those values were found. Any ideas? – Waldo vdM Jun 7 '18 at 12:21
  • Sounds like you want like which(r[]>2) – Spacedman Jun 7 '18 at 13:20
1

Here's how to get the area of a raster that matches some condition. First lets set up some data:

> r = raster(res=c(10,10))
> r
class       : RasterLayer 
dimensions  : 18, 36, 648  (nrow, ncol, ncell)
resolution  : 10, 10  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 

This is a 10-degree square raster in geographic coordinates (lat-long) so it has 648 numbers in it. Lets put random numbers from 0 to 1 in those cells:

> r[]=runif(648)

How many of those cells are over and under 0.99? Compute a test on the values, and add up all the TRUE or FALSE values, since TRUE is 1 and FALSE is 0:

> sum(r[]>0.99)
[1] 8
> sum(r[]<=0.99)
[1] 640

So 8 cells are over 0.99. But the true area of the cells varies as we go towards the poles, because this is a lat-long grid! We can use the area function to get a matching raster with the area in. Try:

> plot(area(r))

to see. What we need to do to compute the total area of those 8 cells is to add up the area raster for those 8 cells.

> sum(area(r)[r[]>.99])
[1] 6604824

which is in square metres.

  • Thanks for taking the time to answer my questions (even if it is really just silly mistakes) but even more so, thank you for looking at the bigger picture and answering the true question about calculating areas (even though I know its only supposed to be one question with one answer). I think I am getting a better grip on the problem and how to solve it. I'm toying with the idea of using mask() and/or reclassify() as alternatives to test the results against each other. – Waldo vdM Jun 7 '18 at 13:47
  • Im just shooting ideas here, but to add to your answer I could also spTransform() to an equal area projection, then just add together the area of any 8 arbitrary cells. – Waldo vdM Jun 7 '18 at 13:52
  • spTransform only works on vector data (doesn't it?). You could warp your raster to a different coordinate system but that's nearly always an irreversible and lossy process that adds errors to your calculation. – Spacedman Jun 7 '18 at 14:55

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.