Hi Brian:
Couple of things: first, viewshed analyses are notoriously complicated
and result in decidedly non-embarassingly parallel problem. A lot of
the most efficient solutions I've seen involve GPU processing.
With that said, you might want to try to port this to my rasterEngine
framework (part of the "spatial.tools" package), where you can define
a local window size and write your function so it processes the output
for that single window (the maximum distance to search for a horizon
would be the window size). rasterEngine will then run that function
in parallel (use sfQuickInit(cpus=however many cpus you have on your
computer)), alleviating some of the common slow-downs with
local-window processing (e.g. it will only read a pixel from your
input ONCE -- a lot of local window analyses that aren't optimized
will end up re-reading regions over and over again).
I have a tutorial at:
http://publish.illinois.edu/jgrn/software-and-datasets/rasterengine-tutorial/
I have some examples of local-window analyses on that page. I haven't
updated the tutorial for use with some of rasterEngine's new features
you may want to also use, e.g. debugging, byte-compiling your
function, and multiple-file outputs, but you can look at those
features in the help (?rasterEngine) if you grab the latest from
r-forge:
http://r-forge.r-project.org/R/?group_id=1492
--j
On Sat, Feb 15, 2014 at 4:10 AM, Florian Betz <FloBetz at web.de> wrote:
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 at coords[1,1] + distance * cos(angle),0)
y <- round(center at 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]]