Thanks all!   Great ideas (and new resources I wasn't aware of).

And thanks so much, Robert, brilliant as always.  This is much faster.  I'm
going run it to make sure I get the same outputs as before, but at first
blush it seems to working great.  I've a lot to learn re: speeding, and
hopefully studying this code will help matters.

To sum up, you basically transformed the non-na cells into a dataframe to
get the cell ID's, which were then used to generate a vector of all the x/y
coordinates from the distance argument (so something like 10x the number of
actual cells), extracted at once, and then (because it extracts in order)
were able to turn that into a matrix, from with you took the max for each
row (which corresponded to a given cell) and generated a raster based on if
that max was greater than the deflection angle.

That was exactly what I was looking for, I'll continue to play with it so I
understand it.  BTW, is there any reason you went with GTiff rather than
ascii?  I can't see how that would matter other than portability to other
systems and/or metadata.

Best,
Brian



On Sat, Feb 15, 2014 at 4:01 PM, Robert J. Hijmans <r.hijmans@gmail.com>wrote:

> 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.perpinan@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 <jgrn@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 <FloBetz@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: jgrn307@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
>

	[[alternative HTML version deleted]]

