Quantifying wind exposure in R

Exposure to storm winds is important for ecological communities - it shapes the pattern of seed dispersal, wind disturbance/blowdown, fire behavior, soil turnover (via tip up mounds), and many other processes.  Quantifying exposure over large areas is tricky, though.  

I put together some code to calculate wind exposure as a function of topographic shading and known storm directions - this is for straight-line winds, primarily over broad extents, and was originally developed in the 1990's for hurricanes.  So it works well on those scales.  This is exceptionally valuable in complex terrain - like Colorado or Alaska.

Top:  Landscape shot.  Middle:  DEM of that landscape.  Bottom:  Calculated wind exposure (red high, blue low).

To calculate exposure, you basically take an angle of wind, and allow that wind to come into your landscape, bend over topographic barriers, and calculate relative exposure.  Areas that are flat have high exposure, of course, as do hills facing the incoming wind direction - if they are not shaded by something upwind.  In real life, of course, you'll have a distribution of wind directions - so calculate exposure from that distribution and weight accordingly (e.g., by relative frequency).  This gives a nice relative score to your landscape.

The code is accomplished via two snippets:


First, the following code quickly calculates - based on user supplied wind direction, deflection angles, and search distances (how far upwind a barrier should matter) - relative exposure:

####################

#    Function

#####################

windout.iter <- function(dem, deflect, angles, max.dist) {
#note that the dem raster must be in planar coordinates

#for smaller datasets:
#dem <- readAll(dem)     #if in memory
res <- res(dem)

# do not ignore the case where x and y resolution are not equal
stopifnot(all.equal(res[1], res[2]))
xr <- res[1]

#number of distances to check, basically goes every other cell for speed.
num.dist <- round(max.dist / xr / 2)
distance <- seq(xr, max.dist, length.out=num.dist)
result <- list()
j <- 1

for (d in deflect) {
    midrow <- cellFromRow(dem,rownr=1)    #note this does the top one
    elev <- extract(dem,midrow)
    coords <- xyFromCell(dem, midrow)

    radangle <- (angles+90) * pi/180  #convert to radians.
    dcosangle <- -cos(radangle) * distance
    dsinangle <- sin(radangle) * distance
    x <- apply(coords[,1,drop=FALSE], 1, function(j) j + dcosangle)
    y <- apply(coords[,2,drop=FALSE], 1, function(j) j + dsinangle)
    xy <- cbind(as.vector(x), as.vector(y))

    comp.elev <- extract(dem, xy)
    comp.elev <- matrix(comp.elev, ncol=num.dist, byrow=TRUE)
    comp.elev <- comp.elev - elev
    comp.elev <- t(t(comp.elev) / distance)
    #notAllNA <- rowSums(is.na(comp.elev)) != num.dist
    ang <- atan(comp.elev) * (180 / pi)

    r <- apply(ang,1,max)

    r <- r<=d
    result[[j]] <- r*1
    j <- j+1
    }

output <-simplify2array(result)
output <- apply(output,1,sum)
output <- output+1
outputs <- list()
outputs[[1]] <- output
outputs[[2]] <- coords
return(outputs)
}
 


This next code then loops through a DEM, first subsetting out an area the size of the max distance (that keeps things fast) and then calculating all wind directions desired and averaging the results.  To do this in practice, one needs to know the distribution of storm-force wind directions, and builds the directions based off that (to weight a given direction, one could either calculate twice or simply duplicate).  The following works at a 5km distance, several deflection angles, and is oriented around S, SE, and SW wind:

########

library(raster)

#load single big DEM
dem <- raster(file.choose())
    plot(dem,maxpixels=10000)
    t <- drawExtent()
    dem <- crop(dem,t)   #If you want to do a focal area

#set projection
proj.def <- "+proj=utm +ellps=WGS84 +zone=8 +units=m"  #Adjust as needed.
dem <- projectRaster(dem, crs=proj.def)

storage <- matrix(nrow=nrow(dem),ncol=ncol(dem))

#set parameters
max.dist <- 5000
deflect <- c(1,3,5,7,9,11,13,14)
angles <- c(135,180,225)

iter <- 1:nrow(dem)    #to 3772 now
r <- res(dem)[1]


for (i in iter) {
    temp.extent <- extent(dem)
    temp.extent@ymax <- temp.extent@ymax-(i*r)    #*res(dem)[1] to avoid top
    temp.extent@ymin <- temp.extent@ymax-(i*r+max.dist+r)

    temp.dem <- crop(dem,temp.extent)
    temp1 <- windout.iter(temp.dem,deflect,angles[1],max.dist)[[1]]
    temp2 <- windout.iter(temp.dem,deflect,angles[2],max.dist)[[1]]
    temp3 <- windout.iter(temp.dem,deflect,angles[3],max.dist)[[1]]

    #temp.coords <- SpatialPoints(temp.coords)
    temp <- apply(cbind(temp1,temp2,temp3),1,mean)
    #t.loc <- cellFromXY(storage,temp.coords)

    storage[i,] <- temp    
    print(i)
    removeTmpFiles(h=0)   #This is needed large processing jobs, which crash otherwise.
    rm(temp.dem)
    gc()
}

gc()
t <- dem   #this creates a place to put the calculated values
t[] <- storage

par(mfrow=c(1,2))   #look at some comparisions
plot(t)
plot(dem)

writeRaster(t,file.choose(),format="GTiff",overwrite=T)


###################################################################3