Hey all, I have a challenge that I'm fairly stuck on.  There is likely a
fast solution, but I don't know what it is.  I'm trying to create a fast
function to calculate exposure to wind based on topography (similar to
Boose et al.).  It's conceptually simple, in that it's similar to a
viewshed - if you're sheltered from a higher elevation down wind, you
aren't exposed.  But, the wind is allowed to "bend" down at a given angle.
So I need to calculate exposure to wind that bends down a bit (like a
viewshed, but the light can bend down a tad).  That bending distance is the
"threshold."

I've got a function that works, but it's very slow.  I know why - it's slow
because it's running cell by cell.  But I'm at a loss as to how to
calculate exposure for the whole raster at a time - mainly because of that
bending.

I've tried translating the raster via the shift function, but haven't
gotten that to work.  Currently, I have a function to extract the
elevations over a line at a given angle/distance.  Then I calculate the
angles, and if any angle is larger than the threshold angle, the point is
sheltered.  If not angles are less, it's not.

There must be a faster way, however.  I imagine there's a clever way to
take advantage of rasters decent computing speed, but my function is the
bottleneck.  Any thoughts?  I've tried working on this for a while, and I'm
at the point where it's just fast enough for what I need, if I run it over
the weekend (64bit, 8mb ram, i7 intel).  Code is pasted below.  I've put
off posting this here hoping to figure it out myself, but I'm not making a
lot of headway, and it's easier to explain by posting code than anything
else.

Thanks in advance!


#Code#

###################
require(rgeos)
library(raster)
library(rgdal)

#parameters for later
angle.list <- c(90,180)   #180is south
dem <- dem    #any dem in UTM coordinates

#function to extract location at given angle and distance.  0 is due
north.  angle in degrees.  I found this online somewhere- would love to
attribute it to the right person but don't remember where I found it.

####
retPoCordinates <- function(center,angle,distance){
    angle <- (angle+90)*.0174532925    #convert to radians.
    x <- round(center@coords[1,1] + distance * cos(angle),0)
    y <- round(center@coords[1,2] + distance * sin(angle),0)
    result <- as.data.frame(cbind(x,y));names(result) <- c("x","y")
    return(result)
    }

################
#Main function
#
################
out <- function(dem, deflect, angle.list,max.dist) {
#note that the dem raster must be in UTM

#direction
num.angle <- length(angle.list)
start.ang <- 1

while (start.ang <= num.angle) {                    #for multiple angles,
single deflection degree

    angle <- angle.list[start.ang]

    num.dist <- round(max.dist/res(dem)[1]/2)            #number of
distances to check, basically goes every other cell for speed.
    distance <- seq(res(dem)[1],max.dist,length.out=num.dist)

    #temp storage
    dist.store <- matrix(nrow=length(distance),ncol=1)
    output <- dem*0
    output <- as.matrix(output)
    output <- as.vector(output)

    #initial elevation vector
    temp.elev <- extract(dem, 1:ncell(dem), df = TRUE)
    temp.elev[is.na(temp.elev)==T] <- 0

    #initial xy coordinates vector
    cells <- temp.elev$ID
    cells <- temp.elev$ID[which(temp.elev[,2]>0)]
    coords <- xyFromCell(dem,cells,spatial=T)

    #for each cell in matrix
    index <- length(cells)
    s <- 1
    d <- res(dem)[1]

    while (s <= index) {                             #bottleneck.

        i <- cells[s]

        #remove NA
        g.elev <- ifelse(is.na(temp.elev[i,2]),0,temp.elev[i,2])

        #if g.elev = 0,
        if (g.elev == 0) {
            temp.out <- -1
            } else {

            #calculate location at distance
            out <- retPoCordinates(coords[s],angle,distance)

            out <- SpatialPoints(out,proj4string=CRS(projection(dem)))
            comp.elev <- extract(dem,out)

            diff <- comp.elev - g.elev
            temp <- diff/distance
            ang <- atan(temp)*57.2958

            #elminate NA's
            ang <- ifelse(is.na(ang)==T,-100,ang)

            #testing
            #print(ang)
            #plot(out,add=T)

            temp.out <- ifelse(max(ang)>deflect,0,1)
            }

        output[i] <- temp.out
        print(s/length(cells))
        s <- s+1
    }

    #save output
    t <- matrix(output,nrow=nrow(dem),ncol=ncol(dem),byrow=T)
    t <- raster(t)
    projection(t) <- projection(dem)
    extent(t) <- extent(dem)
    nam <- paste("angle",angle,deflect,"search",max.dist,sep="_")
    assign(nam,t)

    temp <- get(nam)
    writeRaster(temp,nam,format="ascii",overwrite=T)
    rm(temp)
    gc()


    start.ang <- start.ang + 1
    }
}

#to execute
out(dem,1,angle.list,3000)   #this would calculate the exposure with a
threshold of 1 degree bending (so not much), for the angle list (3 angles),
and search out 3000m

	[[alternative HTML version deleted]]

