[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