[R-sig-Geo] Speeding raster processing function
Jonathan Greenberg
jgrn at illinois.edu
Sat Feb 15 19:52:44 CET 2014
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:
> 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 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]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at 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: jgrn307 at hotmail.com, Gchat: jgrn307, Skype: jgrn3007
More information about the R-sig-Geo
mailing list