Check Regional Information with Focal Statistics

Recently a young friend said that he encountered a very confusing requirement where one of his clients wanted something other than high-resolution data. The client wants to click on an interactive map, and then some statistics in the region around the clicked point will pop up (presented in tables or charts).

In the beginning, the friend immediately thought this should be a spatial interpolation problem. But it seems wrong after careful thinking. If it is a simple interpolation, where does the statistical value of the nearby area come from? In addition, his team leader told him that the spatial map should also look beautiful. This confuses him even more. How does interpolation to a coarser resolution produce a good-looking spatial map? He asked me what it meant. This requirement sounds quite contradictory.

Introduction

It is really interesting. I told him that putting the needs of his client and the requirements of his team leader together constructed a kind of typical spatial analysis called Focal Statistics.

From the docs (Esri 2020):

Focal statistics perform a neighbourhood operation that computes an output raster, where the value for each output cell is a function of the values of all the input cells that are in a specified neighbourhood around that location. The function performed on the input is a statistic, such as the maximum, average, or sum of all values encountered in that neighbourhood .

To illustrate the neighbourhood processing for Focal Statistics, consider calculating a Sum statistic of the neighbourhood around the processing cell (in blue) with the value of 5 in the following diagram in the input raster. A rectangular 3 by 3 cell neighbourhood shape is specified. The sum of the values of the neighbourhood cells (1 + 2 + 3 + 4 + 6 + 7 + 8 + 9 = 40) plus the value of the processing cell (5) equals 45 (40 + 5 = 45). A value of 45 is given to the cell (in red) in the output raster on the right in the same location as the processing cell in the input raster on the left. The a similar calculation is applied to other cells on the input raster to get the final output raster .

It is worth noting:

• A neighbourhood could be defined with a square (e.g., 3 X 3 in the above illustration image), or rectangle, circle, annulus or an irregular shape. Generally, the neighbourhood is also called kernel.
• When a value in a neighbourhood was defined as NoData, it will be ignored in the statistical calculations. In the neighbourhood , all valid values will be applied;
• The available statistics include maximum, mean, median, minimum, percentile, range, standard deviation, sum, etc. Generally, the default statistics type is mean.
• When the processing cell is near the corners and edges of the input raster, the statistics would not apply all of the cells in the predefined neighbourhood (e.g., 3 X 3). Let's still take a 3 X 3 neighbourhood as an example, some possible cases are displayed the following image.

Simple Practice

Focal Statistics are available in quite a few software such as ArcGIS Pro, QGIS, and SAGA-GIS, etc. There are some implementations are available in R or Python, etc. Let's take the python package of xarray-spatial to have a simple demonstration.

import numpy as np
import xarray as xr
from xrspatial.convolution import circle_kernel
from xrspatial.focal import focal_stats

import cartopy.feature as cfeature
import cartopy.crs as ccrs
import matplotlib.pyplot as plt        

First, let's make a circle kernel (i.e., neighborhood)

it is a 2D array where values of 1 indicate the kernel (valid cell that will participate in the statistical calculation).

kernel = circle_kernel(1, 1, 3)
kernel
Out[5]:array([[0., 0., 0., 1., 0., 0., 0.],
       [0., 1., 1., 1., 1., 1., 0.],
       [0., 1., 1., 1., 1., 1., 0.],
       [1., 1., 1., 1., 1., 1., 1.],
       [0., 1., 1., 1., 1., 1., 0.],
       [0., 1., 1., 1., 1., 1., 0.],
       [0., 0., 0., 1., 0., 0., 0.]])        

Read demo data

da_air = xr.tutorial.open_dataset("air_temperature")["air"][0]        

Focal Statistics - Regional (neighborhood) Mean

da_fs = focal_stats(da_air, kernel, stats_funcs=['mean'])
da_fs        

Visual comparison

fig, (ax1, ax2) = plt.subplots(2, 1, 
                               figsize=(12, 9),
                               subplot_kw={'projection': ccrs.LambertConformal()},
                              )
da_air.plot(ax=ax1, 
            vmin=230, vmax=300, 
            transform=ccrs.PlateCarree(),
            cmap='coolwarm',
            extend='both'
           )
ax1.set_title("Raw Data", fontsize=16)
ax1.coastlines()
ax1.axis('off')


da_fs[0].plot(ax=ax2, vmin=230, vmax=300,             
              transform=ccrs.PlateCarree(),
              cmap='coolwarm',
              extend='both'
           )
ax2.set_title("Focal Statistics - Mean", fontsize=16)
ax2.coastlines()
ax2.axis('off')
_ = ax2.text(1, -0.05,                  
             "Produce by chy@CLIMsystems",
             fontsize=9,
             style='italic',
             horizontalalignment='right',
             transform=ax2.transAxes
            )        

References

ESRI, "The types of operations in Spatial Analyst,"?https://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//005900000017000000

University of Washington, "Introduction to Geographic Information Systems in Forest Resources,"?https://courses.washington.edu/gis250/lessons/raster_analysis1/#focal

https://xarray-spatial.org/reference/_autosummary/xrspatial.focal.focal_stats.html#xrspatial.focal.focal_stats

https://search.r-project.org/CRAN/refmans/raster/html/focal.html

https://pro.arcgis.com/en/pro-app/latest/tool-reference/spatial-analyst/focal-statistics.htm

https://pro.arcgis.com/en/pro-app/latest/tool-reference/spatial-analyst/how-focal-statistics-works.htm

https://saga-gis.sourceforge.io/saga_tool_doc/7.7.1/statistics_grid_1.html