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]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
