[R-sig-Geo] Speeding raster processing function
Robert J. Hijmans
r.hijmans at gmail.com
Sun Feb 16 02:01:08 CET 2014
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 at 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 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
>
> _______________________________________________
> 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