Skip to content

Raster data utilities

combine_raster_bands(input_rasters)

Combine multiple rasters into one multiband raster.

The input rasters can be either singleband or multiband. All bands are stacked in the order they are extracted from the input raster list.

All input rasters must have matching spatial metadata (extent, pixel size, CRS).

Parameters:

Name Type Description Default
input_rasters Sequence[DatasetReader]

List of rasters to combine.

 required

Returns:

Type Description
ndarray

The combined raster data.
Profile

The updated raster profile.

Raises:

Type Description
InvalidParameterValueException

Input rasters contains only one raster.
NonMatchingRasterMetadataException

Input rasters have mismatching spatial metadata.
Source code in eis_toolkit/utilities/raster.py 
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
@beartype
def combine_raster_bands(input_rasters: Sequence[rasterio.io.DatasetReader]) -> Tuple[np.ndarray, profiles.Profile]:
    """
    Combine multiple rasters into one multiband raster.

    The input rasters can be either singleband or multiband. All bands are stacked in the order they are
    extracted from the input raster list.

    All input rasters must have matching spatial metadata (extent, pixel size, CRS).

    Args:
        input_rasters: List of rasters to combine.

    Returns:
        The combined raster data.
        The updated raster profile.

    Raises:
        InvalidParameterValueException: Input rasters contains only one raster.
        NonMatchingRasterMetadataException: Input rasters have mismatching spatial metadata.
    """
    if len(input_rasters) < 2:
        raise InvalidParameterValueException(
            f"Expected multiple rasters in the input_rasters list. Rasters: {len(input_rasters)}."
        )

    profiles = []
    bands_arrays = []
    for raster in input_rasters:
        profiles.append(raster.profile)
        for i in range(1, raster.count + 1):
            bands_arrays.append(raster.read(i))

    if not check_raster_grids(profiles, same_extent=True):
        raise NonMatchingRasterMetadataException("Input rasters have mismatching metadata/profiles.")

    out_image = np.stack(bands_arrays, axis=0)

    out_profile = profiles[0].copy()
    out_profile["count"] = len(bands_arrays)

    return out_image, out_profile

profile_from_extent_and_pixel_size(extent, pixel_size, round_strategy='up')

Create a raster profile from given extent and pixel size.

If extent and pixel size do not match exactly, i.e. raster width and height calcalated from bounds and pixel size are not integers, rounding for the width and height is performed.

Parameters:

Name Type Description Default
extent Tuple[Number, Number, Number, Number]

Raster extent in the form (coord_west, coord_east, coord_south, coord_north).

 required
pixel_size Union[Number, Tuple[Number, Number]]

Desired pixel size. If two values are provided, first is used for x and second for y. If one value is provided, the value is used for both directions.

 required
round_strategy Literal[nearest, up, down]

The rounding strategy if extent and pixel size do not match exactly. Defaults to "up".

 'up'

Returns:

Type Description
Profile

Rasterio profile.
Source code in eis_toolkit/utilities/raster.py 
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
@beartype
def profile_from_extent_and_pixel_size(
    extent: Tuple[Number, Number, Number, Number],
    pixel_size: Union[Number, Tuple[Number, Number]],
    round_strategy: Literal["nearest", "up", "down"] = "up",
) -> profiles.Profile:
    """
    Create a raster profile from given extent and pixel size.

    If extent and pixel size do not match exactly, i.e. raster width and height
    calcalated from bounds and pixel size are not integers, rounding for the width and
    height is performed.

    Args:
        extent: Raster extent in the form (coord_west, coord_east, coord_south, coord_north).
        pixel_size: Desired pixel size. If two values are provided, first is used for x and second for y.
            If one value is provided, the value is used for both directions.
        round_strategy: The rounding strategy if extent and pixel size do not match exactly.
            Defaults to "up".

    Returns:
        Rasterio profile.
    """
    if isinstance(pixel_size, Tuple):
        pixel_size_x, pixel_size_y = pixel_size[0], pixel_size[1]
    else:
        pixel_size_x, pixel_size_y = pixel_size, pixel_size

    coord_west, coord_east, coord_south, coord_north = extent
    width_raw = abs(coord_east - coord_west) / pixel_size_x
    height_raw = abs(coord_north - coord_south) / pixel_size_y
    if round_strategy == "down":
        width, height = floor(width_raw), floor(height_raw)
    elif round_strategy == "up":
        width, height = ceil(width_raw), ceil(height_raw)
    else:
        width, height = round(width_raw), round(height_raw)

    raster_meta = {
        "transform": transform.from_bounds(
            coord_west, coord_south, coord_east, coord_north, width=width, height=height
        ),
        "width": width,
        "height": height,
    }
    raster_profile = profiles.Profile(raster_meta)
    return raster_profile

split_raster_bands(raster)

Split multiband raster into singleband rasters.

Parameters:

Name Type Description Default
raster DatasetReader

Input multiband raster.

 required

Returns:

Type Description
Sequence[Tuple[ndarray, Profile]]

Output singleband raster list. List elements are tuples where first element is raster data (2D) and second element is raster profile.

Raises:

Type Description
InvalidParameterValueException

Input raster contains only one band.
Source code in eis_toolkit/utilities/raster.py 
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
@beartype
def split_raster_bands(raster: rasterio.io.DatasetReader) -> Sequence[Tuple[np.ndarray, profiles.Profile]]:
    """
    Split multiband raster into singleband rasters.

    Args:
        raster: Input multiband raster.

    Returns:
        Output singleband raster list. List elements are tuples where first element is raster data (2D) and second \
        element is raster profile.

    Raises:
        InvalidParameterValueException: Input raster contains only one band.
    """
    out_rasters = []
    count = raster.meta["count"]
    if count < 2:
        raise InvalidParameterValueException(f"Expected multiple bands in the input raster. Bands: {count}.")

    for i in range(1, count + 1):
        band_data = raster.read(i)
        band_profile = raster.profile.copy()
        band_profile.update({"count": 1, "dtype": band_data.dtype})
        out_rasters.append((band_data, band_profile))
    return out_rasters

stack_raster_arrays(arrays)

Stack 2D and 3D NumPy arrays (each representing a raster with one or multiple bands) along the bands axis.

Parameters:

Name Type Description Default
arrays Sequence[ndarray]

List of 2D and 3D NumPy arrays. Each 2D array should have shape (height, width). and 3D array shape (bands, height, width).

 required

Returns:

Type Description
ndarray

A single 3D NumPy array where the first dimension size equals the total number of bands.

Raises:

Type Description
InvalidDataShapeException

Input raster arrays have mismatching shapes or all input rasters are not 2D or 3D.
Source code in eis_toolkit/utilities/raster.py 
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
@beartype
def stack_raster_arrays(arrays: Sequence[np.ndarray]) -> np.ndarray:
    """
    Stack 2D and 3D NumPy arrays (each representing a raster with one or multiple bands) along the bands axis.

    Parameters:
        arrays: List of 2D and 3D NumPy arrays. Each 2D array should have shape (height, width).
            and 3D array shape (bands, height, width).

    Returns:
        A single 3D NumPy array where the first dimension size equals the total number of bands.

    Raises:
        InvalidDataShapeException: Input raster arrays have mismatching shapes or all input rasters are not 2D or 3D.
    """
    processed_arrays = []
    for array in arrays:
        # Add a new axis if the array is 2D
        if array.ndim == 2:
            array = array[np.newaxis, :]
        elif array.ndim != 3:
            raise InvalidDataShapeException("All raster arrays must be 2D or 3D for stacking.")
        processed_arrays.append(array)

    shape_set = {arr.shape[1:] for arr in processed_arrays}
    if len(shape_set) != 1:
        raise InvalidDataShapeException(
            "All raster arrays must have the same shape in 2 last dimensions (height, width)."
        )

    # Stack along the first axis
    stacked_array = np.concatenate(processed_arrays, axis=0)

    return stacked_array