None

Note

This tutorial was generated from an IPython notebook that can be downloaded here.

Creating custom focal statistic function

import focal_stats
import rasterio as rio
import matplotlib.pyplot as plt
import numpy as np
import os

os.chdir("../../../")

Loading raster (containing water table depth (Fan et al., 2017)).

with rio.open("data/wtd.tif") as f:
    a = f.read(1).astype(np.float64)
    a[a == -999.9] = np.nan

Inspecting the data

plt.imshow(a, cmap='Blues', vmax=100)
plt.title("Water table depth")
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x7f3f63ec3e90>
../_images/custom_focal_stats_6_1.png

Creating custom focal mean

Firstly, a windowed version of the input raster needs to be defined.

a_windowed = focal_stats.rolling_window(a, window_size=5)
print(a.shape, a_windowed.shape)

(4320, 5400) (4316, 5396, 5, 5)

This windowed version has a slightly different shape on the first two axes. This is because there is no window behaviour defined on the edges. If this is undesired the original array can be padded with the missing number of columns and rows with numpy.pad. Through this function many different edge value assumptions can be made. Here I use the example of continuing with the closest values.

a_padded = np.pad(a, pad_width=2, mode='edge')
a_windowed_padded = focal_stats.rolling_window(a_padded, window_size=5)
print(a.shape, a_windowed_padded.shape)

(4320, 5400) (4320, 5400, 5, 5)

This has the result that the input and output raster share their first two axes.

Now the only thing that needs to happen is a mean operation on the third and fourth axes:

a_mean = a_windowed.mean(axis=(2, 3))

Plotting this shows that the operation generates an image that is very close to the original raster, with some limited smoothing

plt.imshow(a_mean, cmap="Blues", vmax=100)
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x7f3f62503510>
../_images/custom_focal_stats_17_1.png

This can be captured in a custom focal_mean function as follows:

def focal_mean(a, window_size):
    a_windowed = focal_stats.rolling_window(a, window_size=window_size)
    return a_windowed.mean(axis=(2, 3))

Resulting in the same image:

plt.imshow(focal_mean(a, window_size=5), cmap="Blues", vmax=100)
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x7f3f623f4550>
../_images/custom_focal_stats_21_1.png

Note that if a single NaN-value was present in the window, it results in a NaN-value. I dealt with this by inserting 0 in the pixels with NaN-values and using the sum of this array divided by the number of valid values per window (e.g. rolling_sum(~np.isnan(a), window_size=5)).