[R-sig-Geo] Speeding raster processing function

[Florian Betz]

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