myco.geo

Geospatial data read, write and manipulation functions

block_windows(width, height, blockxsize=BLOCKSIZE, blockysize=BLOCKSIZE)

Returns a tuple of iterables with custom window sizes for block-based read/write.

Designed to be used as a custom block size implementation of rasterio's src.block_windows().

Parameters:

Name	Type	Description	Default
width	int	the raster width (xsize) in pixels.	required
height	int	the raster height (ysize) in pixels.	required
blockxsize	int	the size of the x (column) dimension to read per block.	BLOCKSIZE
blockysize	int	the size of the y (row) dimension to read per block.	BLOCKSIZE

Returns:

Name	Type	Description
iter	windows	an iterable with a tuple of block indexes ((0, 0), (0, 1) ... ) and rasterio read Windows (Window(col_off=0, row_off=0, width=blockxsize, height=blockysize), ... ).

Source code in myco/geo.py

def block_windows(
    width: int, height: int, blockxsize: int = BLOCKSIZE, blockysize: int = BLOCKSIZE
) -> iter:
    """Returns a tuple of iterables with custom window sizes for block-based read/write.

    Designed to be used as a custom block size implementation of rasterio's `src.block_windows()`.

    Args:
        width: the raster width (xsize) in pixels.
        height: the raster height (ysize) in pixels.
        blockxsize: the size of the x (column) dimension to read per block.
        blockysize: the size of the y (row) dimension to read per block.

    Returns:
        iter(windows): an iterable with a tuple of block indexes ((0, 0), (0, 1) ... ) and
            rasterio read Windows (Window(col_off=0, row_off=0, width=blockxsize, height=blockysize), ... ).
    """

    # check to see if the blocks are perfect fits for the extent. if not, use custom tile boundary blocks
    xmod = width % blockxsize
    ymod = height % blockysize
    perfectx = xmod == 0
    perfecty = ymod == 0

    # get the total number of blocks to set
    n_xblocks = int(np.ceil(width / blockxsize))
    n_yblocks = int(np.ceil(height / blockysize))
    xblocks = range(n_xblocks)
    yblocks = range(n_yblocks)

    # iterate over each row/col tile combination to create the output idx, window iterator
    windows = []
    for yblock in yblocks:
        for xblock in xblocks:
            idx = (yblock, xblock)
            xoff = blockxsize * xblock
            yoff = blockysize * yblock

            window_width = blockxsize if (xblock != xblocks[-1] or perfectx) else xmod
            window_height = blockysize if (yblock != yblocks[-1] or perfecty) else ymod

            window = rio.windows.Window(xoff, yoff, window_width, window_height)
            windows.append((idx, window))

    return iter(windows)

get_fio_drivers(writeable=False)

Return a list of available fiona drivers

Source code in myco/geo.py

def get_fio_drivers(writeable: bool = False) -> List[str]:
    """Return a list of available fiona drivers"""
    all_drivers = list(fio.drvsupport.supported_drivers.keys())

    if writeable:
        drivers = [
            driver
            for driver in all_drivers
            if "w" in fio.drvsupport.supported_drivers[driver]
        ]
        return drivers

    else:
        return all_drivers

get_rio_drivers(writeable=False)

Return a list of available rasterio drivers

Source code in myco/geo.py

def get_rio_drivers(writeable: bool = False) -> List[str]:
    """Return a list of available rasterio drivers"""
    extensions = rio.drivers.raster_driver_extensions()
    all_drivers = [driver for ext, driver in extensions.items()]
    unique = list(np.unique(all_drivers)) + ["COG"]
    unique.sort()

    if writeable:
        drivers = [
            driver for driver in unique if not rio.drivers.is_blacklisted(driver, "w")
        ]
        return drivers

    else:
        return unique

get_rio_dtypes()

Returns a list of rasterio-format data type strings

Source code in myco/geo.py

def get_rio_dtypes() -> List[str]:
    """Returns a list of rasterio-format data type strings"""
    return list(rio.dtypes.dtype_ranges.keys())

get_rio_resampling()

Returns a list of rasterio-format resampling options

Source code in myco/geo.py

def get_rio_resampling() -> List[str]:
    """Returns a list of rasterio-format resampling options"""
    options = []
    for option in dir(rio.enums.Resampling):
        if "__" not in option:
            options.append(option)

    return options

interior_window_idxs(window_shape, percent=0.5)

Computes the xstart:xend, ystart:yend indices for the interior window of a 2d array.

Parameters:

Name	Type	Description	Default
window_shape	tuple	2d tuple of a window size (e.g. (64, 64))	required
percent	float	the percent of the interior window to include.	0.5

Returns:

Type	Description
tuple	ymin, ymax, xmin, xmax indices for slicing the interior windows.

Source code in myco/geo.py

def interior_window_idxs(window_shape: tuple, percent: float = 0.5) -> tuple:
    """Computes the xstart:xend, ystart:yend indices for the interior window
        of a 2d array.

    Args:
        window_shape: 2d tuple of a window size (e.g. (64, 64))
        percent: the percent of the interior window to include.

    Returns:
        ymin, ymax, xmin, xmax indices for slicing the interior windows.
    """
    ysize, xsize = window_shape
    xrad = (xsize * percent) / 2
    yrad = (ysize * percent) / 2
    xmid = xsize / 2
    ymid = ysize / 2
    xmin = int(xmid - xrad)
    xmax = int(xmid + xrad)
    ymin = int(ymid - yrad)
    ymax = int(ymid + yrad)

    return ymin, ymax, xmin, xmax

is_raster(filepath)

Verify a file is a valid, readable raster dataset

Source code in myco/geo.py

def is_raster(filepath: str) -> bool:
    """Verify a file is a valid, readable raster dataset"""
    try:
        with rio.open(filepath, "r") as src:
            pass
        return True
    except rio.errors.RasterioIOError:
        return False

is_vector(filepath)

Verify a sile is a valid, readable vector dataset

Source code in myco/geo.py

def is_vector(filepath: str) -> bool:
    """Verify a sile is a valid, readable vector dataset"""
    try:
        with fio.open(filepath, "r") as src:
            pass
        return True
    except fio.errors.DriverError:
        return False

point_to_window_bounds(point, xsize, xres, ysize=None, yres=None)

Computes the tile bounds for a window surrounding a point.

Parameters:

Name	Type	Description	Default
point	Point	the point sample centroid.	required
xsize	int	the number of x pixels.	required
xres	float	the size (resolution) of x pixels.	required
ysize	int	the number of y pixels. uses xsize if not passed.	None
yres	float	the size (resolution) of y pixels. uses ysize if not passed.	None

Returns:

Type	Description
Tuple[float, float, float, float]	[xmin, ymin, xmax, ymax] of the bounding box centered on the point.

Source code in myco/geo.py

def point_to_window_bounds(
    point: Point, xsize: int, xres: float, ysize: int = None, yres: float = None
) -> Tuple[float, float, float, float]:
    """Computes the tile bounds for a window surrounding a point.

    Args:
        point: the point sample centroid.
        xsize: the number of x pixels.
        xres: the size (resolution) of x pixels.
        ysize: the number of y pixels.
            uses xsize if not passed.
        yres: the size (resolution) of y pixels.
            uses ysize if not passed.

    Returns:
        [xmin, ymin, xmax, ymax] of the bounding box centered on the point.
    """
    ysize = ysize or xsize
    yres = yres or xres
    width = xsize * xres
    height = ysize * yres
    halfwidth = width / 2.0
    halfheight = height / 2.0
    halfxpix = xres / 2.0
    halfypix = yres / 2.0
    x, y = point.x, point.y
    xmin = x - halfwidth - halfxpix
    ymin = y - halfwidth
    xmax = x + halfwidth - halfxpix
    ymax = y + halfwidth

    return xmin, ymin, xmax, ymax

raster_interior_windows(overlapping_windows, interior_idxs)

Get interior windows for data that match external read blocks

Source code in myco/geo.py

def raster_interior_windows(
    overlapping_windows: List[rio.windows.Window],
    interior_idxs: Tuple[int, int, int, int],
):
    """Get interior windows for data that match external read blocks"""

    # set up the base geometry for interior and exterior windows
    iymin, iymax, ixmin, ixmax = interior_idxs
    iheight = iymax - iymin
    iwidth = ixmax - ixmin

    interior_windows = []
    for window in overlapping_windows:
        wxoff = window.col_off + ixmin
        wyoff = window.row_off + iymin

        interior_window = rio.windows.Window(wxoff, wyoff, iwidth, iheight)
        interior_windows.append(interior_window)

    return interior_windows

raster_overlapping_windows(width, height, window_shape, interior_idxs)

Creates a list of block windows the size of the interior model application window.

Source code in myco/geo.py

def raster_overlapping_windows(
    width: int,
    height: int,
    window_shape: Tuple[int, int],
    interior_idxs: Tuple[int, int, int, int],
):
    """Creates a list of block windows the size of the interior model application window."""

    # set up the base geometry for interior and exterior windows
    iymin, iymax, ixmin, ixmax = interior_idxs
    iheight = iymax - iymin
    iwidth = ixmax - ixmin
    oheight, owidth = window_shape
    left_pad = ixmin
    right_pad = owidth - ixmax
    top_pad = iymin
    bottom_pad = oheight - iymax
    interior_width = width - right_pad - left_pad
    interior_height = height - top_pad - bottom_pad

    # get the windows for the interior subset
    interior_windows = [
        window
        for _, window in block_windows(interior_width, interior_height, iwidth, iheight)
    ]

    # pad them back to the full external width
    exterior_windows = []
    for window in interior_windows:
        wxoff = window.col_off
        wyoff = window.row_off
        wwidth = window.width + left_pad + right_pad
        wheight = window.height + top_pad + bottom_pad

        # adjust edge windows to still able to read a full-width image
        if wheight < oheight:
            yoffset = oheight % wheight
            wyoff -= yoffset
            wheight = oheight

        if wwidth < owidth:
            xoffset = owidth % wwidth
            wxoff -= xoffset
            wwidth = owidth

        if wheight + wyoff > height:
            yoffset = (wheight + wyoff) % height
            wyoff -= yoffset

        if wwidth + wxoff > width:
            xoffset = (wwidth + wxoff) % width
            wxoff -= xoffset

        exterior_window = rio.windows.Window(wxoff, wyoff, wwidth, wheight)
        exterior_windows.append(exterior_window)

    return exterior_windows