5

i'm doing something quite straightforward regarding raster algebra but am struggling either to find the correct function or make a function work correctly;

Basically I have 2 rasters, representing consecutive years, that represent a classification (including NA values);

r <- raster(ncol=10,nrow=10)
r[] <- sample(c(1,2,4,8),size=100,replace=T)
r[runif(10*10) >= 0.50] <- NA

r1 <- raster(ncol=10,nrow=10)
r1[] <- sample(c(1,2,4,8),size=100,replace=T)
r1[runif(10*10) >= 0.50] <- NA

to get my change from one year (r) to the next (r1) I simply subtract r1 from r;

r2 = r - r1
freq(r2,digits=2)

Each new integer codes to a specific change, 0 simply means no change. I don't mind that cells become NA if NA only exists in one raster layer, this is fine. What i do want to do is to examine the 0 values more closely so if the original value in r is 1 and the recent value in r1 is also 1, i want to know this in the resulting r2 - ie quanitfy and analyse the 'no change' cells. Same for r=2 & r1=2, r=4 & r1=4 etc - not just all as 0s but as a collection of new, unused codes to represent what cells stayed the same (and how they stayed the same) from one year to the next.

Basically, i extracted the position of the cells from r2 that equal 0 as a spatial points data frame;

pts <-rasterToPoints(r2,fun=function(x){x==0},spatial=T)
plot(pts)

then i replaced the cells in r2 that equal 0 with values from r (or r1, doesnt matter as they are the same value) using the locations derived above;

r2[r2==0] <- (extract(r,pts))/10
plot(r2)
freq(r2,digits=2)

dividing the replacement values by 10 ensures they do not fall into the same 'bin' as any other mapped code change.

I am sure there is a quicker way? i thought i'd be able to create a raster stack from r and r1, creating an r2, then use some sort of 'where' or 'ifelse' function to perform the calculations above, however everything i do either results in errors or neglect of NA values to make the functions work. Or S4/integer errors.

I'm still not sure this will work on a big sample yet, and my final datasets are very massive.

8

What you want is a conditional calculation: return the value of r whenever r and r1 are equal and otherwise set the output to NA.

The cell-by-cell arithmetic operations seem to be fastest. (They are much faster than, say, using mask or the reclassification functions.) Since they do not appear to offer an actual conditional operator, use two time-honored tricks:

  1. Treat logicals as numbers. FALSE is 0 and TRUE is 1 in arithmetic operations.

  2. Create NA values (or, almost as effectively, infinite values) using invalid arithmetic operations.

One solution is

r3 <- r == r1
r3 <- r3 * r * (1/r3)

It works because when r and r1 are equal, both r3 and 1/r3 equal 1 and the multiplications change nothing: they return the value of r. When r and r1 are not equal, 1/r3 is undefined, producing an infinite result. As a result, freq tabulates only the cell values where r and r1 agree.

On my machine this calculation takes about one second for rasters with 10,000,000 cells. (It's about ten times as long as the simple comparison r - r1.) It will scale in direct proportion to the number of cells until disk paging is invoked, at which point you will be at the mercy of your storage throughput.

(If you can fit all data into RAM, it's even faster to use R's built-in operations on the array of data and then convert it back to a raster object.)

  • this looks like a great suggestion and thanks for the info re processing and memory. However, how does this incorporate those cells that are i need where the result is not a 0? ie where the initial value is 4 (for eg) and the next value is 8, giving a new cell value of -4? every new value represents a type of change, i need them all, but it is the 0s that i need breaking down - all in the same raster – Sam Aug 26 '15 at 7:45
  • i think i fixed it, i'll add an answer of my own to show the full working, using whuber's part solution, thanks again – Sam Aug 26 '15 at 11:55
  • 1
    It might be a little complicated and error-prone to try to mix the two kinds of results in a single raster, Sam. Why not compute one that tabulates the actual differences in r - r1 and a second one, as described here, that tabulates the common values in r and r1? That will include all the information you need and provide added cartographic flexibility (because you can symbolize the two results separately and then overlay them). – whuber Aug 26 '15 at 13:39
  • 1
    yes that did occur to me Whuber. Extensive tests have shown combining them to be reasonably quick on very large rasters and accurate too, but i think i like the idea of having 2 separate rasters for simplicity and, as you say, display flexibility. thanks again :) – Sam Aug 26 '15 at 14:59
1

Since whuber gave me the nice quick piece of code to create a raster with the values that are the same between rasters, i thought i'd finish the entire job;

create rasters and subtract one from the other to get a 'change' raster (full of 0s that need examining as well);

r <- raster(ncol=10,nrow=10)
r[] <- sample(c(1,2,4,8),size=100,replace=T)
r[runif(10*10) >= 0.50] <- NA

r1 <- raster(ncol=10,nrow=10)
r1[] <- sample(c(1,2,4,8),size=100,replace=T)
r1[runif(10*10) >= 0.50] <- NA

r2 = r - r1

then, using whuber's little piece of code, create a new raster where the values are the same (dividing by 10 so the cell values don't match anything calculated previously);

r0 <- r == r1
r0 <- (r0 * r * (1/r0))/10

To combine the original change raster with the new raster that analyses the cells that have no chnage, add them together, first changing NA values to 0;

r0[is.na(r0)] <- 0
rfinal <- r2 + r0

that creates a new raster with all the changes and all of the 0 values changed to a specific code to examine exactly what sort of 'no change' is occurring.

export to csv;

counts <- freq(rfinal,digits=2,useNA='no')
write.csv(counts,"counts.csv")

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.