Check Regional Information via Coarsen
We mentioned that regional information could be viewed with focal statistics (such as mean, max, min, median, etc.). However, the young friend returned and told me his client expressed gratitude for what he had done. At the same time, the client pointed out that the output of focal statistics kept the raw high resolutions and added more layers. The data are too big to use. They do not care about using coarse spatial resolutions, but please support his original idea (click on a map and get regional statistical information).
As he already did the focal statistics, he knows the requirements still could not be done with simple interpolation methods. There must be something similar to the focal statistics. Yes, his intuition is correct. The focal statistics apply a sliding shape to derive regional statistical information for each cell on the input raster, while the current requirement could be implemented just with a fixed shape instead of a sliding one. That is to say, the procedure can be described as follows:
It is a Coarsen procedure. Compared with focal statistics, such an operation will reduce the final output size significantly. Moreover, the calculation is more straightforward.
Practice
Next, let's take the Coarsen function from the excellent python package of xarray and have a simple demonstration.
import numpy as np
import xarray as xr
import cartopy.feature as cfeature
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
Read demo data
da_air = xr.tutorial.open_dataset("air_temperature")["air"][0]
Coarsen
Let's split the demo data up into sub-regions with a square shape of 5 X 5 grids using the function of construct.
The boundary="pad" in the followensured that all regions are the same size even though the data does not evenly divide into these sizes.
sub_regions = da_air.coarsen(lat=5, lon=5, boundary="pad").construct(
lon=("x_coarse", "x_fine"), lat=("y_coarse", "y_fine")
)
sub_regions
[out]: xarray.DataArray 'air' y_coarse: 5 y_fine: 5 x_coarse: 11, x_fine: 5
Now 55 new sub-regions have been created, each of size 5 by 5 grids. We can visualize them in the following image, where we can see how they relate to the original data.
sub_regions.plot(
vmin=230, vmax=300,
cmap='coolwarm',
x="x_fine", y="y_fine",
col="x_coarse",
row="y_coarse",
yincrease=False
);
Calculate regional statistics
As we already split the whole region into 55 new sub-regions, now we are now free to easily apply any custom computation to each coarsened sub-region. This would involve specifying that applied functions should act over the "x_fine" and "y_fine" dimensions, but broadcast over the "x_coarse" and "y_coarse" dimensions.
Here, with the function of reduce, more NumPy functions could be applied to sub-regions. For example, the nanpercentile function could be applied to derive the regional medians.
cda_air = sub_regions.reduce(np.nanpercentile, q=50, dim=["x_fine", "y_fine" ])
cda_air['y_coarse'] = sub_regions.lat.mean(dim='y_fine')
cda_air['x_coarse'] = sub_regions.lon.mean(dim='x_fine')
Visual comparison
It could be found that the output with Coarsen much coarser than those with focal statistics. Moreover, it is easier to implement via programming.
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')
cda_air.plot(ax=ax2, vmin=230, vmax=300,
transform=ccrs.PlateCarree(),
cmap='coolwarm',
extend='both'
)
ax2.set_title("Coarsen - Regional Median", 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
)