Brian, Here is, I think, an improved version. Your code had many steps that could be simplified; but removing the inner loop was key, as you indicated. For the one example I used, this one is > 150 times faster. More could be done (e.g. by using data.table) and it could be made memory safe for large rasters. Robert
library(raster) windout <- function(dem, deflect, angles, max.dist) { #note that the dem raster must be in planar coordinates #for smaller datasets: dem <- readAll(dem) res <- res(dem) # you are ignoring 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() for (angle in angles) { elev <- na.omit(data.frame(ID=1:ncell(dem), value=getValues(dem))) coords <- xyFromCell(dem, elev$ID) radangle <- (angle+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$value comp.elev <- t(t(comp.elev) / distance) notAllNA <- rowSums(is.na(comp.elev)) != num.dist ang <- atan(comp.elev[notAllNA, ]) * (180 / pi) output <- rep(NA, ncell(dem)) output[elev[notAllNA,1]] <- apply(ang, 1, max, na.rm=TRUE) r <- setValues( dem, output <= deflect ) fname <- paste("angle", angle, deflect, "search", max.dist, sep="_") result[as.character(angle)] <- writeRaster(r, fname, format="GTiff", overwrite=TRUE) } result } f <- system.file("external/test.grd", package="raster") r <- raster(f) system.time( x <- windout(r,1,c(90,180),300) ) plot(stack(x)) On Sat, Feb 15, 2014 at 11:58 AM, Oscar Perpiñan <oscar.perpi...@gmail.com> wrote: > Hello, > > Although it is a different problem, maybe our solution for downscaling > of global solar irradiation is useful for you. It is documented in > this arXiv paper: http://arxiv.org/abs/1311.7235. The code is here: > https://github.com/EDMANSolar/downscaling/blob/master/code.R (horizon > blocking with parallel computing begins at line 263). > > Best, > > Oscar. > > ----------------------------------------------------------------- > Oscar Perpiñán Lamigueiro > Grupo de Sistemas Fotovoltaicos (IES-UPM) > Dpto. Ingeniería Eléctrica (ETSIDI-UPM) > URL: http://oscarperpinan.github.io > Twitter: @oscarperpinan > > > 2014-02-15 19:52 GMT+01:00 Jonathan Greenberg <j...@illinois.edu>: >> 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 <flob...@web.de> wrote: >>> Hi Brian, >>> >>> maybe you could also have a look at the wind.shelter function in the RSAGA >>> package? It gives you an index for the exposure/shelter of a surface to the >>> wind and works quite well. >>> Of course, this does not really solve your code problem, but might be a >>> faster alternative... >>> >>> Regards, >>> >>> Flo >>> >>> >>> >>>> 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 >>> >>> >>> _______________________________________________ >>> R-sig-Geo mailing list >>> R-sig-Geo@r-project.org >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> >> >> >> -- >> Jonathan A. Greenberg, PhD >> Assistant Professor >> Global Environmental Analysis and Remote Sensing (GEARS) Laboratory >> Department of Geography and Geographic Information Science >> University of Illinois at Urbana-Champaign >> 259 Computing Applications Building, MC-150 >> 605 East Springfield Avenue >> Champaign, IL 61820-6371 >> Phone: 217-300-1924 >> http://www.geog.illinois.edu/~jgrn/ >> AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007 >> >> _______________________________________________ >> R-sig-Geo mailing list >> R-sig-Geo@r-project.org >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo > > _______________________________________________ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo