[R-sig-Geo] Speeding raster processing function

Oscar Perpiñan oscar.perpinan at gmail.com
Sat Feb 15 20:58:44 CET 2014


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 at 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 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
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo



More information about the R-sig-Geo mailing list