# Focal#

## Importing Packages#

To get started, we’ll import numpy and xarray-spatial, along with datashader and a set of its functions to help us quickly render images.

```
[1]:
```

```
import numpy as np
import datashader as ds
from datashader.transfer_functions import shade
from datashader.transfer_functions import stack
from datashader.transfer_functions import dynspread
from datashader.transfer_functions import set_background
from datashader.colors import Elevation
import xrspatial
```

## Generate Terrain Data#

The rest of the geo-related functions work with raster data - or rasterized data, i.e. data that’s been aggregated into the regular row-column grid pattern of a raster. In the code below, the datashader Canvas object helps us do this aggregation and returns a raster in the format of an xarray DataArray.

To demonstrate using these raster-based functions, let’s generate some fake terrain as an elevation raster (or *digital elevation model* (dem):

```
[2]:
```

```
from xrspatial import generate_terrain
W = 800
H = 600
cvs = ds.Canvas(
plot_width=W, plot_height=H, x_range=(-20e6, 20e6), y_range=(-20e6, 20e6)
)
terrain = generate_terrain(canvas=cvs)
shade(terrain, cmap=["black", "white"], how="linear")
```

```
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Input In [2], in <cell line: 9>()
4 H = 600
6 cvs = ds.Canvas(
7 plot_width=W, plot_height=H, x_range=(-20e6, 20e6), y_range=(-20e6, 20e6)
8 )
----> 9 terrain = generate_terrain(canvas=cvs)
11 shade(terrain, cmap=["black", "white"], how="linear")
TypeError: generate_terrain() got an unexpected keyword argument 'canvas'
```

The grayscale values above show the elevation, scaled linearly in intensity from dark to bright white (with the large black areas being low elevation). This gives us a glimpse of a fair amount of detail, but we could make it more intuitive by shading it like a landscape.

Datashader’s Set1 maps low values to colors representing water and high ones to colors representing mountaintops, with a range of landscape color in between.

We’ll apply this, but first we’ll also apply xarray-spatial’s hillshade to give an illuminated representation.

```
[3]:
```

```
from xrspatial import hillshade
from datashader.colors import Set1
import pandas as pd
cvs = ds.Canvas(plot_width=W, plot_height=H, x_range=(-20, 20), y_range=(-20, 20))
terrain = generate_terrain(canvas=cvs)
terrain_shaded = shade(terrain, cmap=Elevation, alpha=128, how="linear")
illuminated = hillshade(terrain)
illuminated_shaded = shade(illuminated, cmap=["gray", "white"], alpha=255, how="linear")
stack(illuminated_shaded, terrain_shaded)
```

```
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Input In [3], in <cell line: 7>()
3 import pandas as pd
5 cvs = ds.Canvas(plot_width=W, plot_height=H, x_range=(-20, 20), y_range=(-20, 20))
----> 7 terrain = generate_terrain(canvas=cvs)
8 terrain_shaded = shade(terrain, cmap=Elevation, alpha=128, how="linear")
10 illuminated = hillshade(terrain)
TypeError: generate_terrain() got an unexpected keyword argument 'canvas'
```

## Focal Statistics and Convolutions#

Similar to zonal statistics, focal statistics are used to calculate local statistics, but relative to a *focal* point and a given neighborhood. The neighborhood is a user-defined kernel representing the cells which should be used in calculations. Currently, only circle and annulus kernels are implemented, but any custom kernel can be used given that: - The kernel is a `numpy`

array - The kernel’s dimensions are odd, so a `3x1`

kernel is valid while a `3x2`

is not. This is required for
symmetry around the focal point in the current implementation.

The following example uses focal statistics to calculate the topographic position index (TPI), a measure of local topographic position relative to nearby neighbors. The TPI is scale-dependent and will vary based on the inner and outer radii chosen for the annulus kernel used for this calculation. Once calculated, the TPI can be used to classify slope positions and landforms within landscapes, and can also be used as a numeric feature for model inputs.

First, we use focal’s calc_cellsize to get the scaling right.

Next, we multiply by a larger and smaller integer for the outer and inner annulus radii, respectively.

Then, we generate the kernel with annulus_kernel.

And, finally, we use that kernel with focal.apply to calculate the TPI at each point in the terrain.

To visualize all this, we can shade the TPI and stack it in a composite image with the original terrain.

```
[4]:
```

```
from xrspatial import convolution
from xrspatial import focal
cellsize_x, cellsize_y = convolution.calc_cellsize(terrain)
# Use an annulus kernel with a ring at a distance from 25-30 cells away from focal point
outer_radius = cellsize_x * 30
inner_radius = cellsize_x * 25
kernel = convolution.annulus_kernel(cellsize_x, cellsize_y, outer_radius, inner_radius)
tpi = terrain - focal.apply(terrain, kernel)
tpi_terrain = hillshade(terrain - focal.apply(terrain, kernel))
tpi_terrain_shaded = shade(
tpi_terrain, cmap=["white", "black"], alpha=255, how="linear"
)
stack(illuminated_shaded, tpi_terrain_shaded)
```

```
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Input In [4], in <cell line: 4>()
1 from xrspatial import convolution
2 from xrspatial import focal
----> 4 cellsize_x, cellsize_y = convolution.calc_cellsize(terrain)
5 # Use an annulus kernel with a ring at a distance from 25-30 cells away from focal point
6 outer_radius = cellsize_x * 30
NameError: name 'terrain' is not defined
```

### Convolutions#

The `focal.apply`

function can be computationally expensive depending on the size of the kernal and image. Additionally, we’d also like to extend focal’s ability to using custom convolution kernels on the raster.

For example, we may want to apply some kernels from image processing. Below, we’ll use one of them, called the Sobel operator, which is a crude, but computationally inexpensive way to do edge-detection. For the horizontal operator, an approximation of the derivative in the horizontal dimension is calculated, so we’ll set up a kernel below to do that.

(Note: By default, the `convolution`

module will use the local CUDA-enabled GPU unit if it is available. To use only the CPU, you can pass `use_cuda=False`

to `convolution.convolve_2d`

.)

First, we’ll set up a horizontal Sobel kernel manually as a 2D numpy array. Notice the vertical column of zeros and the symmetry on either side of this column.

Next, we’ll apply this kernel to the terrain values to generate a corresponding array of Sobel values.

Finally, we’ll set those values into a DataArray raster with the proper coordinates, dimensions, and attributes from the terrain.

To visualize all this, we’ll apply hillshade and shade to the sobel terrain and stack it in an image with the original.

Notice all the emphasized horizontal edges. (You can also edit the kernel numbers below to make it vertical and see the vertical edges emphasized, instead.)

```
[5]:
```

```
from xrspatial import convolution
from xarray import DataArray
# Use Sobel operator
kernel = np.array([[1, 0, -1], [2, 0, -2], [1, 0, -1]])
display(kernel)
sobel_values = convolution.convolve_2d(terrain.values, kernel)
sobel = DataArray(
sobel_values,
coords=terrain.coords,
dims=terrain.dims,
attrs=terrain.attrs,
)
sobel_terrain = hillshade(sobel)
sobel_terrain_shaded = shade(
sobel_terrain, cmap=["white", "black"], alpha=255, how="linear"
)
stack(illuminated_shaded, sobel_terrain_shaded)
```

```
array([[ 1, 0, -1],
[ 2, 0, -2],
[ 1, 0, -1]])
```

```
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Input In [5], in <cell line: 7>()
5 kernel = np.array([[1, 0, -1], [2, 0, -2], [1, 0, -1]])
6 display(kernel)
----> 7 sobel_values = convolution.convolve_2d(terrain.values, kernel)
8 sobel = DataArray(
9 sobel_values,
10 coords=terrain.coords,
11 dims=terrain.dims,
12 attrs=terrain.attrs,
13 )
14 sobel_terrain = hillshade(sobel)
NameError: name 'terrain' is not defined
```