Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion .github/workflows/test.yml
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ jobs:
strategy:
matrix:
os: ['ubuntu-latest', 'macos-latest', 'windows-latest']
python: ['3.10', 3.11, 3.12, 3.13, 3.14]
python: [3.12, 3.13, 3.14]
env:
OS: ${{ matrix.os }}
PYTHON: ${{ matrix.python }}
Expand Down
136 changes: 51 additions & 85 deletions xrspatial/zonal.py
Original file line number Diff line number Diff line change
Expand Up @@ -1492,10 +1492,10 @@ def regions(
Parameters
----------
raster : xr.DataArray
connections : int, default=4
neighborhood : int, default=4
4 or 8 pixel-based connectivity.
name: str, default='regions'
output xr.DataArray.name property.
name : str, default='regions'
Output xr.DataArray.name property.

Returns
-------
Expand All @@ -1507,90 +1507,56 @@ def regions(

Examples
--------
.. plot::
:include-source:

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

from xrspatial import generate_terrain
from xrspatial.zonal import regions


# Generate Example Terrain
W = 500
H = 300

template_terrain = xr.DataArray(np.zeros((H, W)))
x_range=(-20e6, 20e6)
y_range=(-20e6, 20e6)

terrain_agg = generate_terrain(
template_terrain, x_range=x_range, y_range=y_range
)

# Edit Attributes
terrain_agg = terrain_agg.assign_attrs(
{
'Description': 'Example Terrain',
'units': 'km',
'Max Elevation': '4000',
}
)

terrain_agg = terrain_agg.rename({'x': 'lon', 'y': 'lat'})
terrain_agg = terrain_agg.rename('Elevation')

# Create Regions
regions_agg = regions(terrain_agg)

# Edit Attributes
regions_agg = regions_agg.assign_attrs({'Description': 'Example Regions',
'units': ''})
regions_agg = regions_agg.rename('Region Value')

# Plot Terrain (Values)
terrain_agg.plot(cmap = 'terrain', aspect = 2, size = 4)
plt.title("Terrain (Values)")
plt.ylabel("latitude")
plt.xlabel("longitude")

# Plot Regions
regions_agg.plot(cmap = 'terrain', aspect = 2, size = 4)
plt.title("Regions")
plt.ylabel("latitude")
plt.xlabel("longitude")

.. sourcecode:: python

>>> print(terrain_agg[200:203, 200:202])
<xarray.DataArray 'Elevation' (lat: 3, lon: 2)>
array([[1264.02296597, 1261.947921 ],
[1285.37105519, 1282.48079719],
[1306.02339636, 1303.4069579 ]])
Coordinates:
* lon (lon) float64 -3.96e+06 -3.88e+06
* lat (lat) float64 6.733e+06 6.867e+06 7e+06
Attributes:
res: (80000.0, 133333.3333333333)
Description: Example Terrain
units: km
Max Elevation: 4000

>>> print(regions_agg[200:203, 200:202])
<xarray.DataArray 'Region Value' (lat: 3, lon: 2)>
array([[39557., 39558.],
[39943., 39944.],
[40327., 40328.]])
Coordinates:
* lon (lon) float64 -3.96e+06 -3.88e+06
* lat (lat) float64 6.733e+06 6.867e+06 7e+06
Attributes:
res: (80000.0, 133333.3333333333)
Description: Example Regions
units:
Max Elevation: 4000
>>> import numpy as np
>>> import xarray as xr
>>> from xrspatial.zonal import regions

>>> # Create a raster with distinct value regions
>>> arr = np.array([[1, 1, 0, 2, 2],
... [1, 1, 0, 2, 2],
... [0, 0, 0, 0, 0],
... [3, 3, 0, 3, 3],
... [3, 3, 0, 3, 3]], dtype=np.float64)
>>> raster = xr.DataArray(arr, dims=['y', 'x'])
>>> print(raster.values)
[[1. 1. 0. 2. 2.]
[1. 1. 0. 2. 2.]
[0. 0. 0. 0. 0.]
[3. 3. 0. 3. 3.]
[3. 3. 0. 3. 3.]]

>>> # With 4-connectivity, each group of connected same-value
>>> # pixels becomes a unique region
>>> result = regions(raster, neighborhood=4)
>>> print(result.values)
[[1. 1. 2. 3. 3.]
[1. 1. 2. 3. 3.]
[2. 2. 2. 2. 2.]
[5. 5. 2. 6. 6.]
[5. 5. 2. 6. 6.]]

>>> # Note: The two bottom-corner 3-regions are separate because
>>> # they are not connected (the zero-valued cross separates them)
>>> print(f"Number of unique regions: {len(np.unique(result.values))}")
Number of unique regions: 5

>>> # With 8-connectivity, diagonal neighbors also connect regions
>>> diagonal = np.array([[1, 0, 1],
... [0, 1, 0],
... [1, 0, 1]], dtype=np.float64)
>>> raster_diag = xr.DataArray(diagonal, dims=['y', 'x'])
>>> result_8 = regions(raster_diag, neighborhood=8)
>>> print(result_8.values)
[[1. 2. 1.]
[2. 1. 2.]
[1. 2. 1.]]

>>> # All 1s are connected diagonally into one region,
>>> # all 0s form another region
>>> print(f"Number of unique regions: {len(np.unique(result_8.values))}")
Number of unique regions: 2
"""
if neighborhood not in (4, 8):
raise ValueError("`neighborhood` value must be either 4 or 8)")
Expand Down