Raster

raster

Attributes

Functions:

get_heights_from_chm(coords, coords_crs, chm_path)

Extract heights from a CHM for the given coordinates

Parameters:

Name	Type	Description	Default
coords	ndarray	Coordinates in the form of an array of shape (n, 2) where n is the number of points	required
chm_path	PATH_TYPE	Path to the CHM raster file	required

Returns:

Type	Description
ndarray	np.ndarray: Heights at the given coordinates

Source code in tree_detection_framework/utils/raster.py

def get_heights_from_chm(
    coords: np.ndarray,
    coords_crs: pyproj.CRS,
    chm_path: PATH_TYPE,
) -> np.ndarray:
    """Extract heights from a CHM for the given coordinates

    Args:
        coords (np.ndarray): Coordinates in the form of an array of shape (n, 2) where n is the number of points
        chm_path (PATH_TYPE): Path to the CHM raster file

    Returns:
        np.ndarray: Heights at the given coordinates
    """
    # Open the CHM raster
    with rasterio.open(chm_path) as src:
        chm_crs = src.crs

        # Check if the coords and CHM have different CRS values
        if coords_crs != chm_crs:
            # Create a GeoSeries for coords so that converting CRS becomes easy
            points_gs = gpd.GeoSeries(
                gpd.points_from_xy(coords[:, 0], coords[:, 1]), crs=coords_crs
            ).to_crs(chm_crs)
            coords_transformed = np.column_stack([points_gs.x, points_gs.y])
        else:
            coords_transformed = coords

        # Sample CHM values at coords
        heights = np.array(
            [
                val[0] if val.size > 0 else np.nan
                for val in src.sample(coords_transformed)
            ]
        )

    return heights

get_valid_raster_region(raster_file)

Compute a single-row geodataframe representing the extent of the valid pixels of a raster

Parameters:

Name	Type	Description	Default
raster_file	PATH_TYPE	Path to the raster	required

Returns:

Type	Description
GeoDataFrame	gpd.GeoDataFrame: valid region

Source code in tree_detection_framework/utils/raster.py

def get_valid_raster_region(raster_file: PATH_TYPE) -> gpd.GeoDataFrame:
    """Compute a single-row geodataframe representing the extent of the valid pixels of a raster

    Args:
        raster_file (PATH_TYPE): Path to the raster

    Returns:
        gpd.GeoDataFrame: valid region
    """
    with rasterio.open(raster_file) as src:
        image = src.read(1)

    # Define the valid mask
    mask = np.logical_and(np.isfinite(image), image != src.nodata)

    # The mask is the important part here, the input "image" is all zeros
    # Get the geometry of each disconnected component. Since we don't care about the value, it is
    # recorded.
    geoms = list(
        {"properties": {}, "geometry": s}
        for (s, v) in rasterio.features.shapes(
            np.zeros_like(mask, dtype=np.uint8),
            mask=mask,
            transform=src.transform,
        )
    )

    # Convert to Geodataframe with only a geometry column
    bounds_gdf = gpd.GeoDataFrame.from_features(geoms, crs=src.crs)
    # Ensure there's only one row
    bounds_gdf = bounds_gdf.dissolve()

    return bounds_gdf

load_downsampled_raster_data(dataset_filename, downsample_factor)

Load a raster file spatially downsampled

Parameters:

Name	Type	Description	Default
dataset	PATH_TYPE	Path to the raster	required
downsample_factor	float	Downsample factor of 10 means that pixels are 10 times larger	required

Returns:

Type	Description
	np.array: The downsampled array in the rasterio (c, h, w) convention
	rasterio.DatasetReader: The reader with the transform updated
	rasterio.Transform: The updated transform

Source code in tree_detection_framework/utils/raster.py

def load_downsampled_raster_data(dataset_filename: PATH_TYPE, downsample_factor: float):
    """Load a raster file spatially downsampled

    Args:
        dataset (PATH_TYPE): Path to the raster
        downsample_factor (float): Downsample factor of 10 means that pixels are 10 times larger

    Returns:
        np.array: The downsampled array in the rasterio (c, h, w) convention
        rasterio.DatasetReader: The reader with the transform updated
        rasterio.Transform: The updated transform
    """
    # Open the dataset handler. Note that this doesn't read into memory.
    dataset = rasterio.open(dataset_filename)

    # resample data to target shape
    data = dataset.read(
        out_shape=(
            dataset.count,
            int(dataset.height / downsample_factor),
            int(dataset.width / downsample_factor),
        ),
        resampling=rasterio.enums.Resampling.bilinear,
    )

    # scale image transform
    updated_transform = dataset.transform * dataset.transform.scale(
        (dataset.width / data.shape[-1]), (dataset.height / data.shape[-2])
    )
    # Return the data and the transform
    return data, dataset, updated_transform

plot_from_dataloader(sample)

Plots an image from the dataset.

Parameters:

Name	Type	Description	Default
sample	dict	A dictionary containing the tile to plot. The 'image' key should have a tensor of shape (C, H, W).	required

Returns:

Type	Description
	matplotlib.figure.Figure: A figure containing the plotted image.

Source code in tree_detection_framework/utils/raster.py

def plot_from_dataloader(sample):
    """
    Plots an image from the dataset.

    Args:
        sample (dict): A dictionary containing the tile to plot. The 'image' key should have a tensor of shape (C, H, W).

    Returns:
        matplotlib.figure.Figure: A figure containing the plotted image.
    """
    # Reorder and rescale the image
    image = sample["image"].permute(1, 2, 0).numpy()

    # Create the figure to plot on and return
    fig, ax = plt.subplots()
    # Plot differently based on the number of channels
    n_channels = image.shape[2]

    # Show with a colorbar and default matplotlib mapping for scalar data
    if n_channels == 1:
        cbar = ax.imshow(image)
        plt.colorbar(cbar, ax=ax)
    # Plot as RGB(A)
    elif n_channels in (3, 4):
        if image.dtype != np.uint8:
            # See if this should be interpreted as data 0-255, even if it's float data
            max_val = np.max(image)
            # If the values are greater than 1, assume it's supposed to be unsigned int8
            # and cast to that so it's properly shown
            if max_val > 1:
                image = image.astype(np.uint8)
        # Plot the image
        ax.imshow(image)
    else:
        raise ValueError(f"Cannot plot image with {n_channels} channels")

    return fig

plot_ortho_chm_overlay(ortho_path, chm_path, chm_alpha=0.5, chm_cmap='viridis', chm_vmin=None, chm_vmax=None, figsize=(10, 10), title='CHM Overlay on Ortho')

Plot a CHM raster overlaid on an ortho image, locked to the ortho's extent.

Parameters:

Name	Type	Description	Default
ortho_path	str	Path to the orthophoto raster (RGB).	required
chm_path	str	Path to the CHM raster.	required
chm_alpha	float	Transparency of the CHM overlay. Defaults to 0.5.	0.5
chm_cmap	str	Colormap for the CHM. Defaults to "viridis".	'viridis'
chm_vmin	Optional[float]	Min value for CHM colormap scaling. Defaults to CHM min.	None
chm_vmax	Optional[float]	Max value for CHM colormap scaling. Defaults to CHM max.	None
figsize	tuple	Figure size. Defaults to (10, 10).	(10, 10)
title	str	Plot title. Defaults to "CHM Overlay on Ortho".	'CHM Overlay on Ortho'

Returns:

Type	Description
Figure	matplotlib Figure object.

Source code in tree_detection_framework/utils/raster.py

def plot_ortho_chm_overlay(
    ortho_path: str,
    chm_path: str,
    chm_alpha: float = 0.5,
    chm_cmap: str = "viridis",
    chm_vmin: Optional[float] = None,
    chm_vmax: Optional[float] = None,
    figsize: tuple = (10, 10),
    title: str = "CHM Overlay on Ortho",
) -> plt.Figure:
    """
    Plot a CHM raster overlaid on an ortho image, locked to the ortho's extent.

    Args:
        ortho_path: Path to the orthophoto raster (RGB).
        chm_path:   Path to the CHM raster.
        chm_alpha:  Transparency of the CHM overlay. Defaults to 0.5.
        chm_cmap:   Colormap for the CHM. Defaults to "viridis".
        chm_vmin:   Min value for CHM colormap scaling. Defaults to CHM min.
        chm_vmax:   Max value for CHM colormap scaling. Defaults to CHM max.
        figsize:    Figure size. Defaults to (10, 10).
        title:      Plot title. Defaults to "CHM Overlay on Ortho".

    Returns:
        matplotlib Figure object.
    """
    with rasterio.open(ortho_path) as ortho_src:
        ortho_img = ortho_src.read([1, 2, 3]).transpose(1, 2, 0)
        ortho_bounds = ortho_src.bounds

    with rasterio.open(chm_path) as chm_src:
        chm_img = chm_src.read(1).astype(float)
        chm_bounds = chm_src.bounds
        nodata = chm_src.nodata

    # Mask nodata values so they render as transparent
    if nodata is not None:
        chm_img[chm_img == nodata] = np.nan

    fig, ax = plt.subplots(figsize=figsize)

    ax.imshow(
        ortho_img,
        extent=[
            ortho_bounds.left,
            ortho_bounds.right,
            ortho_bounds.bottom,
            ortho_bounds.top,
        ],
    )
    chm_plot = ax.imshow(
        chm_img,
        extent=[chm_bounds.left, chm_bounds.right, chm_bounds.bottom, chm_bounds.top],
        alpha=chm_alpha,
        cmap=chm_cmap,
        vmin=chm_vmin,
        vmax=chm_vmax,
    )

    # Lock axes to ortho extent
    ax.set_xlim(ortho_bounds.left, ortho_bounds.right)
    ax.set_ylim(ortho_bounds.bottom, ortho_bounds.top)

    fig.colorbar(chm_plot, ax=ax, label="Canopy Height (m)")
    ax.set_title(title)

    plt.close(fig)  # Prevents Jupyter auto-display

    return fig  # Caller can do display(fig) if needed

reproject_raster(input_file, output_file, dst_crs)

summary

Parameters:

Name	Type	Description	Default
input_file	PATH_TYPE	description	required
output_file	PATH_TYPE	description	required
dst_crs	PATH_TYPE	description	required

Returns:

Name	Type	Description
PATH_TYPE	PATH_TYPE	description

Source code in tree_detection_framework/utils/raster.py

def reproject_raster(
    input_file: PATH_TYPE, output_file: PATH_TYPE, dst_crs: PATH_TYPE
) -> PATH_TYPE:
    """_summary_

    Args:
        input_file (PATH_TYPE): _description_
        output_file (PATH_TYPE): _description_
        dst_crs (PATH_TYPE): _description_

    Returns:
        PATH_TYPE: _description_
    """
    # Taken from here: https://rasterio.readthedocs.io/en/latest/topics/reproject.html
    # Open the source raster
    with rasterio.open(input_file, "r") as src:
        # If it is in the desired CRS, then return the input file since this is a no-op
        if dst_crs == src.crs:
            return input_file

        # Calculate the parameters of the transform for the data in the new CRS
        transform, width, height = calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds
        )
        kwargs = src.meta.copy()
        # Create updated metadata for this new file
        kwargs.update(
            {"crs": dst_crs, "transform": transform, "width": width, "height": height}
        )

        # Open the output file
        with rasterio.open(output_file, "w", **kwargs) as dst:
            # Perform reprojection per band
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=transform,
                    dst_crs=dst_crs,
                    resampling=Resampling.nearest,
                )
    # Return the output file path that the reprojected raster was written to.
    return output_file

show_raster(raster_file_path, downsample_factor=10.0, plt_ax=None, CRS=None)

Show a raster, optionally downsampling or reprojecting it

Parameters:

Name	Type	Description	Default
raster_file_path	PATH_TYPE	Path to the raster file	required
downsample_factor	float	How much to downsample the raster before visualization, this makes it faster and consume less memory	10.0
plt_ax	Optional[axes]	Axes to plot on, otherwise the current ones are used. Defaults to None.	None
CRS	Optional[CRS]	The CRS to reproject the data to if set. Defaults to None.	None

Source code in tree_detection_framework/utils/raster.py

def show_raster(
    raster_file_path: PATH_TYPE,
    downsample_factor: float = 10.0,
    plt_ax: Optional[plt.axes] = None,
    CRS: Optional[pyproj.CRS] = None,
):
    """Show a raster, optionally downsampling or reprojecting it

    Args:
        raster_file_path (PATH_TYPE):
            Path to the raster file
        downsample_factor (float):
            How much to downsample the raster before visualization, this makes it faster and consume
            less memory
        plt_ax (Optional[plt.axes], optional):
            Axes to plot on, otherwise the current ones are used. Defaults to None.
        CRS (Optional[pyproj.CRS], optional):
            The CRS to reproject the data to if set. Defaults to None.
    """
    # Check if the file is georeferenced
    with rasterio.open(raster_file_path) as src:
        if src.crs is None:
            crs_available = False
        else:
            crs_available = True

    # Handle cases where no CRS is available
    if not crs_available and CRS is not None:
        print(f"Warning: No CRS found in the raster. Proceeding without reprojection.")
        CRS = None

    # If the CRS is set, ensure the data matches it
    if CRS is not None:
        # Create a temporary file to write to
        temp_output_filename = tempfile.NamedTemporaryFile(suffix=".tif")
        # Get the name of this file
        temp_name = temp_output_filename.name
        # Reproject the raster. If the CRS was the same as requested, the original raster path will
        # be returned. Otherwise, the reprojected raster will be written to the temp file and that
        # path will be returned.
        raster_file_path = reproject_raster(
            input_file=raster_file_path, output_file=temp_name, dst_crs=CRS
        )
    # Load the downsampled image
    img, _, transform = load_downsampled_raster_data(
        raster_file_path, downsample_factor=downsample_factor
    )
    # Plot the image
    rasterio.plot.show(source=img, transform=transform, ax=plt_ax)