# Chapter 9 Geometric operations with rasters

*Last updated: 2024-06-24 11:08:31 *

## Aims

Our aims in this chapter are:

- Make changes in the geometric component of rasters:
- Mosaicking
- Resampling
- Reprojecting

- Apply focal filters on a raster

We will use the following R packages:

`stars`

`units`

## 9.1 Mosaicking rasters

In the next few examples, we will prepare a Digital Elevation Model (DEM) raster of Haifa, by mosaicking, subsetting and reprojecting (Section 9.3).

We start with two \(5°\times5°\) tiles of elevation data from the Shuttle Radar Topography Mission (SRTM) dataset. The tiles are included as two `.tif`

files in the sample files (Appendix A):

As shown in Figure 9.1, the tiles cover the area of northern Israel, including Haifa:

```
plot(dem1, breaks = "equal", col = terrain.colors(10), axes = TRUE)
## downsample set to 5
plot(dem2, breaks = "equal", col = terrain.colors(10), axes = TRUE)
## downsample set to 5
```

Using the `st_bbox`

and `dim`

functions (Section 5.3.8.3), we can see that the two rasters indeed comprise two tiles of the same “large” dataset. First, we can see that the dimensions (number of rows and columns) of the two rasters are identical:

Second, we can see that their extents are aligned. The y-axis (i.e., latitude) extent of the rasters is the same (`30`

-`35`

). The x-axis (i.e., longitude) extents of the rasters are side-by-side (`30`

-`35`

and `35`

-`40`

):

Raster datasets with global coverage (or other large areas) are commonly distributed in the form of tiles. While tile dimensions and placement may differ in each case, the general idea of tiling is the same: reducing the amount of data that the user needs to deal with. Rather than have the user download a huge (e.g., global) raster file and crop their area of interest, which can be challenging and time-consuming due the file size, the user can download only the required tiles, mosaic them, and then crop the specific area of interest. That way, the user doesn’t have to download and deal with huge files where most of the data are irrelevant.

Rasters can be **mosaicked** using the `st_mosaic`

function. The `st_mosaic`

function accepts two or more `stars`

objects—such as `dem1`

and `dem2`

—and returns a combined raster:

Note that the `st_mosaic`

function can be used to combine rasters even if they are not aligned (unlike `dem1`

and `dem2`

, in this example, which are aligned). In case the rasters are not aligned, resampling (Section 9.2) takes place.

The mosaicked DEM is shown in Figure 9.2:

The extent of the resulting raster covers the extents of both inputs:

For the next examples, we “crop” the `dem`

raster according to an extent of \(0.25°\times0.25°\) around Haifa (Figure 9.3).

The `[`

operator can be used to crop the raster (Section 6.2):

The result is shown in Figure 9.4:

Note that in this example, we needed to know in advance the row and column indices of the extent we are interested in cropping. Later on (Section 10.1), we are going to learn how to crop a raster based on an existing vector layer, which is more practical.

## 9.2 Raster resampling

### 9.2.1 The `st_warp`

function

Raster **resampling** is the process of transferring raster values from the original grid to a different grid (Figure 9.5). Resampling is often required for:

- Aligning several input rasters that come from different sources to the same grid, so that they can be subject to spatial operators such as raster algebra (Sections 6.4 and 6.6.1)
- Reducing the resolution of very detailed rasters, so that they are more convenient to work with in terms of processing time and memory use

To demonstrate resampling, we will create a custom `stars`

grid, using the same extent as the `dem`

raster (`st_bbox(dem)`

), but a coarser resolution. The new resolution is \(0.002°\), which is ~2.4 times the original resolution of \(0.00083°\). The grid is created using `st_as_stars`

:

```
grid = st_as_stars(st_bbox(dem), dx = 0.002)
grid
## stars object with 2 dimensions and 1 attribute
## attribute(s):
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## values 0 0 0 0 0 0
## dimension(s):
## from to offset delta refsys x/y
## x 1 251 34.74 0.002 WGS 84 [x]
## y 1 251 33.04 -0.002 WGS 84 [y]
```

Recall that we already used the `st_as_stars`

function to convert a `matrix`

to a `stars`

raster (Section 6.3.3). What we do here is another mode of operation of the `st_as_stars`

function, creating an *empty* raster given a `bbox`

and the resolution `dx`

.

We can resample a raster using the `st_warp`

function. The first two parameters of `st_warp`

are:

- The raster to be resampled, where the values come from
- The raster defining the new grid

Here is the expression to resample the values of the `dem`

raster into the new raster `grid`

:

Which one of the rasters in Figure 9.5 represents the role of

`dem`

and`grid`

, in our case?

Try resampling

`dem`

to a different grid with an even coarser resulution, such as`0.02`

, and plot the result

The original DEM is shown along with the resampled one in Figure 9.6.

```
plot(dem, breaks = "equal", col = terrain.colors(10), main = "input (delta=0.0008)", key.pos = 4)
plot(dem1, breaks = "equal", col = terrain.colors(10), main = "output (delta=0.002)", key.pos = 4)
```

When resampling, the computer needs to decide which pixel value(s) to transfer to each of the “new” pixels, defined by the new grid. There are several possible options, known as resampling algorithms, or resampling methods. In the next three sections, we are going to demonstrate three common ones:

### 9.2.2 Nearest neighbor resampling

The `st_warp`

expression from Section 9.2.1 actually used the **nearest neighbor** resampling method (`method="near"`

), which is the default:

To understand what actually happens in nearest neighbor resampling, let’s take a look at a small part of the raster “before” and “after” images (Figure 9.7). If you look closely, you can see that the values of the original raster are passed to the resampled raster. What happens when there is more than one pixel of the original raster coinciding with a single pixel in the new grid? How can we decide which value is going to be passed? In nearest neighbor resampling, the new raster pixels get the value from the *nearest* pixel of the original raster. Note that some of the values may be “lost” this way, since they were not passed on to the new raster^{34}.

### 9.2.3 Bilinear resampling

**Bilinear** resampling is another resampling method. In bilinear resampling, each new raster cell gets a *weighted average* of *four* nearest cells from the input, rather than just one. Bilinear resampling is specified with `method="bilinear"`

in `st_warp`

. Two technical things to note about `st_warp`

when using a method other than `method="near"`

:

`use_gdal=TRUE`

needs to be specified, otherwise the`method`

argument is ignored`no_data_value`

needs to be set to a value outside of the valid range (such as`no_data_value=-9999`

for a DEM raster); this has no effect on the final result, however the specified value is used in a temporary raster which is created when using the`st_warp`

function

For example:

With bilinear resampling, the output raster is “smoothed”, containing new values which are averages of (some of) the values in the original raster (Figure 9.8).

### 9.2.4 Average resampling

Another useful method is the **average** resampling method, where each new cell gets the average of *all* overlapping input cells:

The result of average resampling is shown in Figure 9.9:

In addition to `"near"`

, `"bilinear"`

, and `"average"`

, the `st_warp`

function supports other resampling methods, including: `"cubic"`

, `"cubicspline"`

, `"lanczos"`

, `"mode"`

, `"max"`

, `"min"`

, `"med"`

, `"q1"`

, and `"q3"`

.

Bilinear resampling may be preferred when the result is primarily used for visualization, because the result appears smoother. When using the result for further analysis, however, the original values should be preserved. Generally, either nearest neighbor, mode, or average methods are recommended (Table 9.1), depending on the resolutions of the input and template, and the raster data type. Specifically, keep in mind that when the input raster is categorical, such as a raster with land cover classes `1`

, `2`

, `3`

, etc., nearest neighbor or mode resampling are the only valid methods, while numeric calculations (such as average) on category IDs makes no sense.

Resolution | Type | Resampling method |
---|---|---|

template ≤ input | Categorical / Continuous | Nearest neighbor |

template > input | Categorical | Mode |

template > input | Continuous | Average |

## 9.3 Raster reprojection

Raster **reprojection** is more complex than vector layer reprojection (Section 7.9.2). In addition to transforming (pixel) coordinates, like in vector layer reprojection, raster reprojection requires a resampling step in order to “arrange” the transformed values back into a regular grid (Figure 9.14).

In terms of code, the `st_warp`

function, which we used for resampling (Section 9.2), is used for raster reprojection too. The only difference is that, in raster reprojection, the destination “grid” is in a different CRS than the input raster.

For example, the following expression reprojects the DEM of Haifa from WGS84 (`4326`

) to UTM (`32636`

), using the nearest neighbor resampling method. Note that, in this example, we are not passing a `stars`

object with the destination grid. Instead, we are letting the function to automatically generate the new grid, only specifying the destination CRS (`crs=32636`

) and resolution (`cellsize=90`

):

The original raster (in WGS84) and the reprojected one (in UTM) are shown in Figure 9.10.

Note that the coordinate units in the reprojected raster (Figure 9.10) are no longer degrees, but meters. Also, the area contaning non-missing values is slightly rotated compared to the input, because the WGS84 and UTM systems are not parallel, at least in this particular location. A zoomed-in view (Figure 9.11) of the original raster and the new grid demonstrates the latter more clearly. Compare this situation with the previous examples where the destination grid was in the same CRS (Figure 9.7–9.9), therefore the output raster pixels may be shifted, and/or “stretched” or “shrinked” (due to different origins and resolutions, respectively), however they are always parallel to the input pixels.

As another example, let’s reproject the `MOD13A3_2000_2019.tif`

, from its rarely used Sinusoidal projection, to a projection more widely used and more suitable for the specific region, such as ITM (Table 7.4). First, let’s import the raster from the GeoTIFF file:

Then, we can reproject the raster using `st_warp`

. This time we specify just the destination CRS (`crs=2039`

), letting the function automatically choose an optimal resolution:

The original and reprojected rasters are shown in Figure 9.12:

```
plot(r[,,,1,drop=TRUE], key.pos = 4, axes = TRUE, col = hcl.colors(11, "Spectral"), main = "input (Sinusoidal)")
plot(r_itm[,,,1,drop=TRUE], key.pos = 4, axes = TRUE, col = hcl.colors(11, "Spectral"), main = "output (ITM)")
```

To see the process of reprojection more clearly, let’s examine a small subset of the NDVI raster:

The original and reprojected raster subsets are shown in Figure 9.13.

What happens in the reprojection can be thought of as a two-step process (Figure 9.14). In the first step, the pixel outlines are reprojected as if they were polygons (Section 7.9.2), which results in an *irregular* grid. The grid is then resampled (Section 9.2) to form a regular grid, so that it can be represented by a raster once again.

## 9.4 Focal filters

### 9.4.1 Introduction

So far, we only dealt with arithmetic operations that address the values of each pixel in isolation from neighboring pixels, such as in raster algebra (Sections 6.4 and 6.6.1). Another class of raster operations is where the calculation of each pixels depends on values of *neighboring* pixels.

The most prominent example of a raster calculation based on neighboring cells is *moving window* calculations, also known as **focal filters**. With a moving window calculation, raster values are transformed based on the values from a neighborhood surrounding each pixel. The functions applied on the neighborhood are varied, from simple functions such as `mean`

for a low-pass filter (Section 9.4.2) to more complex functions, such as those used to calculate topographic slope and aspect (Section 9.4.4–9.4.5).

For example, a \(3\times3\) `mean`

filter applied on a raster results in a new raster, where the values are averages of \(3\times3\) neighborhoods centered on that pixel. In Figure 9.15, the highlghted value in the output on the right (`18.7`

) is the average of the highlighted \(3\times3\) neighborhood of the input on the left:

### 9.4.2 Low pass filter

#### 9.4.2.1 What is a low pass filter?

The purpose of the `mean`

, or low-pass, filter (Figure 9.15) is to produce a **smoothed** image, where extreme values (possibly noise) are cacelled out. For example, the result of applying a \(3\times3\) `mean`

filter on a uniform image with one extreme value is shown in Figure 9.16.

How was the value of

`5.4`

obtained, in the result shown in Figure 9.16?

The `stars`

package does not, at the moment, contain a function for focal filtering^{35}. As an exercise, we are going to create our own function for focal filtering. For simplicity, we will restrict ourselves to the case where the focal “window” is \(3\times3\) pixels—which is the most common case.

Our working plan is as follows:

- Write a function named
`get_neighbors`

that accepts a position of the focal cell (row & column) in a`matrix`

, and returns the 9 values of its \(3\times3\) pixel neighborhood - Write another function named
`focal2`

that accepts a two-dimensional`stars`

object and a function`fun`

, iterates over the cells and applies`fun`

on all \(3\times3\) pixel neighborhoods extracted using`get_heighbors`

#### 9.4.2.2 The `get_neighbors`

function

We start with a function that accepts a position (row & column) in a `matrix`

and returns a numeric vector with the \(3\times3\) neighborhood, hereby named `get_neighbors`

. The function accepts a matrix `m`

and a position `pos`

. The `pos`

argument is a vector of length two, of the form `c(row,column)`

. The function extracts the 9 values in the respective \(3\times3\) neighborhood and returns them as a vector of length 9:

```
get_neighbors = function(m, pos) {
i = (pos[1]-1):(pos[1]+1)
j = (pos[2]-1):(pos[2]+1)
as.vector(t(m[i, j]))
}
```

For example, suppose we have a \(5\times5\) matrix `m`

:

```
m = matrix(1:25, ncol = 5, nrow = 5)
m
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 6 11 16 21
## [2,] 2 7 12 17 22
## [3,] 3 8 13 18 23
## [4,] 4 9 14 19 24
## [5,] 5 10 15 20 25
```

Using `get_neighbors`

, we can get the values of just about any \(3\times3\) neighborhood except for the outermost rows and columns (see below). For example, the following expression returns the values of the neighborhood centered on row `3`

, column `4`

:

How does the `get_neighbors`

function work? The function first calculates the required range of rows and columns:

Then, the function extracts the corresponding `matrix`

subset:

transposes it:

and converts to a vector:

Transposing is necessary so that the `matrix`

values are returned *by row*, rather than the default *by column* (Section 5.1.4.1).

Note that our function is not designed to operate on the matrix edges, where the \(3\times3\) neighborhood is incomplete. For example, the following expression produces an error:

#### 9.4.2.3 The `focal2`

function

Now, let’s see how we can use the `get_neighbors`

function to apply a focal filter on a raster. We will use the `dem.tif`

small DEM for demonstration:

First, we create a copy of the raster, named `template`

. The `template`

raster will be used as a “template” when converting the filtered `matrix`

back to a `stars`

object:

Next, we extract the raster values as a `matrix`

:

Note that the function relies on the fact that the `stars`

object has just two dimensions (x and y), in which case `input`

is going to be a `matrix`

. The `matrix`

is transposed, using `t`

, to maintain the right orientation of the values matrix (Figure 6.11). This is important when using functions that distinguish between the north-south and east-west directions, such as *topographic aspect* (Section 9.4.5).

Next, we create *another* `matrix`

to hold the output values. The values are initially set to `NA`

:

Now comes the actual computation. We are using two `for`

loops to go over all raster cells, excluding the first and last rows and columns. For each cell, we:

- extract the \(3\times3\) neighborhood
`[i,j]`

, - apply a function—such as
`mean`

, in this case—on the vector of extracted values, and - place the result into the corresponding cell
`[i,j]`

in the output.

The complete code of the `for`

loops is as follows:

```
for(i in 2:(nrow(input) - 1)) {
for(j in 2:(ncol(input) - 1)) {
v = get_neighbors(input, c(i, j))
output[i, j] = mean(v, na.rm = TRUE)
}
}
```

Note that the function starts at row `i=2`

and ends at row `i=nrow(input)-1`

. Similarly, it starts at column `j=2`

and ends at `j=ncol(input)-1`

.

In the end, when both `for`

loops have been completed, the `output`

matrix contains the new, filtered, raster values. What is left to be done is to put the `matrix`

of new values into the template, to get back a `stars`

object. We are using `t`

once again, to transform the `matrix`

back into the “stars” arrangement (Figure 6.11):

The original image (`x`

) and the filtered image (`template`

) are shown in Figure 9.17. The figure highlights one \(3\times3\) neighborhood in the input, and the corresponding average of that neighborhood in the output^{36}.

```
col = terrain.colors(5)
plot(x, text_values = TRUE, col = col, key.pos = 4, reset = FALSE, main = "input")
plot(st_geometry(st_as_sf(x, na.rm = FALSE)), border = "grey", add = TRUE)
plot(st_geometry(st_as_sf(x[,4:6,6:8])), lwd = 2, add = TRUE)
plot(round(template, 1), text_values = TRUE, col = col, key.pos = 4, reset = FALSE, main = "output")
plot(st_geometry(st_as_sf(template, na.rm = FALSE)), border = "grey", add = TRUE)
plot(st_geometry(st_as_sf(template[,5,7])), lwd = 2, add = TRUE)
```

Why do the outermost rows and columns in Figure 9.17 appear empty? Which value do these pixels contain, and where did it come from?

Let’s manually check the calculation of the \(3\times3\) neighborhood highlighted in Figure 9.17:

Why do you think we got

`NaN`

values in some of the cells? Can you simulate the situation to see how an`NaN`

value is produced? What can we do to get`NA`

, instead, in those pixes that get`NaN`

?

Wrapping up our code in a function, called `focal2`

, can be done as follows. The input raster and the function are replaced with parameters named `r`

and `fun`

, respectively:

```
focal2 = function(r, fun) {
template = r
input = t(template[[1]])
output = matrix(NA, nrow = nrow(input), ncol = ncol(input))
for(i in 2:(nrow(input) - 1)) {
for(j in 2:(ncol(input) - 1)) {
v = get_neighbors(input, c(i, j))
output[i, j] = fun(v)
}
}
template[[1]] = t(output)
return(template)
}
```

How can we pass additional parameters to the function we use, such as `na.rm=TRUE`

for `mean`

? The special dots `...`

argument is used for that. Now, any additional argument(s) passed to `focal2`

(such as `na.rm=TRUE`

) will be passed on as additional argument(s) to `fun`

:

```
focal2 = function(r, fun, ...) {
template = r
input = t(template[[1]])
output = matrix(NA, nrow = nrow(input), ncol = ncol(input))
for(i in 2:(nrow(input) - 1)) {
for(j in 2:(ncol(input) - 1)) {
v = get_neighbors(input, c(i, j))
output[i, j] = fun(v, ...)
}
}
template[[1]] = t(output)
return(template)
}
```

Now that we have a custom focal filter function `focal2`

, let’s try to apply a different filter, such as a *maximum* filter:

```
x_max = focal2(x, max, na.rm = TRUE)
## Warning in fun(v, ...): no non-missing arguments to max; returning -Inf
## Warning in fun(v, ...): no non-missing arguments to max; returning -Inf
```

The reason for the warnings produced by the above expression is that `max`

applied on an empty vector gives `-Inf`

.

```
t(x_max[[1]])
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] NA NA NA NA NA NA NA
## [2,] NA -Inf -Inf 6 9 9 NA
## [3,] NA 61 61 9 10 16 NA
## [4,] NA 132 132 132 11 27 NA
## [5,] NA 254 254 254 146 27 NA
## [6,] NA 254 340 340 340 163 NA
## [7,] NA 254 383 448 448 448 NA
## [8,] NA 253 383 448 448 448 NA
## [9,] NA 253 383 448 448 448 NA
## [10,] NA NA NA NA NA NA NA
```

To get `NA`

instead, we can use a slightly more complex function that first checks if particular neighborhood contains any non-`NA`

values, any only then applies `max`

:

Recall that we used the same principle when applying the `min`

and `max`

functions with `st_apply`

(Section 6.6.1.4).

The resulting raster `x_max`

is shown in Figure 9.18. Indeed, every pixel value in the output raster is the maximal value among the values in its \(3\times3\) neighborhood:

```
col = terrain.colors(10)
plot(x, text_values = TRUE, breaks = "equal", col = col, key.pos = 4, reset = FALSE, main = "input")
plot(st_geometry(st_as_sf(x, na.rm = FALSE)), border = "grey", add = TRUE)
plot(st_geometry(st_as_sf(x[,4:6,6:8])), lwd = 2, add = TRUE)
plot(x_max, breaks = "equal", text_values = TRUE, col = col, key.pos = 4, reset = FALSE, main = "output")
plot(st_geometry(st_as_sf(x_max, na.rm = FALSE)), border = "grey", add = TRUE)
plot(st_geometry(st_as_sf(x_max[,5,7])), lwd = 2, add = TRUE)
```

Let’s try the `focal2`

function on another, bigger raster. For example, we can apply a low pass (i.e., `mean`

) filter on the first layer of the MODIS NDVI raster, as follows:

The original image and the filtered result are shown in Figure 9.19:

```
plot(r_itm1, col = hcl.colors(11, "Spectral"), main = "input")
plot(r_itm1_mean, col = hcl.colors(11, "Spectral"), main = "output")
```

Why are there

`NA`

areas in the raster, even though we used`na.rm=TRUE`

?

### 9.4.3 Maximum filter

For another example, let’s reconstruct the `l_rec`

raster (Section 6.5):

```
l = read_stars("data/landsat_04_10_2000.tif")
red = l[,,,3, drop = TRUE]
nir = l[,,,4, drop = TRUE]
ndvi = (nir - red) / (nir + red)
names(ndvi) = "NDVI"
l_rec = ndvi
l_rec[l_rec < 0.2] = 0
l_rec[l_rec >= 0.2] = 1
```

Given a raster with `0`

and `1`

values, such as `l_rec`

, we may want to convert all `0`

cells neighboring to a `1`

cell to become `1`

. That way, for instance, the areas of the planted forests we see in the center of the image will be come more continuous, which will make it easier to transform those areas into polygons (Section 10.3.2). This can be achieved with a focal filter and the `max`

function:

The original raster and the filtered result are shown in Figure 9.20:

```
plot(l_rec, col = c("grey90", "darkgreen"), main = "input")
plot(l_rec_focal, col = c("grey90", "darkgreen"), main = "output")
```

### 9.4.4 Topographic slope

So far we saw examples where the focal function is a simple built-in one, such as `mean`

or `max`

. In some cases, it is required to use a more complex function. For example, topographic indices such as slope and aspect employ complex functions where direction matters: each of the nine pixels in the \(3\times3\) neignborhood is treated differently.
For example, to calculate topographic slope based on elevation values in a \(3\times3\) neignborhood, the following function can be used. Note that the slope calulation also depends on raster *resolution*, which is passed as an additional parameter named `res`

:

```
slope = function(x, res) {
dzdx = ((x[3] + 2*x[6] + x[9]) - (x[1] + 2*x[4] + x[7])) / (8 * res)
dzdy = ((x[7] + 2*x[8] + x[9]) - (x[1] + 2*x[2] + x[3])) / (8 * res)
atan(sqrt(dzdx^2 + dzdy^2)) * (180 / pi)
}
```

We will not go into details on how the function works. You may refer to the *How slope works* article in the ArcGIS documentation for an explanation.

For example:

The `slope`

function can be passed to `focal2`

to apply the slope calculation on the entire `dem`

raster:

It is also convenient to set raster units with `set_units`

(Section 8.3.2.2). The units of slope are decimal degrees (`°`

):

The resulting topographic slope raster is shown in the left panel in Figure 9.21.

### 9.4.5 Topographic aspect

Another function, as shown below, can be used to calculate topographic aspect:

```
aspect = function(x) {
dzdx = ((x[3] + 2*x[6] + x[9]) - (x[1] + 2*x[4] + x[7])) / 8
dzdy = ((x[7] + 2*x[8] + x[9]) - (x[1] + 2*x[2] + x[3])) / 8
a = (180 / pi) * atan2(dzdy, -dzdx)
if(is.na(a)) return(NA)
if (a < 0) a = 90 - a else
if (a > 90) a = 360 - a + 90 else
a = 90 - a
return(a)
}
```

For details on how the function works, see the *How aspect works* article in the ArcGIS documentation.

For example:

Again, the `aspect`

function can be passed to `focal2`

to apply it on the entire raster:

```
dem_aspect = focal2(dem, aspect)
names(dem_aspect) = "aspect"
dem_aspect[[1]] = set_units(dem_aspect[[1]], "degree")
```

The resulting topographic aspect raster is shown in the right panel in Figure 9.21:

```
plot(dem_slope, breaks = "equal", col = rev(hcl.colors(11, "Spectral")), key.pos = 4)
plot(dem_aspect, breaks = "equal", col = hcl.colors(11, "Spectral"), key.pos = 4)
```

Note that our custom `focal2`

+`get_neighbors`

functions are quite minimal, and can be improved in several ways:

- Being able to set neighbor sizes other than \(3\times3\)
- Dealing with the first/last rows and and columns
- Dealing with rasters that have more than two dimensions (separate filter per dimension? a three-dimensional filter?)
- Making the calculation more efficient (e.g., using C/C++ code inside R or using parallel computation)