raster

Functions for spatial processing of raster TIFF files.

default_summary_func(v)

Sum an array and round to the nearest integer.

Args:

Returns:

Source code in rastertoolkit/raster.py

def default_summary_func(
    v: np.ndarray,
) -> int:
    """
    Sum an array and round to the nearest integer.

    Args:

    Returns:
    """
    return int(np.round(np.sum(v), 0))

extract_xy_info_from_raster(raster)

Extracts x, y, dx, and dy from a TiffPage object.

Parameters:

Name	Type	Description	Default
raster	TiffPage	Single tiff layer.	required

Returns:

Name	Type	Description
tuple	tuple[float, float, float, float]	A tuple of x, y, dx, and dy.

Source code in rastertoolkit/raster.py

def extract_xy_info_from_raster(
    raster: TiffPage,
) -> tuple[float, float, float, float]:
    """
    Extracts x, y, dx, and dy from a TiffPage object.

    Args:
        raster (TiffPage): Single tiff layer.

    Returns:
        tuple: A tuple of x, y, dx, and dy.
    """

    # Extract data from raster
    tags = get_tiff_tags(raster)
    if "ModelTiepointTag" in tags:
        point = tags["ModelTiepointTag"]
        scale = tags["ModelPixelScaleTag"]
        x0, y0 = point[3], point[4]
        dx, dy = scale[0], -scale[1]
    elif "ModelTransformationTag" in tags:
        vector = tags["ModelTransformationTag"]
        x0, y0 = vector[3], vector[7]
        dx, dy = vector[0], vector[5]
    else:
        raise ValueError('Invalid GeoTIFF tags.')

    # Make sure values are in range
    if not (-180 < x0 < 180):
        raise ValueError(f"Tie point x coordinate (longitude) has invalid range: {x0}.")
    if not (-85 < y0 < 85):
        raise ValueError(f"Tie point y coordinate (latitude) has invalid range: {y0}.")
    if not (0 < dx < 1):
        raise ValueError(f"Pixel dx scale has invalid range: {dx}.")
    if not (-1 < dy < 0):
        raise ValueError(f"Pixel dy scale has invalid range: {dy}.")
    nodata = getattr(raster, "nodata", None)
    if nodata is not None and nodata >= 0:
        raise ValueError(f"Invalid no-data attribute; must be negative: {nodata}.")

    return x0, y0, dx, dy

get_tiff_tags(raster)

Reads tags from a TiffPage object.

Parameters:

Name	Type	Description	Default
raster	TiffPage	Single tiff layer.	required

Returns:

Name	Type	Description
dict	dict[str, Any]	A dictionary of tag names and values.

Source code in rastertoolkit/raster.py

def get_tiff_tags(
    raster: TiffPage,
) -> dict[str, Any]:
    """
    Reads tags from a TiffPage object.

    Args:
        raster (TiffPage): Single tiff layer.

    Returns:
        dict: A dictionary of tag names and values.
    """
    return {tag_obj.name: tag_obj.value for tag_obj in raster.tags}

init_sparse_matrix(raster, band)

Initialize a matrix from a raster TiffPage object with values > 0

Parameters:

Name	Type	Description	Default
raster	TiffPage	Single tiff page.	required

Returns:

Source code in rastertoolkit/raster.py

def init_sparse_matrix(
    raster: TiffPage,
    band: int,
) -> np.ndarray:
    """
    Initialize a matrix from a raster TiffPage object with values > 0

    Args:
        raster (TiffPage): Single tiff page.

    Returns:
    """

    # Extract data from raster
    x0, y0, dx, dy = extract_xy_info_from_raster(raster)

    dat_mat = raster.asarray()
    if (dat_mat.ndim == 2 and band > 0):
        raise ValueError('Invalid raster band.')
    if (dat_mat.ndim == 3):
        if (dat_mat.shape[-1] <= band):
            raise ValueError('Invalid raster band.')
        dat_mat = dat_mat[:, :, band]

    xy_ints = np.argwhere(dat_mat[:, :] > 0)
    sparse_data = np.zeros((xy_ints.shape[0], 3), dtype=float)

    # Construct sparse matrix of (long, lat, data)
    sparse_data[:, 0] = x0 + dx * xy_ints[:, 1] + dx / 2.0
    sparse_data[:, 1] = y0 + dy * xy_ints[:, 0] + dy / 2.0
    sparse_data[:, 2] = dat_mat[xy_ints[:, 0], xy_ints[:, 1]]

    return sparse_data

interpolate_at_weight_data(shape, weight_clip, value_clip, data_bool)

Interpolate at weight data.

Parameters:

Name	Type	Description	Default
shape	ShapeView	Shape object.	required
weight_clip	ndarray	Clipped weight data.	required
value_clip	ndarray	Clipped value data.	required
data_bool	bool	Boolean indicating if the data is interior to the shape.	required

Returns:

Name	Type	Description
float	float	The interpolated value at weight data.

Source code in rastertoolkit/raster.py

def interpolate_at_weight_data(
    shape: ShapeView,
    weight_clip: np.ndarray,
    value_clip: np.ndarray,
    data_bool: bool,
) -> float:
    """
    Interpolate at weight data.

    Args:
        shape (ShapeView): Shape object.
        weight_clip (np.ndarray): Clipped weight data.
        value_clip (np.ndarray): Clipped value data.
        data_bool (bool): Boolean indicating if the data is interior to the shape.

    Returns:
        float: The interpolated value at weight data.
    """
    # Calculate total weight
    weight = np.sum(weight_clip[data_bool, 2])

    # Prep interpolate coordinates and value arguments
    value_args = [value_clip[:, 0:2], value_clip[:, 2]]

    if weight > 0:
        # Interpolate at weight, assign -1 for problems
        val_est = interpolate.griddata(*value_args, weight_clip[:, 0:2], fill_value=-1)
        if -1 in val_est:
            err_dex = (val_est == -1)
            # Use the nearest value for problems
            val_rev = interpolate.griddata(
                *value_args, weight_clip[err_dex, 0:2], method="nearest"
            )
            val_est[err_dex] = val_rev
        # Use weight values
        final_val = np.sum(weight_clip[data_bool, 2] * val_est[data_bool]) / weight
    else:
        # No weight data inside shape, interpolate at shape boundary, assign -1 for problems
        val_est = interpolate.griddata(*value_args, shape.points[:, 0:2], fill_value=-1)
        if -1 in val_est:
            err_dex = val_est == -1
            # Use the nearest value for problems
            val_rev = interpolate.griddata(
                *value_args, shape.points[err_dex, 0:2], method="nearest"
            )
            val_est[err_dex] = val_rev

        # Average values at shape perimeter
        final_val = np.mean(val_est)

    return final_val

is_interior(shape, data_clip)

Check if the data is interior to the shape.

Parameters:

Name	Type	Description	Default
shape	ShapeView	Shape object.	required
data_clip	ndarray	Clipped data.	required

Returns:

Name	Type	Description
bool	bool	True if the data is interior to the shape.

Source code in rastertoolkit/raster.py

def is_interior(
    shape: ShapeView,
    data_clip: np.ndarray,
) -> bool:
    """
    Check if the data is interior to the shape.

    Args:
        shape (ShapeView): Shape object.
        data_clip (np.ndarray): Clipped data.

    Returns:
        bool: True if the data is interior to the shape.
    """
    # Track booleans (indicates if lat/long is interior)
    data_bool = np.zeros(data_clip.shape[0], dtype=bool)

    # Iterate over parts of shapefile
    for path_shp, area_prt in zip(shape.paths, shape.areas):
        # Union of positive areas; intersection with negative areas
        if area_prt > 0:
            data_bool = np.logical_or(
                data_bool, path_shp.contains_points(data_clip[:, :2])
            )
        else:
            data_bool = np.logical_and(
                data_bool, np.logical_not(path_shp.contains_points(data_clip[:, :2]))
            )

    return data_bool

print_status(shape, data_dict, k1, shape_count)

Print status message.

Args:

Returns:

Source code in rastertoolkit/raster.py

def print_status(
    shape: ShapeView,
    data_dict: dict,
    k1: int,
    shape_count: int,
) -> None:
    """
    Print status message.

    Args:

    Returns:
    """
    perc = round(100 * (k1 + 1) / shape_count)
    print(
        k1 + 1,
        "of",
        shape_count,
        f"({perc}%)",
        shape.name,
        shape.center,
        data_dict[shape.name],
    )

raster_clip(raster_file, shape_stem, raster_page=0, raster_band=0, shape_attr='DOTNAME', attr_filter=None, summary_func=None, include_latlon=False, quiet=False)

Extracts data from a raster based on shapes.

Parameters:

Name	Type	Description	Default
raster_file	str	Local path to a raster file.	required
shape_stem	str	Local path stem referencing a set of shape files.	required
raster_page	int	Index of the page to extract from the raster.	0
raster_band	int	Index of the band to extract from the raster.	0
shape_attr	str	The shape attribute name to be used as the output dictionary key. This attribute is used as the shape name.	'DOTNAME'
attr_filter	str	String that the shape name must start with to be included. Excluded shapes are ignored.	None
summary_func	Callable	Aggregation function to be used for summarizing clipped data for each shape.	None
include_latlon	bool	Flag to include lat/lon in the dictionary entry.	False
quiet	bool	Flag to control whether status messages are printed.	False

Returns:

Name	Type	Description
dict	dict[str, Union[float, int]]	A dictionary with dot names as keys and calculated aggregations as values.

Source code in rastertoolkit/raster.py

def raster_clip(
    raster_file: Union[str, Path],
    shape_stem: Union[str, Path],
    raster_page: int = 0,
    raster_band: int = 0,
    shape_attr: str = "DOTNAME",
    attr_filter: Union[str, None] = None,
    summary_func: Callable = None,
    include_latlon: bool = False,
    quiet: bool = False,
) -> dict[str, Union[float, int]]:
    """
    Extracts data from a raster based on shapes.

    Args:
        raster_file (str): Local path to a raster file.
        shape_stem (str): Local path stem referencing a set of shape files.
        raster_page (int): Index of the page to extract from the raster.
        raster_band (int): Index of the band to extract from the raster.
        shape_attr (str): The shape attribute name to be used as the output dictionary key.
            This attribute is used as the shape name.
        attr_filter (str): String that the shape name must start with to be included.
            Excluded shapes are ignored.
        summary_func (Callable): Aggregation function to be used for summarizing clipped
            data for each shape.
        include_latlon (bool, optional): Flag to include lat/lon in the dictionary entry.
        quiet (bool, optional): Flag to control whether status messages are printed.

    Returns:
        dict: A dictionary with dot names as keys and calculated aggregations as values.
    """
    if not Path(raster_file).is_file():
        raise FileNotFoundError(f"Raster file not found: {raster_file}")

    print("Loading data...")

    # Load data, init sparse matrix
    shapes = ShapeView.from_file(shape_stem, shape_attr, attr_filter)
    raster = TiffFile(raster_file)
    sparse_data = init_sparse_matrix(raster.pages[raster_page], raster_band)

    # Output dictionary
    data_dict = dict()
    shape_len = len(shapes)
    print("Clipping:")

    fts = {}
    # Init the futures executor
    executor = ThreadPoolExecutor(max_workers=(os.cpu_count() - 1))

    # Iterate over shapes in shapefile
    for k1, shp in enumerate(shapes):
        fts[k1] = executor.submit(
            raster_clip_single,
            shp=shp,
            sparse_data=sparse_data,
            k1=k1,
            shape_len=shape_len,
            summary_func=summary_func,
            include_latlon=include_latlon,
            quiet=quiet,
        )

    for k1, ft in fts.items():
        data_dict.update(ft.result())

    executor.shutdown(wait=True)

    return data_dict

raster_clip_single(shp, sparse_data, k1, shape_len, summary_func, include_latlon, quiet)

Extracts data from a raster based on shapes.

Parameters:

Name	Type	Description	Default
shp	ShapeView	Shape object.	required
sparse_data	ndarray	Sparse matrix of raster data.	required
k1	int	Index of the shape.	required
shape_len	int	Total number of shapes.	required
summary_func	Callable	Aggregation function to be used for summarizing clipped data for each shape.	required
include_latlon	bool	Flag to include lat/lon in the dictionary entry.	required
quiet	bool	Flag to control whether status messages are printed.	required

Returns:

Name	Type	Description
dict	dict[str, Union[float, int]]	A dictionary with dot names as keys and calculated aggregations as values.

Source code in rastertoolkit/raster.py

def raster_clip_single(
    shp: ShapeView,
    sparse_data: np.ndarray,
    k1: int,
    shape_len: int,
    summary_func: Callable,
    include_latlon: bool,
    quiet: bool,
) -> dict[str, Union[float, int]]:
    """
    Extracts data from a raster based on shapes.

    Args:
        shp (ShapeView): Shape object.
        sparse_data (np.ndarray): Sparse matrix of raster data.
        k1 (int): Index of the shape.
        shape_len (int): Total number of shapes.
        summary_func (Callable): Aggregation function to be used for summarizing clipped
            data for each shape.
        include_latlon (bool): Flag to include lat/lon in the dictionary entry.
        quiet (bool): Flag to control whether status messages are printed.

    Returns:
        dict: A dictionary with dot names as keys and calculated aggregations as values.
    """
    data_dict = {}
    show_status = not quiet or k1 % 1000 == 0 or k1 in [0, shape_len - 1]
    # Null shape; error in shapefile
    shp.validate()

    # Subset matrix for clipping
    data_clip = subset_matrix_for_clipping(shp, sparse_data)

    if data_clip.shape[0] == 0:
        data_dict[shp.name] = summary_entry(shp, {"pop": 0}, include_latlon)
        if show_status:
            print_status(shp, data_dict, k1, shape_len)
        return data_dict

    # Pop values
    value = data_clip[is_interior(shp, data_clip), 2]

    # Entry dictionary
    summary_func = summary_func or default_summary_func
    entry = {"pop": summary_func(value)}

    # Set entry and print status
    data_dict[shp.name] = summary_entry(shp, entry, include_latlon)
    if show_status:
        print_status(shp, data_dict, k1, shape_len)

    return data_dict

raster_clip_weighted(raster_weight, raster_value, shape_stem, raster_weight_page=0, raster_weight_band=0, raster_value_page=0, raster_value_band=0, shape_attr='DOTNAME', attr_filter=None, weight_summary_func=None, include_latlon=False)

Extracts data from a raster based on shapes.

Parameters:

Name	Type	Description	Default
raster_weight	str	Local path to a raster file used for weights.	required
raster_value	str	Local path to a raster file used for values.	required
shape_stem	str	Local path stem referencing a set of shape files.	required
raster_weight_page	int	Index of the page to extract from the weights raster.	0
raster_weight_band	int	Index of the band to extract from the weights raster.	0
raster_value_page	int	Index of the page to extract from the value raster.	0
raster_value_band	int	Index of the band to extract from the value raster.	0
shape_attr	str	The shape attribute name to be used as the output dictionary key. This attribute is used as the shape name.	'DOTNAME'
attr_filter	str	String that the shape name must start with to be included. Excluded shapes are ignored.	None
weight_summary_func	Callable	Aggregation function to be used for summarizing clipped data for each shape.	None
include_latlon	bool	Flag to include lat/lon in the dictionary entry.	False

Returns:

Name	Type	Description
dict	dict[str, Union[float, int]]	A dictionary with dot names as keys and calculated aggregations as values.

Source code in rastertoolkit/raster.py

def raster_clip_weighted(
    raster_weight: Union[str, Path],
    raster_value: Union[str, Path],
    shape_stem: Union[str, Path],
    raster_weight_page: int = 0,
    raster_weight_band: int = 0,
    raster_value_page: int = 0,
    raster_value_band: int = 0,
    shape_attr: str = "DOTNAME",
    attr_filter: Union[str, None] = None,
    weight_summary_func: Callable = None,
    include_latlon: bool = False,
) -> dict[str, Union[float, int]]:
    """
    Extracts data from a raster based on shapes.

    Args:
        raster_weight (str): Local path to a raster file used for weights.
        raster_value (str): Local path to a raster file used for values.
        shape_stem (str): Local path stem referencing a set of shape files.
        raster_weight_page (int): Index of the page to extract from the weights raster.
        raster_weight_band (int): Index of the band to extract from the weights raster.
        raster_value_page (int): Index of the page to extract from the value raster.
        raster_value_band (int): Index of the band to extract from the value raster.
        shape_attr (str): The shape attribute name to be used as the output dictionary key.
            This attribute is used as the shape name.
        attr_filter (str): String that the shape name must start with to be included.
            Excluded shapes are ignored.
        weight_summary_func (Callable): Aggregation function to be used for summarizing
            clipped data for each shape.
        include_latlon (bool, optional): Flag to include lat/lon in the dictionary entry.

    Returns:
        dict: A dictionary with dot names as keys and calculated aggregations as values.
    """
    if not Path(raster_weight).is_file():
        raise FileNotFoundError(f"Population raster file not found: {raster_weight}")
    if not Path(raster_value).is_file():
        raise FileNotFoundError(f"Values raster file not found: {raster_value}")

    print("Loading data...")

    # Load data shape and rasters
    shapes = ShapeView.from_file(shape_stem, shape_attr, attr_filter)
    raster_weights = TiffFile(raster_weight)
    raster_values = TiffFile(raster_value)

    # Init sparse matrices
    sparse_weight = init_sparse_matrix(raster_weights.pages[raster_weight_page], raster_weight_band)
    sparse_val = init_sparse_matrix(raster_values.pages[raster_value_page], raster_value_band)

    # Output dictionary
    data_dict = {}
    shape_len = len(shapes)
    print("Clipping:")

    fts = {}
    # Init the futures executor
    executor = ThreadPoolExecutor(max_workers=(os.cpu_count() - 1))

    # Iterate over shapes in shapefile
    for k1, shp in enumerate(shapes):
        fts[k1] = executor.submit(
            raster_clip_weighted_single,
            shp=shp,
            sparse_weight=sparse_weight,
            sparse_val=sparse_val,
            k1=k1,
            shape_len=shape_len,
            weight_summary_func=weight_summary_func,
            include_latlon=include_latlon,
        )

    for k1, ft in fts.items():
        data_dict.update(ft.result())

    executor.shutdown(wait=True)

    return data_dict

raster_clip_weighted_single(shp, sparse_weight, sparse_val, k1, shape_len, weight_summary_func, include_latlon)

Extracts weighted data from a raster based on a single shape.

Parameters:

Name	Type	Description	Default
shp	ShapeView	Shape object.	required
sparse_weight	ndarray	Sparse matrix of weight (population) raster data.	required
sparse_val	ndarray	Sparse matrix of value raster data.	required
k1	int	Index of the shape.	required
shape_len	int	Total number of shapes.	required
weight_summary_func	Callable	Aggregation function to be used for summarizing clipped data for each shape.	required
include_latlon	bool	Flag to include lat/lon in the dictionary entry.	required

Returns:

Name	Type	Description
dict	dict[str, Union[float, int]]	A dictionary with dot names as keys and calculated aggregations as values.

Source code in rastertoolkit/raster.py

def raster_clip_weighted_single(
    shp: ShapeView,
    sparse_weight: np.ndarray,
    sparse_val: np.ndarray,
    k1: int,
    shape_len: int,
    weight_summary_func: Callable,
    include_latlon: bool,
) -> dict[str, Union[float, int]]:
    """
    Extracts weighted data from a raster based on a single shape.

    Args:
        shp (ShapeView): Shape object.
        sparse_weight (np.ndarray): Sparse matrix of weight (population) raster data.
        sparse_val (np.ndarray): Sparse matrix of value raster data.
        k1 (int): Index of the shape.
        shape_len (int): Total number of shapes.
        weight_summary_func (Callable): Aggregation function to be used for summarizing
            clipped data for each shape.
        include_latlon (bool): Flag to include lat/lon in the dictionary entry.

    Returns:
        dict: A dictionary with dot names as keys and calculated aggregations as values.
    """
    data_dict = {}
    # Null shape; error in shapefile
    shp.validate()

    # Subset matrices for clipping
    weight_clip = subset_matrix_for_clipping(shape=shp, sparse_data=sparse_weight)
    val_clip = subset_matrix_for_clipping(shape=shp, sparse_data=sparse_val, pad=1)

    # Track booleans (indicates if lat/long is interior)
    data_bool = is_interior(shp, weight_clip)

    # Interpolate at population data
    final_val = interpolate_at_weight_data(shp, weight_clip, val_clip, data_bool)

    # Pop values
    values = weight_clip[data_bool, 2]

    # Entry dictionary
    weight_summary_func = weight_summary_func or default_summary_func
    entry = {"pop": weight_summary_func(values), "val": final_val}

    # Set entry and print status
    data_dict[shp.name] = summary_entry(shp, entry, include_latlon)
    print_status(shp, data_dict, k1, shape_len)

    return data_dict

subset_matrix_for_clipping(shape, sparse_data, pad=0)

Subset the matrix for clipping

Parameters:

Name	Type	Description	Default
shape	ShapeView	Shape object.	required
sparse_data	ndarray	Sparse matrix of raster data.	required
pad	int	Padding for clipping.	0

Returns:

Type	Description
ndarray	np.ndarray: A subset of the matrix for clipping.

Source code in rastertoolkit/raster.py

def subset_matrix_for_clipping(
    shape: ShapeView,
    sparse_data: np.ndarray,
    pad: int = 0,
) -> np.ndarray:
    """
    Subset the matrix for clipping

    Args:
        shape (ShapeView): Shape object.
        sparse_data (np.ndarray): Sparse matrix of raster data.
        pad (int): Padding for clipping.

    Returns:
        np.ndarray: A subset of the matrix for clipping.
    """
    clip_bool1 = np.logical_and(
        sparse_data[:, 0] > shape.xy_min[0] - pad,
        sparse_data[:, 1] > shape.xy_min[1] - pad,
    )
    clip_bool2 = np.logical_and(
        sparse_data[:, 0] < shape.xy_max[0] + pad,
        sparse_data[:, 1] < shape.xy_max[1] + pad,
    )
    data_clip = sparse_data[np.logical_and(clip_bool1, clip_bool2), :]

    return data_clip

summary_entry(shape, entry, include_latlon)

Summarize the entry for the shape.

Parameters:

Name	Type	Description	Default
shape	ShapeView	Shape object.	required
entry	Union[dict, float, int]	Entry for the shape.	required
include_latlon	bool	Flag to include lat/lon in the dictionary entry.	required

Returns:

Type	Description
Union[dict, float, int]	Union[dict, float, int]: The summarized entry for the shape.

Source code in rastertoolkit/raster.py

def summary_entry(
    shape: ShapeView,
    entry: Union[dict, float, int],
    include_latlon: bool,
) -> Union[dict, float, int]:
    """
    Summarize the entry for the shape.

    Args:
        shape (ShapeView): Shape object.
        entry (Union[dict, float, int]): Entry for the shape.
        include_latlon (bool): Flag to include lat/lon in the dictionary entry.

    Returns:
        Union[dict, float, int]: The summarized entry for the shape.
    """
    if include_latlon:
        if not isinstance(entry, dict) or len(entry) == 0:
            raise TypeError(
                f"Invalid entry; expected a non-empty dict when include_latlon is True, got {entry!r}."
            )
        lon = shape.center[0] if shape else np.nan
        lat = shape.center[1] if shape else np.nan
        final_entry = {"lat": lat, "lon": lon}
        final_entry.update(entry)
    else:
        if isinstance(entry, dict) and len(entry) == 1:
            final_entry = list(entry.values())[0]
        else:
            final_entry = entry

    return final_entry