[R-sig-Geo] Speeding up raster calculations

[Jonathan Greenberg]

Franz:

Here is an example of a weighted by distance from center of the window
function using rasterEngine.  I'd recommend grabbing the latest
version from r-forge (1.3.4 at the time of this email).  I'm
calculating a distance matrix within the function that can be used to
create the weights.  You'd probably want to make the distance matrix
once outside of the function and feed it in as an args=.

library(spatial.tools)

focal_function_with_distance <- function(resource_map,inverse_distance_weight)
{

# You can actually do this calculation once outside of the function,
and feed it in as a weight.
# Create a vector of x-distances from the center of the window:
x_vector <- -floor(dim(resource_map)[1]/2):floor(dim(resource_map)[1]/2)
# Create a vector of y-distance from the center of the window:
y_vector <- -floor(dim(resource_map)[2]/2):floor(dim(resource_map)[2]/2)
# Now turn them into matrices:
x_matrix <- matrix(data=x_vector,ncol=dim(resource_map)[1],nrow=dim(resource_map)[2],byrow=TRUE)
y_matrix <- matrix(data=y_vector,ncol=dim(resource_map)[1],nrow=dim(resource_map)[2],byrow=FALSE)
# Now create a window of the same size as the local window that
contains the distance to the
# center of the window:
distance_matrix <- sqrt(x_matrix^2 + y_matrix^2)

# Some random weighting function:
# Avoid dividing by 0:
weighted_matrix <- 1/(distance_matrix+0.00001)^inverse_distance_weight

# Calculate the mean of the weight times the local window:
weighted_mean <- mean(resource_map*array(weighted_matrix,dim=dim(resource_map)))

# Return the value:
return(weighted_mean)
}



tahoe_lidar_highesthit <-
  raster(system.file("external/tahoe_lidar_highesthit.tif",
package="spatial.tools"))

# Turn on the parallel engine:
sfQuickInit()
# Use rasterEngine.  Turn debugmode=TRUE if you want to debug your
function on a single window.
# This uses a 5 x 5 window:

weighted_by_distance_raster <-
rasterEngine(resource_map=tahoe_lidar_highesthit,fun=focal_function_with_distance,
args=list(inverse_distance_weight=2),window_dims=c(5,5),debugmode=FALSE)

# When you are done it is usually a good idea to shut down the engine
unless you plan
# on doing more rasterEngine runs:
sfQuickStop()

On Mon, Mar 3, 2014 at 8:22 PM, Jonathan Greenberg <jgrn at illinois.edu> wrote:
> Franz:
>
> I'm not sure what you mean that you can't use focal_hpc for your
> problem -- it seems like your problem would be fine for use with that
> function -- its a bit easier to use rasterEngine, by the way.  You can
> use > 1 input raster into your functions if you are using
> rasterEngine, and you can pass non-raster parameters to the function
> as well.  If you had a simple set of data to test, perhaps I can help
> you write the function.
>
> As a heads up, focal window algorithms are slow, and get slower the
> larger the window size.  rasterEngine/focal_hpc in spatial.tools will
> use parallel processing to speed it up, however.
>
> --j
>
> --j
>
> On Mon, Mar 3, 2014 at 4:27 PM, Franz Grassmann
> <franz.grassmann at gmail.com> wrote:
>> Dear list members!
>>
>> Maybe, someone has a good idea how to speed up my raster-calculation, or knows the right way how to deal with my "problem".
>>
>> For my master thesis, I use the focal-function and its focalWeight ("raster"-package) for some raster calculations. For every cell of my map, I sum up the "resources" of the surrounding cells by a certain weight that declines with distance (Its about foraging bees...). I got a map that holds the resources, and a certain focalWeight for every species of my model.
>>
>> The command code looks similar to this:
>>
>> result.raster <- focal(x = resource.map,
>>                           fun = function(x) {sum(x)/sum(fW)},
>>                           w = fW)
>>
>> It takes very long to create my raster, because the fW is quiet large for some species. I can't use the focal_hpc()-command, because it doesn't support the kind of focalWeight I use. For my analysis I also have to run this calculations for a number of species, seasons and scenarios.
>>
>> My questions: Did I choose the right way to solve my tasks, or is there another, faster way for my calculations? I hope I didn't choose the wrong package or program...
>> I am no programmer, but I am very interested in this topic an would welcome every tip on how to solve this, or hint where to read/learn more about this.
>>
>> Best regards,
>> Franz Grassmann
>> _______________________________________________
>> 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



-- 
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