io

BackgroundBlockWriter

Bases: BackgroundWriter

Class to write data to multiple files in the background using gdal bindings.

Source code in src/dolphin/io/_writers.py

class BackgroundBlockWriter(BackgroundWriter):
    """Class to write data to multiple files in the background using `gdal` bindings."""

    def __init__(
        self,
        *,
        max_queue: int = 0,
        debug: bool = False,
        keep_bits: int | None = None,
        **kwargs,
    ):
        super().__init__(nq=max_queue, name="Writer")
        self.keep_bits = keep_bits
        if debug:
            #  background thread. Just synchronously write data
            self.notify_finished()
            self.queue_write = self.write  # type: ignore[assignment]

    def write(
        self,
        data: ArrayLike,
        filename: Filename,
        row_start: int,
        col_start: int,
        band: int | None = None,
    ):
        """Write out an ndarray to a subset of the pre-made `filename`.

        Parameters
        ----------
        data : ArrayLike
            2D or 3D data array to save.
        filename : Filename
            list of output files to save to, or (if cur_block is 2D) a single file.
        row_start : int
            Row index to start writing at.
        col_start : int
            Column index to start writing at.
        band : int, optional
            Band index in the file to write. Defaults to None, which uses first band,
            or whole raster for 3D data.

        Raises
        ------
        ValueError
            If length of `output_files` does not match length of `cur_block`.

        """
        from dolphin.io import write_block

        if np.issubdtype(data.dtype, np.floating) and self.keep_bits is not None:
            round_mantissa(data, keep_bits=self.keep_bits)

        write_block(data, filename, row_start, col_start, band=band)

write(data, filename, row_start, col_start, band=None)

Write out an ndarray to a subset of the pre-made filename.

Parameters:

Name	Type	Description	Default
data	ArrayLike	2D or 3D data array to save.	required
filename	Filename	list of output files to save to, or (if cur_block is 2D) a single file.	required
row_start	int	Row index to start writing at.	required
col_start	int	Column index to start writing at.	required
band	int	Band index in the file to write. Defaults to None, which uses first band, or whole raster for 3D data.	None

Raises:

Type	Description
ValueError	If length of output_files does not match length of cur_block.

Source code in src/dolphin/io/_writers.py

def write(
    self,
    data: ArrayLike,
    filename: Filename,
    row_start: int,
    col_start: int,
    band: int | None = None,
):
    """Write out an ndarray to a subset of the pre-made `filename`.

    Parameters
    ----------
    data : ArrayLike
        2D or 3D data array to save.
    filename : Filename
        list of output files to save to, or (if cur_block is 2D) a single file.
    row_start : int
        Row index to start writing at.
    col_start : int
        Column index to start writing at.
    band : int, optional
        Band index in the file to write. Defaults to None, which uses first band,
        or whole raster for 3D data.

    Raises
    ------
    ValueError
        If length of `output_files` does not match length of `cur_block`.

    """
    from dolphin.io import write_block

    if np.issubdtype(data.dtype, np.floating) and self.keep_bits is not None:
        round_mantissa(data, keep_bits=self.keep_bits)

    write_block(data, filename, row_start, col_start, band=band)

BackgroundRasterWriter

Bases: BackgroundWriter, DatasetWriter

Class to write data to files in a background thread.

Source code in src/dolphin/io/_writers.py

class BackgroundRasterWriter(BackgroundWriter, DatasetWriter):
    """Class to write data to files in a background thread."""

    def __init__(
        self,
        filename: Filename,
        *,
        max_queue: int = 0,
        debug: bool = False,
        keep_bits: int | None = None,
        **kwargs,
    ):
        super().__init__(nq=max_queue, name="Writer")
        if debug:
            #  background thread. Just synchronously write data
            self.notify_finished()
            self.queue_write = self.write  # type: ignore[assignment]

        if Path(filename).exists():
            self._raster = RasterWriter(filename, keep_bits=keep_bits)
        else:
            self._raster = RasterWriter.create(filename, keep_bits=keep_bits, **kwargs)
        self.filename = filename
        self.ndim = 2

    def write(self, key: tuple[Index, ...], value: np.ndarray):
        """Write out an ndarray to a subset of the pre-made `filename`.

        Parameters
        ----------
        key : tuple[Index,...]
            The key of the data to write.

        value : np.ndarray
            The block of data to write.

        """
        self._raster[key] = value

    def __setitem__(self, key: tuple[Index, ...], value: np.ndarray, /) -> None:
        self.queue_write(key, value)

    def close(self):
        """Close the underlying dataset and stop the background thread."""
        self._raster.close()
        self.notify_finished()

    @property
    def closed(self) -> bool:
        """bool : True if the dataset is closed."""  # noqa: D403
        return self._raster.closed

    @property
    def shape(self):
        return self._raster.shape

    @property
    def dtype(self) -> np.dtype:
        return self._raster.dtype

    def __enter__(self):
        return self

    def __exit__(self, exc_type, exc_value, traceback):  # type: ignore[no-untyped-def]
        self.close()

closed property

bool : True if the dataset is closed.

close()

Close the underlying dataset and stop the background thread.

Source code in src/dolphin/io/_writers.py

def close(self):
    """Close the underlying dataset and stop the background thread."""
    self._raster.close()
    self.notify_finished()

write(key, value)

Write out an ndarray to a subset of the pre-made filename.

Parameters:

Name	Type	Description	Default
key	tuple[Index, ...]	The key of the data to write.	required
value	ndarray	The block of data to write.	required

Source code in src/dolphin/io/_writers.py

def write(self, key: tuple[Index, ...], value: np.ndarray):
    """Write out an ndarray to a subset of the pre-made `filename`.

    Parameters
    ----------
    key : tuple[Index,...]
        The key of the data to write.

    value : np.ndarray
        The block of data to write.

    """
    self._raster[key] = value

BackgroundReader

Bases: BackgroundWorker

Base class for reading data in a background thread (pre-fetching).

After instantiating an object, a client sends it data selection parameters (slices, indices, etc.) via the queue_read method and retrieves the result with the get_data method. In order to get useful concurrency, that usually means you'll want to queue the read for the next data block before starting work on the current block. The reader remains active until notify_finished is called and all blocks have been retrieved. Subclasses must define the read method.

Parameters:

Name	Type	Description	Default
nq	int	Number of read results that can be stored before blocking, <= 0 for unbounded. Default is 1.	1
timeout	float	Interval in seconds used to check for finished notification once write queue is empty.	_DEFAULT_TIMEOUT

Source code in src/dolphin/io/_background.py

class BackgroundReader(BackgroundWorker):
    """Base class for reading data in a background thread (pre-fetching).

    After instantiating an object, a client sends it data selection parameters
    (slices, indices, etc.) via the `queue_read` method and retrieves the result
    with the `get_data` method.  In order to get useful concurrency, that
    usually means you'll want to queue the read for the next data block before
    starting work on the current block.  The reader remains active until
    `notify_finished` is called and all blocks have been retrieved.  Subclasses
    must define the `read` method.

    Parameters
    ----------
    nq : int
        Number of read results that can be stored before blocking, <= 0 for
        unbounded.  Default is 1.
    timeout : float
        Interval in seconds used to check for finished notification once write
        queue is empty.

    """

    def __init__(self, nq=1, timeout=_DEFAULT_TIMEOUT, **kwargs):
        super().__init__(
            num_results_queue=nq,
            timeout=timeout,
            store_results=True,
            # If we're reading data, we don't care about the result queue
            drop_unfinished_results=True,
            **kwargs,
        )

    # rename queue_work -> queue_read
    def queue_read(self, *args, **kw):
        """Add selection parameters (slices, etc.) to the read queue to be processed.

        Same input interface as `read`.
        """
        self.queue_work(*args, **kw)

    # rename get_result -> get_data
    def get_data(self):
        """Retrieve the least-recently read chunk of data.

        Blocks until a result is available.
        Same output interface as `read`.
        """
        return self.get_result()

    # rename process -> read
    def process(self, *args, **kw):
        return self.read(*args, **kw)

    @abc.abstractmethod
    def read(self, *args, **kw):
        """User-defined method for reading a chunk of data."""

get_data()

Retrieve the least-recently read chunk of data.

Blocks until a result is available. Same output interface as read.

Source code in src/dolphin/io/_background.py

def get_data(self):
    """Retrieve the least-recently read chunk of data.

    Blocks until a result is available.
    Same output interface as `read`.
    """
    return self.get_result()

queue_read(*args, **kw)

Add selection parameters (slices, etc.) to the read queue to be processed.

Same input interface as read.

Source code in src/dolphin/io/_background.py

def queue_read(self, *args, **kw):
    """Add selection parameters (slices, etc.) to the read queue to be processed.

    Same input interface as `read`.
    """
    self.queue_work(*args, **kw)

read(*args, **kw) abstractmethod

User-defined method for reading a chunk of data.

Source code in src/dolphin/io/_background.py

@abc.abstractmethod
def read(self, *args, **kw):
    """User-defined method for reading a chunk of data."""

BackgroundStackWriter

Bases: BackgroundWriter, DatasetStackWriter

Class to write 3D data to multiple files in a background thread.

Will create/overwrite the files in file_list if they exist.

Source code in src/dolphin/io/_writers.py

class BackgroundStackWriter(BackgroundWriter, DatasetStackWriter):
    """Class to write 3D data to multiple files in a background thread.

    Will create/overwrite the files in `file_list` if they exist.
    """

    def __init__(
        self,
        file_list: Sequence[Filename],
        *,
        like_filename: Filename | None = None,
        max_queue: int = 0,
        debug: bool = False,
        keep_bits: int | None = None,
        **file_creation_kwargs,
    ):
        from dolphin.io import write_arr

        super().__init__(nq=max_queue, name="GdalStackWriter")
        if debug:
            # Stop background thread. Just synchronously write data
            self.notify_finished()
            self.queue_write = self.write  # type: ignore[assignment]

        for fn in file_list:
            write_arr(
                arr=None,
                output_name=fn,
                like_filename=like_filename,
                **file_creation_kwargs,
            )

        self.file_list = file_list
        self.keep_bits = keep_bits

        with rasterio.open(self.file_list[0]) as src:
            self.shape = (len(self.file_list), *src.shape)
            self.dtype = src.dtypes[0]

    def write(
        self, data: ArrayLike, row_start: int, col_start: int, band: int | None = None
    ):
        """Write out an ndarray to a subset of the pre-made `filename`.

        Parameters
        ----------
        data : ArrayLike
            3D data array to save.
        row_start : int
            Row index to start writing at.
        col_start : int
            Column index to start writing at.
        band : int, optional
            Band index in the file to write. Defaults to None, which uses first band,
            or whole raster for 3D data.

        Raises
        ------
        ValueError
            If length of `output_files` does not match length of `cur_block`.

        """
        from dolphin.io import write_block

        _do_round = (
            np.issubdtype(data.dtype, np.floating) and self.keep_bits is not None
        )
        if data.ndim == 2:
            data = data[None, ...]
        if data.shape[0] != len(self.file_list):
            raise ValueError(f"{data.shape = }, but {len(self.file_list) = }")
        for fn, layer in zip(self.file_list, data, strict=False):
            if _do_round:
                assert self.keep_bits is not None
                round_mantissa(layer, keep_bits=self.keep_bits)
            write_block(layer, fn, row_start, col_start, band=band)

    def __setitem__(self, key, value):
        # Unpack the slices
        band, rows, cols = _unpack_3d_slices(key)
        # Check we asked to write to all the files
        if band not in (slice(None), slice(None, None, None), ...):
            self.notify_finished()
            raise NotImplementedError("Can only write to all files at once.")
        self.queue_write(value, rows.start, cols.start)

    @property
    def closed(self) -> bool:
        """bool : True if the dataset is closed."""  # noqa: D403
        return self._thread.is_alive() is False

    def close(self) -> None:
        """Close the underlying dataset and stop the background thread."""
        self.notify_finished()

closed property

bool : True if the dataset is closed.

close()

Close the underlying dataset and stop the background thread.

Source code in src/dolphin/io/_writers.py

def close(self) -> None:
    """Close the underlying dataset and stop the background thread."""
    self.notify_finished()

write(data, row_start, col_start, band=None)

Write out an ndarray to a subset of the pre-made filename.

Parameters:

Name	Type	Description	Default
data	ArrayLike	3D data array to save.	required
row_start	int	Row index to start writing at.	required
col_start	int	Column index to start writing at.	required
band	int	Band index in the file to write. Defaults to None, which uses first band, or whole raster for 3D data.	None

Raises:

Type	Description
ValueError	If length of output_files does not match length of cur_block.

Source code in src/dolphin/io/_writers.py

def write(
    self, data: ArrayLike, row_start: int, col_start: int, band: int | None = None
):
    """Write out an ndarray to a subset of the pre-made `filename`.

    Parameters
    ----------
    data : ArrayLike
        3D data array to save.
    row_start : int
        Row index to start writing at.
    col_start : int
        Column index to start writing at.
    band : int, optional
        Band index in the file to write. Defaults to None, which uses first band,
        or whole raster for 3D data.

    Raises
    ------
    ValueError
        If length of `output_files` does not match length of `cur_block`.

    """
    from dolphin.io import write_block

    _do_round = (
        np.issubdtype(data.dtype, np.floating) and self.keep_bits is not None
    )
    if data.ndim == 2:
        data = data[None, ...]
    if data.shape[0] != len(self.file_list):
        raise ValueError(f"{data.shape = }, but {len(self.file_list) = }")
    for fn, layer in zip(self.file_list, data, strict=False):
        if _do_round:
            assert self.keep_bits is not None
            round_mantissa(layer, keep_bits=self.keep_bits)
        write_block(layer, fn, row_start, col_start, band=band)

BackgroundWorker

Bases: ABC

Base class for doing work in a background thread.

After instantiating an object, a client sends it work with the queue_work method and retrieves the result with the get_result method (hopefully after doing something else useful in between). The worker remains active until notify_finished is called. Subclasses must define the process method.

Parameters:

Name	Type	Description	Default
num_work_queue	int	Max number of work items to queue before blocking, <= 0 for unbounded.	0
num_results_queue	int	Max number of results to generate before blocking, <= 0 for unbounded.	0
store_results	bool	Whether to store return values of process method. If True then get_result must be called once for every queue_work call.	True
timeout	float	Interval in seconds used to check for finished notification once work queue is empty.	_DEFAULT_TIMEOUT

Notes

The usual caveats about Python threading apply. It's typically a poor choice for concurrency unless the global interpreter lock (GIL) has been released, which can happen in IO calls and compiled extensions.

Source code in src/dolphin/io/_background.py

class BackgroundWorker(abc.ABC):
    """Base class for doing work in a background thread.

    After instantiating an object, a client sends it work with the `queue_work`
    method and retrieves the result with the `get_result` method (hopefully
    after doing something else useful in between).  The worker remains active
    until `notify_finished` is called.  Subclasses must define the `process`
    method.

    Parameters
    ----------
    num_work_queue : int
        Max number of work items to queue before blocking, <= 0 for unbounded.
    num_results_queue : int
        Max number of results to generate before blocking, <= 0 for unbounded.
    store_results : bool
        Whether to store return values of `process` method.  If True then
        `get_result` must be called once for every `queue_work` call.
    timeout : float
        Interval in seconds used to check for finished notification once work
        queue is empty.

    Notes
    -----
    The usual caveats about Python threading apply.  It's typically a poor
    choice for concurrency unless the global interpreter lock (GIL) has been
    released, which can happen in IO calls and compiled extensions.

    """

    def __init__(
        self,
        num_work_queue=0,
        num_results_queue=0,
        store_results=True,
        drop_unfinished_results=False,
        timeout=_DEFAULT_TIMEOUT,
        name="BackgroundWorker",
    ):
        self.name = name
        self.store_results = store_results
        self.timeout = timeout
        self._finished_event = Event()
        self._work_queue = Queue(num_work_queue)
        if self.store_results:
            self._results_queue = Queue(num_results_queue)
        self._thread = Thread(target=self._consume_work_queue, name=name)
        self._thread.start()
        self._drop_unfinished_results = drop_unfinished_results

    def _consume_work_queue(self):
        while True:
            if not main_thread().is_alive():
                break

            logger.debug(f"{self.name} getting work")
            if self._finished_event.is_set():
                do_exit = self._drop_unfinished_results or (
                    self._work_queue.unfinished_tasks == 0
                )
                if do_exit:
                    break
                else:
                    # Keep going even if finished event is set
                    logger.debug(
                        f"{self.name} Finished... but waiting for work queue to empty,"
                        f" {self._work_queue.qsize()} items left,"
                        f" {self._work_queue.unfinished_tasks} unfinished"
                    )
            try:
                args, kw = self._work_queue.get(timeout=self.timeout)
                logger.debug(f"{self.name} processing")
                result = self.process(*args, **kw)
                self._work_queue.task_done()
                # Notify the queue that processing is done
                logger.debug(f"{self.name} got result")
            except Empty:
                logger.debug(f"{self.name} timed out, checking if done")
                continue

            if self.store_results:
                logger.debug(f"{self.name} saving result in queue")
                while True:
                    try:
                        self._results_queue.put(result, timeout=2)
                        break
                    except Full:
                        logger.debug(f"{self.name} result queue full, waiting...")
                        continue

    @abc.abstractmethod
    def process(self, *args, **kw):
        """User-defined task to operate in background thread."""

    def queue_work(self, *args, **kw):
        """Add a job to the work queue to be executed.

        Blocks if work queue is full.
        Same input interface as `process`.
        """
        if self._finished_event.is_set():
            msg = "Attempted to queue_work after notify_finished!"
            raise RuntimeError(msg)
        self._work_queue.put((args, kw))

    def get_result(self):
        """Get the least-recent value from the result queue.

        Blocks until a result is available.
        Same output interface as `process`.
        """
        while True:
            try:
                result = self._results_queue.get(timeout=self.timeout)
                self._results_queue.task_done()
                break
            except Empty as e:
                logger.debug(f"{self.name} get_result timed out, checking if done")
                if self._finished_event.is_set():
                    msg = "Attempted to get_result after notify_finished!"
                    raise RuntimeError(msg) from e
                continue
        return result

    def notify_finished(self, timeout=None):
        """Signal that all work has finished.

        Indicate that no more work will be added to the queue, and block until
        all work has been processed.
        If `store_results=True` also block until all results have been retrieved.
        """
        self._finished_event.set()
        if self.store_results and not self._drop_unfinished_results:
            self._results_queue.join()
        self._thread.join(timeout)

    def __del__(self):
        logger.debug(f"{self.name} notifying of exit")
        self.notify_finished()

get_result()

Get the least-recent value from the result queue.

Blocks until a result is available. Same output interface as process.

Source code in src/dolphin/io/_background.py

def get_result(self):
    """Get the least-recent value from the result queue.

    Blocks until a result is available.
    Same output interface as `process`.
    """
    while True:
        try:
            result = self._results_queue.get(timeout=self.timeout)
            self._results_queue.task_done()
            break
        except Empty as e:
            logger.debug(f"{self.name} get_result timed out, checking if done")
            if self._finished_event.is_set():
                msg = "Attempted to get_result after notify_finished!"
                raise RuntimeError(msg) from e
            continue
    return result

notify_finished(timeout=None)

Signal that all work has finished.

Indicate that no more work will be added to the queue, and block until all work has been processed. If store_results=True also block until all results have been retrieved.

Source code in src/dolphin/io/_background.py

def notify_finished(self, timeout=None):
    """Signal that all work has finished.

    Indicate that no more work will be added to the queue, and block until
    all work has been processed.
    If `store_results=True` also block until all results have been retrieved.
    """
    self._finished_event.set()
    if self.store_results and not self._drop_unfinished_results:
        self._results_queue.join()
    self._thread.join(timeout)

process(*args, **kw) abstractmethod

User-defined task to operate in background thread.

Source code in src/dolphin/io/_background.py

@abc.abstractmethod
def process(self, *args, **kw):
    """User-defined task to operate in background thread."""

queue_work(*args, **kw)

Add a job to the work queue to be executed.

Blocks if work queue is full. Same input interface as process.

Source code in src/dolphin/io/_background.py

def queue_work(self, *args, **kw):
    """Add a job to the work queue to be executed.

    Blocks if work queue is full.
    Same input interface as `process`.
    """
    if self._finished_event.is_set():
        msg = "Attempted to queue_work after notify_finished!"
        raise RuntimeError(msg)
    self._work_queue.put((args, kw))

BackgroundWriter

Bases: BackgroundWorker

Base class for writing data in a background thread.

After instantiating an object, a client sends it data with the queue_write method. The writer remains active until notify_finished is called. Subclasses must define the write method.

Parameters:

Name	Type	Description	Default
nq	int	Number of write jobs that can be queued before blocking, <= 0 for unbounded. Default is 1.	1
timeout	float	Interval in seconds used to check for finished notification once write queue is empty.	_DEFAULT_TIMEOUT

Source code in src/dolphin/io/_background.py

class BackgroundWriter(BackgroundWorker):
    """Base class for writing data in a background thread.

    After instantiating an object, a client sends it data with the `queue_write`
    method.  The writer remains active until `notify_finished` is called.
    Subclasses must define the `write` method.

    Parameters
    ----------
    nq : int
        Number of write jobs that can be queued before blocking, <= 0 for
        unbounded.  Default is 1.
    timeout : float
        Interval in seconds used to check for finished notification once write
        queue is empty.

    """

    def __init__(self, nq=1, timeout=_DEFAULT_TIMEOUT, **kwargs):
        super().__init__(
            num_work_queue=nq,
            store_results=False,
            timeout=timeout,
            **kwargs,
        )

    # rename queue_work -> queue_write
    def queue_write(self, *args, **kw):
        """Add data to the queue to be written.

        Blocks if write queue is full.
        Same interfaces as `write`.
        """
        self.queue_work(*args, **kw)

    # rename process -> write
    def process(self, *args, **kw):
        self.write(*args, **kw)

    @property
    def num_queued(self):
        """Number of items waiting in the queue to be written."""
        return self._work_queue.qsize()

    @abc.abstractmethod
    def write(self, *args, **kw):
        """User-defined method for writing data."""

num_queued property

Number of items waiting in the queue to be written.

queue_write(*args, **kw)

Add data to the queue to be written.

Blocks if write queue is full. Same interfaces as write.

Source code in src/dolphin/io/_background.py

def queue_write(self, *args, **kw):
    """Add data to the queue to be written.

    Blocks if write queue is full.
    Same interfaces as `write`.
    """
    self.queue_work(*args, **kw)

write(*args, **kw) abstractmethod

User-defined method for writing data.

Source code in src/dolphin/io/_background.py

@abc.abstractmethod
def write(self, *args, **kw):
    """User-defined method for writing data."""

BinaryReader dataclass

Bases: DatasetReader

A flat binary file for storing array data.

See Also

HDF5Dataset RasterReader

Notes

This class does not store an open file object. Instead, the file is opened on-demand for reading or writing and closed immediately after each read/write operation. This allows multiple spawned processes to write to the file in coordination (as long as a suitable mutex is used to guard file access.)

Source code in src/dolphin/io/_readers.py

@dataclass
class BinaryReader(DatasetReader):
    """A flat binary file for storing array data.

    See Also
    --------
    HDF5Dataset
    RasterReader

    Notes
    -----
    This class does not store an open file object. Instead, the file is opened on-demand
    for reading or writing and closed immediately after each read/write operation. This
    allows multiple spawned processes to write to the file in coordination (as long as a
    suitable mutex is used to guard file access.)

    """

    filename: Path
    """pathlib.Path : The file path."""

    shape: tuple[int, ...]
    """tuple of int : Tuple of array dimensions."""

    dtype: np.dtype
    """numpy.dtype : Data-type of the array's elements."""

    nodata: Optional[float] = None
    """Optional[float] : Value to use for nodata pixels."""

    def __post_init__(self):
        self.filename = Path(self.filename)
        if not self.filename.exists():
            msg = f"File {self.filename} does not exist."
            raise FileNotFoundError(msg)
        self.dtype = np.dtype(self.dtype)

    @property
    def ndim(self) -> int:  # type: ignore[override]
        """int : Number of array dimensions."""  # noqa: D403
        return len(self.shape)

    def __getitem__(self, key: tuple[Index, ...], /) -> np.ndarray:
        with self.filename.open("rb") as f:  # noqa: SIM117
            # Memory-map the entire file.
            with mmap.mmap(f.fileno(), 0, access=mmap.ACCESS_READ) as mm:
                # In order to safely close the memory-map, there can't be any dangling
                # references to it, so we return a copy (not a view) of the requested
                # data and decref the array object.
                arr = np.frombuffer(mm, dtype=self.dtype).reshape(self.shape)
                data = arr[key].copy()
                del arr
        return _mask_array(data, self.nodata) if self.nodata is not None else data

    def __array__(self) -> np.ndarray:
        return self[:,]

    @classmethod
    def from_gdal(
        cls, filename: Filename, band: int = 1, nodata: Optional[float] = None
    ) -> BinaryReader:
        """Create a BinaryReader from a GDAL-readable file.

        Parameters
        ----------
        filename : Filename
            Path to the file to read.
        band : int, optional
            Band to read from the file, by default 1
        nodata : float, optional
            Value to use for nodata pixels, by default None
            If None passed, will search for a nodata value in the file.

        Returns
        -------
        BinaryReader
            The BinaryReader object.

        """
        with rio.open(filename) as src:
            dtype = src.dtypes[band - 1]
            shape = src.shape
            nodata = src.nodatavals[band - 1]
        return cls(
            Path(filename),
            shape=shape,
            dtype=dtype,
            nodata=nodata or nodata,
        )

dtype instance-attribute

numpy.dtype : Data-type of the array's elements.

filename instance-attribute

pathlib.Path : The file path.

ndim property

int : Number of array dimensions.

nodata = None class-attribute instance-attribute

Optional[float] : Value to use for nodata pixels.

shape instance-attribute

tuple of int : Tuple of array dimensions.

from_gdal(filename, band=1, nodata=None) classmethod

Create a BinaryReader from a GDAL-readable file.

Parameters:

Name	Type	Description	Default
filename	Filename	Path to the file to read.	required
band	int	Band to read from the file, by default 1	1
nodata	float	Value to use for nodata pixels, by default None If None passed, will search for a nodata value in the file.	None

Returns:

Type	Description
BinaryReader	The BinaryReader object.

Source code in src/dolphin/io/_readers.py

@classmethod
def from_gdal(
    cls, filename: Filename, band: int = 1, nodata: Optional[float] = None
) -> BinaryReader:
    """Create a BinaryReader from a GDAL-readable file.

    Parameters
    ----------
    filename : Filename
        Path to the file to read.
    band : int, optional
        Band to read from the file, by default 1
    nodata : float, optional
        Value to use for nodata pixels, by default None
        If None passed, will search for a nodata value in the file.

    Returns
    -------
    BinaryReader
        The BinaryReader object.

    """
    with rio.open(filename) as src:
        dtype = src.dtypes[band - 1]
        shape = src.shape
        nodata = src.nodatavals[band - 1]
    return cls(
        Path(filename),
        shape=shape,
        dtype=dtype,
        nodata=nodata or nodata,
    )

BinaryStackReader dataclass

Bases: BaseStackReader

Source code in src/dolphin/io/_readers.py

@dataclass
class BinaryStackReader(BaseStackReader):
    @classmethod
    def from_file_list(
        cls, file_list: Iterable[Filename], shape_2d: tuple[int, int], dtype: np.dtype
    ) -> BinaryStackReader:
        """Create a BinaryStackReader from a list of files.

        Parameters
        ----------
        file_list : Iterable[Filename]
            Iterable of paths to the files to read.
        shape_2d : tuple[int, int]
            Shape of each file.
        dtype : np.dtype
            Data type of each file.

        Returns
        -------
        BinaryStackReader
            The BinaryStackReader object.

        """
        files = list(file_list)
        readers = [BinaryReader(Path(f), shape=shape_2d, dtype=dtype) for f in files]
        return cls(file_list=files, readers=readers, num_threads=1)

    @classmethod
    def from_gdal(
        cls,
        file_list: Sequence[Filename],
        band: int = 1,
        num_threads: int = 1,
        nodata: Optional[float] = None,
    ) -> BinaryStackReader:
        """Create a BinaryStackReader from a list of GDAL-readable files.

        Parameters
        ----------
        file_list : Sequence[Filename]
            List of paths to the files to read.
        band : int, optional
            Band to read from the file, by default 1
        num_threads : int, optional (default 1)
            Number of threads to use for reading.
        nodata : float, optional
            Manually set value to use for nodata pixels, by default None
            If None passed, will search for a nodata value in the file.

        Returns
        -------
        BinaryStackReader
            The BinaryStackReader object.

        """
        readers = []
        dtypes = set()
        shapes = set()
        for f in file_list:
            with rio.open(f) as src:
                dtypes.add(src.dtypes[band - 1])
                shapes.add(src.shape)
            if len(dtypes) > 1:
                msg = "All files must have the same data type."
                raise ValueError(msg)
            if len(shapes) > 1:
                msg = "All files must have the same shape."
                raise ValueError(msg)
            readers.append(BinaryReader.from_gdal(f, band=band))
        return cls(
            file_list=file_list,
            readers=readers,
            num_threads=num_threads,
            nodata=nodata,
        )

from_file_list(file_list, shape_2d, dtype) classmethod

Create a BinaryStackReader from a list of files.

Parameters:

Name	Type	Description	Default
file_list	Iterable[Filename]	Iterable of paths to the files to read.	required
shape_2d	tuple[int, int]	Shape of each file.	required
dtype	dtype	Data type of each file.	required

Returns:

Type	Description
BinaryStackReader	The BinaryStackReader object.

Source code in src/dolphin/io/_readers.py

@classmethod
def from_file_list(
    cls, file_list: Iterable[Filename], shape_2d: tuple[int, int], dtype: np.dtype
) -> BinaryStackReader:
    """Create a BinaryStackReader from a list of files.

    Parameters
    ----------
    file_list : Iterable[Filename]
        Iterable of paths to the files to read.
    shape_2d : tuple[int, int]
        Shape of each file.
    dtype : np.dtype
        Data type of each file.

    Returns
    -------
    BinaryStackReader
        The BinaryStackReader object.

    """
    files = list(file_list)
    readers = [BinaryReader(Path(f), shape=shape_2d, dtype=dtype) for f in files]
    return cls(file_list=files, readers=readers, num_threads=1)

from_gdal(file_list, band=1, num_threads=1, nodata=None) classmethod

Create a BinaryStackReader from a list of GDAL-readable files.

Parameters:

Name	Type	Description	Default
file_list	Sequence[Filename]	List of paths to the files to read.	required
band	int	Band to read from the file, by default 1	1
num_threads	int, optional (default 1)	Number of threads to use for reading.	1
nodata	float	Manually set value to use for nodata pixels, by default None If None passed, will search for a nodata value in the file.	None

Returns:

Type	Description
BinaryStackReader	The BinaryStackReader object.

Source code in src/dolphin/io/_readers.py

@classmethod
def from_gdal(
    cls,
    file_list: Sequence[Filename],
    band: int = 1,
    num_threads: int = 1,
    nodata: Optional[float] = None,
) -> BinaryStackReader:
    """Create a BinaryStackReader from a list of GDAL-readable files.

    Parameters
    ----------
    file_list : Sequence[Filename]
        List of paths to the files to read.
    band : int, optional
        Band to read from the file, by default 1
    num_threads : int, optional (default 1)
        Number of threads to use for reading.
    nodata : float, optional
        Manually set value to use for nodata pixels, by default None
        If None passed, will search for a nodata value in the file.

    Returns
    -------
    BinaryStackReader
        The BinaryStackReader object.

    """
    readers = []
    dtypes = set()
    shapes = set()
    for f in file_list:
        with rio.open(f) as src:
            dtypes.add(src.dtypes[band - 1])
            shapes.add(src.shape)
        if len(dtypes) > 1:
            msg = "All files must have the same data type."
            raise ValueError(msg)
        if len(shapes) > 1:
            msg = "All files must have the same shape."
            raise ValueError(msg)
        readers.append(BinaryReader.from_gdal(f, band=band))
    return cls(
        file_list=file_list,
        readers=readers,
        num_threads=num_threads,
        nodata=nodata,
    )

BlockIndices dataclass

Class holding slices for 2D array access.

Source code in src/dolphin/io/_blocks.py

@dataclass(frozen=True)
class BlockIndices:
    """Class holding slices for 2D array access."""

    row_start: int
    row_stop: Optional[int]  # Can be None if we want slice(0, None)
    col_start: int
    col_stop: Optional[int]

    @classmethod
    def from_slices(cls, row_slice: slice, col_slice: slice) -> BlockIndices:
        return cls(
            row_start=row_slice.start,
            row_stop=row_slice.stop,
            col_start=col_slice.start,
            col_stop=col_slice.stop,
        )

    @property
    def row_slice(self) -> slice:
        return slice(self.row_start, self.row_stop)

    @property
    def col_slice(self) -> slice:
        return slice(self.col_start, self.col_stop)

    # Allows user to unpack like (row_slice, col_slice) = block_index
    def __iter__(self):
        return iter((self.row_slice, self.col_slice))

BlockProcessor

Bases: Protocol

Protocol for a block-wise processing function.

Reads a block of data from each reader, processes it, and returns the result as an array-like object.

Source code in src/dolphin/io/_process.py

class BlockProcessor(Protocol):
    """Protocol for a block-wise processing function.

    Reads a block of data from each reader, processes it, and returns the result
    as an array-like object.
    """

    def __call__(
        self,
        readers: Sequence[StackReader],
        rows: slice,
        cols: slice,
    ) -> tuple[ArrayLike, slice, slice]: ...

DatasetReader

Bases: Protocol

An array-like interface for reading input datasets.

DatasetReader defines the abstract interface that types must conform to in order to be read by functions which iterate in blocks over the input data. Such objects must export NumPy-like dtype, shape, and ndim attributes, and must support NumPy-style slice-based indexing.

Note that this protocol allows objects to be passed to dask.array.from_array which needs .shape, .ndim, .dtype and support numpy-style slicing.

Source code in src/dolphin/io/_readers.py

@runtime_checkable
class DatasetReader(Protocol):
    """An array-like interface for reading input datasets.

    `DatasetReader` defines the abstract interface that types must conform to in order
    to be read by functions which iterate in blocks over the input data.
    Such objects must export NumPy-like `dtype`, `shape`, and `ndim` attributes,
    and must support NumPy-style slice-based indexing.

    Note that this protocol allows objects to be passed to `dask.array.from_array`
    which needs `.shape`, `.ndim`, `.dtype` and support numpy-style slicing.
    """

    dtype: np.dtype
    """numpy.dtype : Data-type of the array's elements."""

    shape: tuple[int, ...]
    """tuple of int : Tuple of array dimensions."""

    ndim: int
    """int : Number of array dimensions."""

    def __getitem__(self, key: tuple[Index, ...], /) -> ArrayLike:
        """Read a block of data."""
        ...

dtype instance-attribute

numpy.dtype : Data-type of the array's elements.

ndim instance-attribute

int : Number of array dimensions.

shape instance-attribute

tuple of int : Tuple of array dimensions.

__getitem__(key)

Read a block of data.

Source code in src/dolphin/io/_readers.py

def __getitem__(self, key: tuple[Index, ...], /) -> ArrayLike:
    """Read a block of data."""
    ...

DatasetStackWriter

Bases: Protocol

An array-like interface for writing output datasets.

DatasetWriter defines the abstract interface that types must conform to in order to be used by functions which write outputs in blocks. Such objects must export NumPy-like dtype, shape, and ndim attributes, and must support NumPy-style slice-based indexing for setting data..

Source code in src/dolphin/io/_writers.py

@runtime_checkable
class DatasetStackWriter(Protocol):
    """An array-like interface for writing output datasets.

    `DatasetWriter` defines the abstract interface that types must conform to in order
    to be used by functions which write outputs in blocks.
    Such objects must export NumPy-like `dtype`, `shape`, and `ndim` attributes,
    and must support NumPy-style slice-based indexing for setting data..
    """

    dtype: np.dtype
    """numpy.dtype : Data-type of the array's elements."""

    shape: tuple[int, ...]
    """tuple of int : Tuple of array dimensions."""

    ndim: int = 3
    """int : Number of array dimensions."""

    def __setitem__(self, key: tuple[Index, ...], value: np.ndarray, /) -> None:
        """Write a block of data."""
        ...

dtype instance-attribute

numpy.dtype : Data-type of the array's elements.

ndim = 3 class-attribute instance-attribute

int : Number of array dimensions.

shape instance-attribute

tuple of int : Tuple of array dimensions.

__setitem__(key, value)

Write a block of data.

Source code in src/dolphin/io/_writers.py

def __setitem__(self, key: tuple[Index, ...], value: np.ndarray, /) -> None:
    """Write a block of data."""
    ...

DatasetWriter

Bases: Protocol

An array-like interface for writing output datasets.

DatasetWriter defines the abstract interface that types must conform to in order to be used by functions which write outputs in blocks. Such objects must export NumPy-like dtype, shape, and ndim attributes, and must support NumPy-style slice-based indexing for setting data..

Source code in src/dolphin/io/_writers.py

@runtime_checkable
class DatasetWriter(Protocol):
    """An array-like interface for writing output datasets.

    `DatasetWriter` defines the abstract interface that types must conform to in order
    to be used by functions which write outputs in blocks.
    Such objects must export NumPy-like `dtype`, `shape`, and `ndim` attributes,
    and must support NumPy-style slice-based indexing for setting data..
    """

    dtype: np.dtype
    """numpy.dtype : Data-type of the array's elements."""

    shape: tuple[int, ...]
    """tuple of int : Tuple of array dimensions."""

    ndim: int
    """int : Number of array dimensions."""

    def __setitem__(self, key: tuple[Index, ...], value: np.ndarray, /) -> None:
        """Write a block of data."""
        ...

dtype instance-attribute

numpy.dtype : Data-type of the array's elements.

ndim instance-attribute

int : Number of array dimensions.

shape instance-attribute

tuple of int : Tuple of array dimensions.

__setitem__(key, value)

Write a block of data.

Source code in src/dolphin/io/_writers.py

def __setitem__(self, key: tuple[Index, ...], value: np.ndarray, /) -> None:
    """Write a block of data."""
    ...

EagerLoader

Bases: BackgroundReader

Class to pre-fetch data chunks in a background thread.

Source code in src/dolphin/io/_readers.py

class EagerLoader(BackgroundReader):
    """Class to pre-fetch data chunks in a background thread."""

    def __init__(
        self,
        reader: DatasetReader,
        block_shape: tuple[int, int],
        overlaps: tuple[int, int] = (0, 0),
        skip_empty: bool = True,
        nodata_value: Optional[float] = None,
        nodata_mask: Optional[ArrayLike] = None,
        queue_size: int = 1,
        timeout: float = _DEFAULT_TIMEOUT,
    ):
        super().__init__(nq=queue_size, timeout=timeout, name="EagerLoader")
        self.reader = reader
        # Set up the generator of ((row_start, row_end), (col_start, col_end))
        # convert the slice generator to a list so we have the size
        nrows, ncols = self.reader.shape[-2:]
        self.slices = list(
            iter_blocks(
                arr_shape=(nrows, ncols),
                block_shape=block_shape,
                overlaps=overlaps,
            )
        )
        if nodata_value is None:
            nodata_value = getattr(reader, "nodata", None)
        self._queue_size = queue_size
        self._skip_empty = skip_empty
        self._nodata_mask = nodata_mask
        self._block_shape = block_shape
        self._nodata = nodata_value
        if self._nodata is None:
            self._nodata = np.nan

    def read(self, rows: slice, cols: slice) -> tuple[np.ndarray, tuple[slice, slice]]:
        logger.debug(f"EagerLoader reading {rows}, {cols}")
        cur_block = self.reader[..., rows, cols]
        return cur_block, (rows, cols)

    def iter_blocks(
        self, **tqdm_kwargs
    ) -> Generator[tuple[np.ndarray, tuple[slice, slice]], None, None]:
        # Queue up all slices to the work queue
        queued_slices = []
        for rows, cols in self.slices:
            # Skip queueing a read if all nodata
            if self._skip_empty and self._nodata_mask is not None:
                logger.debug("Checking nodata mask")
                if self._nodata_mask[rows, cols].all():
                    logger.debug("Skipping!")
                    continue
            self.queue_read(rows, cols)
            queued_slices.append((rows, cols))

        logger.info(f"Processing {self._block_shape} sized blocks...")
        for _ in trange(len(queued_slices), **tqdm_kwargs):
            cur_block, (rows, cols) = self.get_data()
            logger.debug(f"got data for {rows, cols}: {cur_block.shape}")

            # Otherwise look at the actual block we loaded
            if self._skip_empty:
                if isinstance(cur_block, np.ma.MaskedArray) and cur_block.mask.all():
                    continue
                if np.isnan(self._nodata):
                    block_is_nodata = np.isnan(cur_block)
                else:
                    block_is_nodata = cur_block == self._nodata
                if np.all(block_is_nodata):
                    logger.debug(
                        f"Skipping block {rows}, {cols} since it was all nodata"
                    )
                    continue
            yield cur_block, (rows, cols)

        self.notify_finished()

HDF5Reader dataclass

Bases: DatasetReader

A Dataset in an HDF5 file.

Attributes:

Name	Type	Description
filename	Path | str	Location of HDF5 file.
dset_name	str	Path to the dataset within the file.
chunks	(tuple[int, ...], optional)	Chunk shape of the dataset, or None if file is unchunked.
keep_open	bool, optional (default False)	If True, keep the HDF5 file handle open for faster reading.

See Also

BinaryReader RasterReader

Notes

If keep_open=True, this class does not store an open file object. Otherwise, the file is opened on-demand for reading or writing and closed immediately after each read/write operation. If passing the HDF5Reader to multiple spawned processes, it is recommended to set keep_open=False .

Source code in src/dolphin/io/_readers.py

@dataclass
class HDF5Reader(DatasetReader):
    """A Dataset in an HDF5 file.

    Attributes
    ----------
    filename : pathlib.Path | str
        Location of HDF5 file.
    dset_name : str
        Path to the dataset within the file.
    chunks : tuple[int, ...], optional
        Chunk shape of the dataset, or None if file is unchunked.
    keep_open : bool, optional (default False)
        If True, keep the HDF5 file handle open for faster reading.


    See Also
    --------
    BinaryReader
    RasterReader

    Notes
    -----
    If `keep_open=True`, this class does not store an open file object.
    Otherwise, the file is opened on-demand for reading or writing and closed
    immediately after each read/write operation.
    If passing the `HDF5Reader` to multiple spawned processes, it is recommended
    to set `keep_open=False` .

    """

    filename: Path
    """pathlib.Path : The file path."""

    dset_name: str
    """str : The path to the dataset within the file."""

    nodata: Optional[float] = None
    """Optional[float] : Value to use for nodata pixels.

    If None, looks for `_FillValue` or `missing_value` attributes on the dataset.
    """

    keep_open: bool = False
    """bool : If True, keep the HDF5 file handle open for faster reading."""

    keepdims: bool = True
    """bool : Maintain the dimension of the point array. If set to False, will
    skip `squeeze` on outputs with one dimension size of 1. Default is True."""

    def __post_init__(self):
        filename = Path(self.filename)

        hf = h5py.File(filename, "r")
        close_on_exit = True
        try:
            dset = hf[self.dset_name]
            self.shape = dset.shape
            self.dtype = dset.dtype
            self.chunks = dset.chunks
            if self.nodata is None:
                self.nodata = dset.attrs.get("_FillValue", None)
                if self.nodata is None:
                    self.nodata = dset.attrs.get("missing_value", None)
            if self.keep_open:
                self._hf = hf
                self._dset = dset
                close_on_exit = False
        finally:
            if close_on_exit:
                hf.close()

    def close(self):
        """Close the underlying HDF5 file handle if it was kept open."""
        if self.keep_open and hasattr(self, "_hf"):
            self._hf.close()

    def __del__(self):
        self.close()

    def __enter__(self):
        return self

    def __exit__(self, *args):
        self.close()

    @property
    def ndim(self) -> int:  # type: ignore[override]
        """Int : Number of array dimensions."""
        return len(self.shape)

    def __array__(self) -> np.ndarray:
        return self[:,]

    def __getitem__(self, key: tuple[Index, ...], /) -> np.ndarray:
        if self.keep_open:
            data = self._dset[key]
        else:
            with h5py.File(self.filename, "r") as f:
                data = f[self.dset_name][key]
        return _mask_array(data, self.nodata) if self.nodata is not None else data

dset_name instance-attribute

str : The path to the dataset within the file.

filename instance-attribute

pathlib.Path : The file path.

keep_open = False class-attribute instance-attribute

bool : If True, keep the HDF5 file handle open for faster reading.

keepdims = True class-attribute instance-attribute

bool : Maintain the dimension of the point array. If set to False, will skip squeeze on outputs with one dimension size of 1. Default is True.

ndim property

Int : Number of array dimensions.

nodata = None class-attribute instance-attribute

Optional[float] : Value to use for nodata pixels.

If None, looks for _FillValue or missing_value attributes on the dataset.

close()

Close the underlying HDF5 file handle if it was kept open.

Source code in src/dolphin/io/_readers.py

def close(self):
    """Close the underlying HDF5 file handle if it was kept open."""
    if self.keep_open and hasattr(self, "_hf"):
        self._hf.close()

HDF5StackReader dataclass

Bases: BaseStackReader

A stack of datasets in an HDF5 file.

See Also

BinaryStackReader StackReader

Notes

If keep_open=True, this class stores an open file object. Otherwise, the file is opened on-demand for reading or writing and closed immediately after each read/write operation. If passing the HDF5StackReader to multiple spawned processes, it is recommended to set keep_open=False.

Source code in src/dolphin/io/_readers.py

@dataclass
class HDF5StackReader(BaseStackReader):
    """A stack of datasets in an HDF5 file.

    See Also
    --------
    BinaryStackReader
    StackReader

    Notes
    -----
    If `keep_open=True`, this class stores an open file object.
    Otherwise, the file is opened on-demand for reading or writing and closed
    immediately after each read/write operation.
    If passing the `HDF5StackReader` to multiple spawned processes, it is recommended
    to set `keep_open=False`.

    """

    @classmethod
    def from_file_list(
        cls,
        file_list: Iterable[Filename],
        dset_names: str | Sequence[str],
        keep_open: bool = False,
        num_threads: int = 1,
        nodata: Optional[float] = None,
    ) -> HDF5StackReader:
        """Create a HDF5StackReader from a list of files.

        Parameters
        ----------
        file_list : Iterable[Filename]
            Iterable of paths to the files to read.
        dset_names : str | Sequence[str]
            Name of the dataset to read from each file.
            If a single string, will be used for all files.
        keep_open : bool, optional (default False)
            If True, keep the HDF5 file handles open for faster reading.
        num_threads : int, optional (default 1)
            Number of threads to use for reading.
        nodata : float, optional
            Manually set value to use for nodata pixels, by default None
            If None passed, will search for a nodata value in the file.

        Returns
        -------
        HDF5StackReader
            The HDF5StackReader object.

        """
        files = list(file_list)
        if isinstance(dset_names, str):
            dset_names = [dset_names] * len(files)

        readers = [
            HDF5Reader(Path(f), dset_name=dn, keep_open=keep_open, nodata=nodata)
            for (f, dn) in zip(files, dset_names, strict=False)
        ]
        # Check if nodata values were found in the files
        nds = {r.nodata for r in readers}
        if len(nds) == 1:
            nodata = nds.pop()

        return cls(files, readers, num_threads=num_threads, nodata=nodata)

from_file_list(file_list, dset_names, keep_open=False, num_threads=1, nodata=None) classmethod

Create a HDF5StackReader from a list of files.

Parameters:

Name	Type	Description	Default
file_list	Iterable[Filename]	Iterable of paths to the files to read.	required
dset_names	str | Sequence[str]	Name of the dataset to read from each file. If a single string, will be used for all files.	required
keep_open	bool, optional (default False)	If True, keep the HDF5 file handles open for faster reading.	False
num_threads	int, optional (default 1)	Number of threads to use for reading.	1
nodata	float	Manually set value to use for nodata pixels, by default None If None passed, will search for a nodata value in the file.	None

Returns:

Type	Description
HDF5StackReader	The HDF5StackReader object.

Source code in src/dolphin/io/_readers.py

@classmethod
def from_file_list(
    cls,
    file_list: Iterable[Filename],
    dset_names: str | Sequence[str],
    keep_open: bool = False,
    num_threads: int = 1,
    nodata: Optional[float] = None,
) -> HDF5StackReader:
    """Create a HDF5StackReader from a list of files.

    Parameters
    ----------
    file_list : Iterable[Filename]
        Iterable of paths to the files to read.
    dset_names : str | Sequence[str]
        Name of the dataset to read from each file.
        If a single string, will be used for all files.
    keep_open : bool, optional (default False)
        If True, keep the HDF5 file handles open for faster reading.
    num_threads : int, optional (default 1)
        Number of threads to use for reading.
    nodata : float, optional
        Manually set value to use for nodata pixels, by default None
        If None passed, will search for a nodata value in the file.

    Returns
    -------
    HDF5StackReader
        The HDF5StackReader object.

    """
    files = list(file_list)
    if isinstance(dset_names, str):
        dset_names = [dset_names] * len(files)

    readers = [
        HDF5Reader(Path(f), dset_name=dn, keep_open=keep_open, nodata=nodata)
        for (f, dn) in zip(files, dset_names, strict=False)
    ]
    # Check if nodata values were found in the files
    nds = {r.nodata for r in readers}
    if len(nds) == 1:
        nodata = nds.pop()

    return cls(files, readers, num_threads=num_threads, nodata=nodata)

RasterReader dataclass

Bases: DatasetReader

A single raster band of a GDAL-compatible dataset.

See Also

BinaryReader HDF5

Notes

If keep_open=True, this class does not store an open file object. Otherwise, the file is opened on-demand for reading or writing and closed immediately after each read/write operation. If passing the RasterReader to multiple spawned processes, it is recommended to set keep_open=False .

Source code in src/dolphin/io/_readers.py

@dataclass
class RasterReader(DatasetReader):
    """A single raster band of a GDAL-compatible dataset.

    See Also
    --------
    BinaryReader
    HDF5

    Notes
    -----
    If `keep_open=True`, this class does not store an open file object.
    Otherwise, the file is opened on-demand for reading or writing and closed
    immediately after each read/write operation.
    If passing the `RasterReader` to multiple spawned processes, it is recommended
    to set `keep_open=False` .

    """

    filename: Filename
    """Filename : The file path."""

    band: int
    """int : Band index (1-based)."""

    driver: str
    """str : Raster format driver name."""

    crs: rio.crs.CRS
    """rio.crs.CRS : The dataset's coordinate reference system."""

    transform: rio.transform.Affine
    """
    rasterio.transform.Affine : The dataset's georeferencing transformation matrix.

    This transform maps pixel row/column coordinates to coordinates in the dataset's
    coordinate reference system.
    """

    shape: tuple[int, int]
    dtype: np.dtype

    nodata: Optional[float] = None
    """Optional[float] : Value to use for nodata pixels."""

    keep_open: bool = False
    """bool : If True, keep the rasterio file handle open for faster reading."""

    chunks: Optional[tuple[int, int]] = None
    """Optional[tuple[int, int]] : Chunk shape of the dataset, or None if unchunked."""

    keepdims: bool = True
    """bool : Maintain the dimension of the point array. If set to False, will
    skip `squeeze` on outputs with one dimension size of 1. Default is True."""

    @classmethod
    def from_file(
        cls,
        filename: Filename,
        band: int = 1,
        nodata: Optional[float] = None,
        keepdims: bool = True,
        keep_open: bool = False,
        **options,
    ) -> RasterReader:
        with rio.open(filename, "r", **options) as src:
            shape = (src.height, src.width)
            dtype = np.dtype(src.dtypes[band - 1])
            driver = src.driver
            crs = src.crs
            nodata = nodata or src.nodatavals[band - 1]
            transform = src.transform
            chunks = src.block_shapes[band - 1]

            return cls(
                filename=filename,
                band=band,
                driver=driver,
                crs=crs,
                transform=transform,
                shape=shape,
                dtype=dtype,
                nodata=nodata,
                keepdims=keepdims,
                keep_open=keep_open,
                chunks=chunks,
            )

    def __post_init__(self):
        if self.keep_open:
            self._src = rio.open(self.filename, "r")

    def close(self):
        """Close the underlying rasterio dataset if it was kept open."""
        if self.keep_open and hasattr(self, "_src"):
            self._src.close()

    def __del__(self):
        self.close()

    def __enter__(self):
        return self

    def __exit__(self, *args):
        self.close()

    @property
    def ndim(self) -> int:  # type: ignore[override]
        """Int : Number of array dimensions."""
        return 2

    def __array__(self) -> np.ndarray:
        return self[:, :]

    def __getitem__(self, key: tuple[Index, ...], /) -> np.ndarray:
        if key is ... or key == ():
            key = (slice(None), slice(None))

        if not isinstance(key, tuple):
            msg = "Index must be a tuple of slices or integers."
            raise TypeError(msg)

        r_slice, c_slice = _ensure_slices(*key[-2:])
        window = rasterio.windows.Window.from_slices(
            r_slice,
            c_slice,
            height=self.shape[0],
            width=self.shape[1],
        )
        if self.keep_open:
            out = self._src.read(self.band, window=window)
        else:
            with rio.open(self.filename) as src:
                out = src.read(self.band, window=window)
        out_masked = _mask_array(out, self.nodata) if self.nodata is not None else out
        # Note that Rasterio doesn't use the `step` of a slice, so we need to
        # manually slice the output array.
        r_step, c_step = r_slice.step or 1, c_slice.step or 1
        if self.keepdims:
            return out_masked[::r_step, ::c_step]
        else:
            return out_masked[::r_step, ::c_step].squeeze()

band instance-attribute

int : Band index (1-based).

chunks = None class-attribute instance-attribute

Optional[tuple[int, int]] : Chunk shape of the dataset, or None if unchunked.

crs instance-attribute

rio.crs.CRS : The dataset's coordinate reference system.

driver instance-attribute

str : Raster format driver name.

filename instance-attribute

Filename : The file path.

keep_open = False class-attribute instance-attribute

bool : If True, keep the rasterio file handle open for faster reading.

keepdims = True class-attribute instance-attribute

bool : Maintain the dimension of the point array. If set to False, will skip squeeze on outputs with one dimension size of 1. Default is True.

ndim property

Int : Number of array dimensions.

nodata = None class-attribute instance-attribute

Optional[float] : Value to use for nodata pixels.

transform instance-attribute

rasterio.transform.Affine : The dataset's georeferencing transformation matrix.

This transform maps pixel row/column coordinates to coordinates in the dataset's coordinate reference system.

close()

Close the underlying rasterio dataset if it was kept open.

Source code in src/dolphin/io/_readers.py

def close(self):
    """Close the underlying rasterio dataset if it was kept open."""
    if self.keep_open and hasattr(self, "_src"):
        self._src.close()

RasterStackReader dataclass

Bases: BaseStackReader

A stack of datasets for any GDAL-readable rasters.

See Also

BinaryStackReader HDF5StackReader

Notes

If keep_open=True, this class stores an open file object. Otherwise, the file is opened on-demand for reading or writing and closed immediately after each read/write operation.

Source code in src/dolphin/io/_readers.py

@dataclass
class RasterStackReader(BaseStackReader):
    """A stack of datasets for any GDAL-readable rasters.

    See Also
    --------
    BinaryStackReader
    HDF5StackReader

    Notes
    -----
    If `keep_open=True`, this class stores an open file object.
    Otherwise, the file is opened on-demand for reading or writing and closed
    immediately after each read/write operation.

    """

    @classmethod
    def from_file_list(
        cls,
        file_list: Iterable[Filename],
        bands: int | Sequence[int] = 1,
        keepdims: bool = True,
        keep_open: bool = False,
        num_threads: int = 1,
        nodata: Optional[float] = None,
    ) -> RasterStackReader:
        """Create a RasterStackReader from a list of files.

        Parameters
        ----------
        file_list : Iterable[Filename]
            Iterable of paths to the files to read.
        bands : int | Sequence[int]
            Band to read from each file.
            If a single int, will be used for all files.
            Default = 1.
        keepdims : bool
            Maintain the dimension of the point array. If set to False, will
            skip `squeeze` on outputs with one dimension size of 1.
            Default is False.
        keep_open : bool, optional (default False)
            If True, keep the rasterio file handles open for faster reading.
        num_threads : int, optional (default 1)
            Number of threads to use for reading.
        nodata : float, optional
            Manually set value to use for nodata pixels, by default None

        Returns
        -------
        RasterStackReader
            The RasterStackReader object.

        """
        files = list(file_list)
        if isinstance(bands, int):
            bands = [bands] * len(files)

        readers = [
            RasterReader.from_file(f, band=b, keep_open=keep_open, keepdims=keepdims)
            for (f, b) in zip(files, bands, strict=False)
        ]
        # Check if nodata values were found in the files
        nds = {r.nodata for r in readers}
        if len(nds) == 1:
            nodata = nds.pop()
        return cls(
            files,
            readers,
            num_threads=num_threads,
            nodata=nodata,
            keepdims=keepdims,
        )

from_file_list(file_list, bands=1, keepdims=True, keep_open=False, num_threads=1, nodata=None) classmethod

Create a RasterStackReader from a list of files.

Parameters:

Name	Type	Description	Default
file_list	Iterable[Filename]	Iterable of paths to the files to read.	required
bands	int | Sequence[int]	Band to read from each file. If a single int, will be used for all files. Default = 1.	1
keepdims	bool	Maintain the dimension of the point array. If set to False, will skip squeeze on outputs with one dimension size of 1. Default is False.	True
keep_open	bool, optional (default False)	If True, keep the rasterio file handles open for faster reading.	False
num_threads	int, optional (default 1)	Number of threads to use for reading.	1
nodata	float	Manually set value to use for nodata pixels, by default None	None

Returns:

Type	Description
RasterStackReader	The RasterStackReader object.

Source code in src/dolphin/io/_readers.py

@classmethod
def from_file_list(
    cls,
    file_list: Iterable[Filename],
    bands: int | Sequence[int] = 1,
    keepdims: bool = True,
    keep_open: bool = False,
    num_threads: int = 1,
    nodata: Optional[float] = None,
) -> RasterStackReader:
    """Create a RasterStackReader from a list of files.

    Parameters
    ----------
    file_list : Iterable[Filename]
        Iterable of paths to the files to read.
    bands : int | Sequence[int]
        Band to read from each file.
        If a single int, will be used for all files.
        Default = 1.
    keepdims : bool
        Maintain the dimension of the point array. If set to False, will
        skip `squeeze` on outputs with one dimension size of 1.
        Default is False.
    keep_open : bool, optional (default False)
        If True, keep the rasterio file handles open for faster reading.
    num_threads : int, optional (default 1)
        Number of threads to use for reading.
    nodata : float, optional
        Manually set value to use for nodata pixels, by default None

    Returns
    -------
    RasterStackReader
        The RasterStackReader object.

    """
    files = list(file_list)
    if isinstance(bands, int):
        bands = [bands] * len(files)

    readers = [
        RasterReader.from_file(f, band=b, keep_open=keep_open, keepdims=keepdims)
        for (f, b) in zip(files, bands, strict=False)
    ]
    # Check if nodata values were found in the files
    nds = {r.nodata for r in readers}
    if len(nds) == 1:
        nodata = nds.pop()
    return cls(
        files,
        readers,
        num_threads=num_threads,
        nodata=nodata,
        keepdims=keepdims,
    )

RasterWriter dataclass

Bases: DatasetWriter, AbstractContextManager['RasterWriter']

A single raster band in a GDAL-compatible dataset containing one or more bands.

Raster provides a convenient interface for using SNAPHU to unwrap ground-projected interferograms in raster formats supported by the Geospatial Data Abstraction Library (GDAL). It acts as a thin wrapper around a Rasterio dataset and a band index, providing NumPy-like access to the underlying raster data.

Data access is performed lazily -- the raster contents are not stored in memory unless/until they are explicitly accessed by an indexing operation.

Raster objects must be closed after use in order to ensure that any data written to them is flushed to disk and any associated file objects are closed. The Raster class implements Python's context manager protocol, which can be used to reliably ensure that the raster is closed upon exiting the context manager.

Source code in src/dolphin/io/_writers.py

@dataclass
class RasterWriter(DatasetWriter, AbstractContextManager["RasterWriter"]):
    """A single raster band in a GDAL-compatible dataset containing one or more bands.

    `Raster` provides a convenient interface for using SNAPHU to unwrap ground-projected
    interferograms in raster formats supported by the Geospatial Data Abstraction
    Library (GDAL). It acts as a thin wrapper around a Rasterio dataset and a band
    index, providing NumPy-like access to the underlying raster data.

    Data access is performed lazily -- the raster contents are not stored in memory
    unless/until they are explicitly accessed by an indexing operation.

    `Raster` objects must be closed after use in order to ensure that any data written
    to them is flushed to disk and any associated file objects are closed. The `Raster`
    class implements Python's context manager protocol, which can be used to reliably
    ensure that the raster is closed upon exiting the context manager.
    """

    filename: Filename
    """str or Path : Path to the file to write."""
    band: int = 1
    """int : Band index in the file to write."""
    keep_bits: int | None = None
    """int : For floating point rasters, the number of mantissa bits to keep."""

    def __post_init__(self) -> None:
        # Open the dataset.
        self.dataset = rasterio.open(self.filename, mode="r+")

        # Check that `band` is a valid band index in the dataset.
        nbands = self.dataset.count
        if not (1 <= self.band <= nbands):
            errmsg = (
                f"band index {self.band} out of range: dataset contains {nbands} raster"
                " band(s)"
            )
            raise IndexError(errmsg)

        self.ndim = 2

    @classmethod
    def create(
        cls,
        fp: Filename,
        width: int | None = None,
        height: int | None = None,
        dtype: DTypeLike | None = None,
        driver: str | None = None,
        crs: str | Mapping[str, str] | rasterio.crs.CRS | None = None,
        transform: rasterio.transform.Affine | None = None,
        *,
        like_filename: Filename | None = None,
        keep_bits: int | None = None,
        **kwargs: Any,
    ) -> Self:
        """Create a new single-band raster dataset.

        If another raster is passed via the `like` argument, the new dataset will
        inherit the shape, data-type, driver, coordinate reference system (CRS), and
        geotransform of the reference raster. Driver-specific dataset creation options
        such as chunk size and compression flags may also be inherited.

        All other arguments take precedence over `like` and may be used to override
        attributes of the reference raster when creating the new raster.

        Parameters
        ----------
        fp : str or path-like
            File system path or URL of the local or remote dataset.
        width, height : int or None, optional
            The numbers of columns and rows of the raster dataset. Required if `like` is
            not specified. Otherwise, if None, the new dataset is created with the same
            width/height as `like`. Defaults to None.
        dtype : data-type or None, optional
            Data-type of the raster dataset's elements. Must be convertible to a
            `numpy.dtype` object and must correspond to a valid GDAL datatype. Required
            if `like` is not specified. Otherwise, if None, the new dataset is created
            with the same data-type as `like`. Defaults to None.
        driver : str or None, optional
            Raster format driver name. If None, the method will attempt to infer the
            driver from the file extension. Defaults to None.
        crs : str, dict, rasterio.crs.CRS, or None; optional
            The coordinate reference system. If None, the CRS of `like` will be used, if
            available, otherwise the raster will not be georeferenced. Defaults to None.
        transform : rasterio.transform.Affine or None, optional
            Affine transformation mapping the pixel space to geographic space. If None,
            the geotransform of `like` will be used, if available, otherwise the default
            transform will be used. Defaults to None.
        like_filename : Raster or None, optional
            An optional reference raster. If not None, the new raster will be created
            with the same metadata (shape, data-type, driver, CRS/geotransform, etc) as
            the reference raster. All other arguments will override the corresponding
            attribute of the reference raster. Defaults to None.
        keep_bits : int, optional
            Number of bits to preserve in mantissa. Defaults to None.
            Lower numbers will truncate the mantissa more and enable more compression.
        **kwargs : dict, optional
            Additional driver-specific creation options passed to `rasterio.open`.

        """
        if like_filename is not None:
            with rasterio.open(like_filename) as dataset:
                kwargs = dataset.profile | kwargs

        if width is not None:
            kwargs["width"] = width
        if height is not None:
            kwargs["height"] = height
        if dtype is not None:
            kwargs["dtype"] = np.dtype(dtype)
        if driver is not None:
            kwargs["driver"] = driver
        if crs is not None:
            kwargs["crs"] = crs
        if transform is not None:
            kwargs["transform"] = transform

        # Always create a single-band dataset, even if `like` was part of a multi-band
        # dataset.
        kwargs["count"] = 1

        # Create the new single-band dataset.
        with rasterio.open(fp, mode="w+", **kwargs):
            pass

        return cls(fp, band=1, keep_bits=keep_bits)

    @property
    def dtype(self) -> np.dtype:
        return np.dtype(self.dataset.dtypes[self.band - 1])

    @property
    def height(self) -> int:
        """int : The number of rows in the raster."""  # noqa: D403
        return self.dataset.height  # type: ignore[no-any-return]

    @property
    def width(self) -> int:
        """int : The number of columns in the raster."""  # noqa: D403
        return self.dataset.width  # type: ignore[no-any-return]

    @property
    def shape(self):
        return self.height, self.width

    @property
    def closed(self) -> bool:
        """bool : True if the dataset is closed."""  # noqa: D403
        return self.dataset.closed  # type: ignore[no-any-return]

    def close(self) -> None:
        """Close the underlying dataset.

        Has no effect if the dataset is already closed.
        """
        if not self.closed:
            self.dataset.close()

    def __exit__(self, exc_type, exc_value, traceback):  # type: ignore[no-untyped-def]
        self.close()

    def _window_from_slices(self, key: slice | tuple[slice, ...]) -> Window:
        if isinstance(key, slice):
            row_slice = key
            col_slice = slice(None)
        else:
            row_slice, col_slice = key

        return Window.from_slices(
            row_slice, col_slice, height=self.height, width=self.width
        )

    def __repr__(self) -> str:
        clsname = type(self).__name__
        return f"{clsname}(dataset={self.dataset!r}, band={self.band!r})"

    def __setitem__(self, key: tuple[Index, ...], value: np.ndarray, /) -> None:
        if np.issubdtype(value.dtype, np.floating) and self.keep_bits is not None:
            round_mantissa(value, keep_bits=self.keep_bits)
        with rasterio.open(
            self.filename,
            "r+",
        ) as dataset:
            if len(key) == 2:
                rows, cols = key
            elif len(key) == 3:
                _, rows, cols = _unpack_3d_slices(key)
            else:
                raise ValueError(
                    f"Invalid key for {self.__class__!r}.__setitem__: {key!r}"
                )
            try:
                window = Window.from_slices(
                    rows,
                    cols,
                    height=dataset.height,
                    width=dataset.width,
                )
            except rasterio.errors.WindowError as e:
                raise ValueError(f"Error creating window: {key = }, {value = }") from e
            return dataset.write(value, self.band, window=window)

band = 1 class-attribute instance-attribute

int : Band index in the file to write.

closed property

bool : True if the dataset is closed.

filename instance-attribute

str or Path : Path to the file to write.

height property

int : The number of rows in the raster.

keep_bits = None class-attribute instance-attribute

int : For floating point rasters, the number of mantissa bits to keep.

width property

int : The number of columns in the raster.

close()

Close the underlying dataset.

Has no effect if the dataset is already closed.

Source code in src/dolphin/io/_writers.py

def close(self) -> None:
    """Close the underlying dataset.

    Has no effect if the dataset is already closed.
    """
    if not self.closed:
        self.dataset.close()

create(fp, width=None, height=None, dtype=None, driver=None, crs=None, transform=None, *, like_filename=None, keep_bits=None, **kwargs) classmethod

Create a new single-band raster dataset.

If another raster is passed via the like argument, the new dataset will inherit the shape, data-type, driver, coordinate reference system (CRS), and geotransform of the reference raster. Driver-specific dataset creation options such as chunk size and compression flags may also be inherited.

All other arguments take precedence over like and may be used to override attributes of the reference raster when creating the new raster.

Parameters:

Name	Type	Description	Default
fp	str or path - like	File system path or URL of the local or remote dataset.	required
width	int or None	The numbers of columns and rows of the raster dataset. Required if like is not specified. Otherwise, if None, the new dataset is created with the same width/height as like. Defaults to None.	None
height	int or None	The numbers of columns and rows of the raster dataset. Required if like is not specified. Otherwise, if None, the new dataset is created with the same width/height as like. Defaults to None.	None
dtype	data - type or None	Data-type of the raster dataset's elements. Must be convertible to a numpy.dtype object and must correspond to a valid GDAL datatype. Required if like is not specified. Otherwise, if None, the new dataset is created with the same data-type as like. Defaults to None.	None
driver	str or None	Raster format driver name. If None, the method will attempt to infer the driver from the file extension. Defaults to None.	None
crs	str, dict, rasterio.crs.CRS, or None; optional	The coordinate reference system. If None, the CRS of like will be used, if available, otherwise the raster will not be georeferenced. Defaults to None.	None
transform	Affine or None	Affine transformation mapping the pixel space to geographic space. If None, the geotransform of like will be used, if available, otherwise the default transform will be used. Defaults to None.	None
like_filename	Raster or None	An optional reference raster. If not None, the new raster will be created with the same metadata (shape, data-type, driver, CRS/geotransform, etc) as the reference raster. All other arguments will override the corresponding attribute of the reference raster. Defaults to None.	None
keep_bits	int	Number of bits to preserve in mantissa. Defaults to None. Lower numbers will truncate the mantissa more and enable more compression.	None
**kwargs	dict	Additional driver-specific creation options passed to rasterio.open.	{}

Source code in src/dolphin/io/_writers.py

@classmethod
def create(
    cls,
    fp: Filename,
    width: int | None = None,
    height: int | None = None,
    dtype: DTypeLike | None = None,
    driver: str | None = None,
    crs: str | Mapping[str, str] | rasterio.crs.CRS | None = None,
    transform: rasterio.transform.Affine | None = None,
    *,
    like_filename: Filename | None = None,
    keep_bits: int | None = None,
    **kwargs: Any,
) -> Self:
    """Create a new single-band raster dataset.

    If another raster is passed via the `like` argument, the new dataset will
    inherit the shape, data-type, driver, coordinate reference system (CRS), and
    geotransform of the reference raster. Driver-specific dataset creation options
    such as chunk size and compression flags may also be inherited.

    All other arguments take precedence over `like` and may be used to override
    attributes of the reference raster when creating the new raster.

    Parameters
    ----------
    fp : str or path-like
        File system path or URL of the local or remote dataset.
    width, height : int or None, optional
        The numbers of columns and rows of the raster dataset. Required if `like` is
        not specified. Otherwise, if None, the new dataset is created with the same
        width/height as `like`. Defaults to None.
    dtype : data-type or None, optional
        Data-type of the raster dataset's elements. Must be convertible to a
        `numpy.dtype` object and must correspond to a valid GDAL datatype. Required
        if `like` is not specified. Otherwise, if None, the new dataset is created
        with the same data-type as `like`. Defaults to None.
    driver : str or None, optional
        Raster format driver name. If None, the method will attempt to infer the
        driver from the file extension. Defaults to None.
    crs : str, dict, rasterio.crs.CRS, or None; optional
        The coordinate reference system. If None, the CRS of `like` will be used, if
        available, otherwise the raster will not be georeferenced. Defaults to None.
    transform : rasterio.transform.Affine or None, optional
        Affine transformation mapping the pixel space to geographic space. If None,
        the geotransform of `like` will be used, if available, otherwise the default
        transform will be used. Defaults to None.
    like_filename : Raster or None, optional
        An optional reference raster. If not None, the new raster will be created
        with the same metadata (shape, data-type, driver, CRS/geotransform, etc) as
        the reference raster. All other arguments will override the corresponding
        attribute of the reference raster. Defaults to None.
    keep_bits : int, optional
        Number of bits to preserve in mantissa. Defaults to None.
        Lower numbers will truncate the mantissa more and enable more compression.
    **kwargs : dict, optional
        Additional driver-specific creation options passed to `rasterio.open`.

    """
    if like_filename is not None:
        with rasterio.open(like_filename) as dataset:
            kwargs = dataset.profile | kwargs

    if width is not None:
        kwargs["width"] = width
    if height is not None:
        kwargs["height"] = height
    if dtype is not None:
        kwargs["dtype"] = np.dtype(dtype)
    if driver is not None:
        kwargs["driver"] = driver
    if crs is not None:
        kwargs["crs"] = crs
    if transform is not None:
        kwargs["transform"] = transform

    # Always create a single-band dataset, even if `like` was part of a multi-band
    # dataset.
    kwargs["count"] = 1

    # Create the new single-band dataset.
    with rasterio.open(fp, mode="w+", **kwargs):
        pass

    return cls(fp, band=1, keep_bits=keep_bits)

S3Path

Bases: GeneralPath

A convenience class to handle paths on S3.

This class relies on pathlib.Path for operations using urllib to parse the url.

If passing a url with a trailing slash, the slash will be preserved when converting back to string.

Note that pure path manipulation functions do not require boto3, but functions which interact with S3 (e.g. exists(), .read_text()) do.

Attributes:

Name	Type	Description
bucket	str	Name of bucket in the url
path	Path	The URL path after s3:///
key	str	Alias of path converted to a string

Examples:

>>> from dolphin.io import S3Path
>>> s3_path = S3Path("s3://bucket/path/to/file.txt")
>>> str(s3_path)
's3://bucket/path/to/file.txt'
>>> s3_path.parent
S3Path("s3://bucket/path/to/")
>>> str(s3_path.parent)
's3://bucket/path/to/'

Source code in src/dolphin/io/_paths.py

class S3Path(GeneralPath):
    """A convenience class to handle paths on S3.

    This class relies on `pathlib.Path` for operations using `urllib` to parse the url.

    If passing a url with a trailing slash, the slash will be preserved
    when converting back to string.

    Note that pure path manipulation functions do *not* require `boto3`,
    but functions which interact with S3 (e.g. `exists()`, `.read_text()`) do.

    Attributes
    ----------
    bucket : str
        Name of bucket in the url
    path : pathlib.Path
        The URL path after s3://<bucket>/
    key : str
        Alias of `path` converted to a string

    Examples
    --------
    >>> from dolphin.io import S3Path
    >>> s3_path = S3Path("s3://bucket/path/to/file.txt")
    >>> str(s3_path)
    's3://bucket/path/to/file.txt'
    >>> s3_path.parent
    S3Path("s3://bucket/path/to/")
    >>> str(s3_path.parent)
    's3://bucket/path/to/'

    """

    def __init__(self, s3_url: Union[str, "S3Path"], unsigned: bool = False):
        """Create an S3Path.

        Parameters
        ----------
        s3_url : str or S3Path
            The S3 url to parse.
        unsigned : bool, optional
            If True, disable signing requests to S3.

        """
        # Names come from the urllib.parse.ParseResult
        if isinstance(s3_url, S3Path):
            self._scheme: str = s3_url._scheme
            self._netloc: str = s3_url._netloc
            self.bucket: str = s3_url.bucket
            self.path: Path = s3_url.path
            self._trailing_slash: str = s3_url._trailing_slash
        else:
            parsed: ParseResult = urlparse(s3_url)
            self._scheme = parsed.scheme
            self._netloc = self.bucket = parsed.netloc
            self._parsed = parsed
            self.path = Path(parsed.path)
            self._trailing_slash = "/" if s3_url.endswith("/") else ""

        if self._scheme != "s3":
            raise ValueError(f"{s3_url} is not an S3 url")

        self._unsigned = unsigned

    @classmethod
    def from_bucket_key(cls, bucket: str, key: str):
        """Create a `S3Path` from the bucket name and key/prefix.

        Matches API of some Boto3 functions which use this format.

        Parameters
        ----------
        bucket : str
            Name of S3 bucket.
        key : str
            S3 url of path after the bucket.

        """
        return cls(f"s3://{bucket}/{key}")

    def get_path(self):
        # For S3 paths, we need to add the double slash and netloc back to the front
        return f"{self._scheme}://{self._netloc}{self.path.as_posix()}{self._trailing_slash}"

    @property
    def key(self) -> str:
        """Name of key/prefix within the bucket with leading slash removed."""
        return f"{str(self.path.as_posix()).lstrip('/')}{self._trailing_slash}"

    @property
    def parent(self):
        parent_path = self.path.parent
        # Since this is a parent, it will will always end in a slash
        if self._scheme == "s3":
            # For S3 paths, we need to add the scheme and netloc back to the front
            return S3Path(f"{self._scheme}://{self._netloc}{parent_path.as_posix()}/")
        else:
            # For local paths, we can just convert the path to a string
            return S3Path(str(parent_path) + "/")

    @property
    def suffix(self):
        return self.path.suffix

    def resolve(self) -> S3Path:
        """Resolve the path to an absolute path- S3 paths are always absolute."""
        return self

    def _get_client(self):
        import boto3
        from botocore import UNSIGNED
        from botocore.config import Config

        if self._unsigned:
            return boto3.client("s3", config=Config(signature_version=UNSIGNED))
        else:
            return boto3.client("s3")

    def exists(self) -> bool:
        """Whether this path exists on S3."""
        client = self._get_client()
        resp = client.list_objects_v2(
            Bucket=self.bucket,
            Prefix=self.key,
            MaxKeys=1,
        )
        return resp.get("KeyCount") == 1

    def read_text(self) -> str:
        """Download/read the S3 file as text."""
        return self._download_as_bytes().decode()

    def read_bytes(self) -> bytes:
        """Download/read the S3 file as bytes."""
        return self._download_as_bytes()

    def _download_as_bytes(self) -> bytes:
        """Download file to a `BytesIO` buffer to read as bytes."""
        client = self._get_client()

        bio = BytesIO()
        client.download_fileobj(self.bucket, self.key, bio)
        bio.seek(0)
        out = bio.read()
        bio.close()
        return out

    def __truediv__(self, other):
        new = copy.deepcopy(self)
        new.path = self.path / other
        new._trailing_slash = "/" if str(other).endswith("/") else ""
        return new

    def __eq__(self, other):
        if isinstance(other, S3Path):
            return self.get_path() == other.get_path()
        elif isinstance(other, str):
            return self.get_path() == other
        else:
            return False

    def __hash__(self):
        return hash(self.get_path())

    def __repr__(self):
        return f'S3Path("{self.get_path()}")'

    def __str__(self):
        return self.get_path()

    def to_gdal(self):
        """Convert this S3Path to a GDAL URL."""
        return f"/vsis3/{self.bucket}/{self.key}"

key property

Name of key/prefix within the bucket with leading slash removed.

__init__(s3_url, unsigned=False)

Create an S3Path.

Parameters:

Name	Type	Description	Default
s3_url	str or S3Path	The S3 url to parse.	required
unsigned	bool	If True, disable signing requests to S3.	False

Source code in src/dolphin/io/_paths.py

def __init__(self, s3_url: Union[str, "S3Path"], unsigned: bool = False):
    """Create an S3Path.

    Parameters
    ----------
    s3_url : str or S3Path
        The S3 url to parse.
    unsigned : bool, optional
        If True, disable signing requests to S3.

    """
    # Names come from the urllib.parse.ParseResult
    if isinstance(s3_url, S3Path):
        self._scheme: str = s3_url._scheme
        self._netloc: str = s3_url._netloc
        self.bucket: str = s3_url.bucket
        self.path: Path = s3_url.path
        self._trailing_slash: str = s3_url._trailing_slash
    else:
        parsed: ParseResult = urlparse(s3_url)
        self._scheme = parsed.scheme
        self._netloc = self.bucket = parsed.netloc
        self._parsed = parsed
        self.path = Path(parsed.path)
        self._trailing_slash = "/" if s3_url.endswith("/") else ""

    if self._scheme != "s3":
        raise ValueError(f"{s3_url} is not an S3 url")

    self._unsigned = unsigned

exists()

Whether this path exists on S3.

Source code in src/dolphin/io/_paths.py

def exists(self) -> bool:
    """Whether this path exists on S3."""
    client = self._get_client()
    resp = client.list_objects_v2(
        Bucket=self.bucket,
        Prefix=self.key,
        MaxKeys=1,
    )
    return resp.get("KeyCount") == 1

from_bucket_key(bucket, key) classmethod

Create a S3Path from the bucket name and key/prefix.

Matches API of some Boto3 functions which use this format.

Parameters:

Name	Type	Description	Default
bucket	str	Name of S3 bucket.	required
key	str	S3 url of path after the bucket.	required

Source code in src/dolphin/io/_paths.py

@classmethod
def from_bucket_key(cls, bucket: str, key: str):
    """Create a `S3Path` from the bucket name and key/prefix.

    Matches API of some Boto3 functions which use this format.

    Parameters
    ----------
    bucket : str
        Name of S3 bucket.
    key : str
        S3 url of path after the bucket.

    """
    return cls(f"s3://{bucket}/{key}")

read_bytes()

Download/read the S3 file as bytes.

Source code in src/dolphin/io/_paths.py

def read_bytes(self) -> bytes:
    """Download/read the S3 file as bytes."""
    return self._download_as_bytes()

read_text()

Download/read the S3 file as text.

Source code in src/dolphin/io/_paths.py

def read_text(self) -> str:
    """Download/read the S3 file as text."""
    return self._download_as_bytes().decode()

resolve()

Resolve the path to an absolute path- S3 paths are always absolute.

Source code in src/dolphin/io/_paths.py

def resolve(self) -> S3Path:
    """Resolve the path to an absolute path- S3 paths are always absolute."""
    return self

to_gdal()

Convert this S3Path to a GDAL URL.

Source code in src/dolphin/io/_paths.py

def to_gdal(self):
    """Convert this S3Path to a GDAL URL."""
    return f"/vsis3/{self.bucket}/{self.key}"

StackReader

Bases: DatasetReader, Protocol

An array-like interface for reading a 3D stack of input datasets.

StackReader defines the abstract interface that types must conform to in order to be valid inputs to be read in functions like dolphin.ps.create_ps. It is a specialization of [DatasetReader][] that requires a 3D shape.

Source code in src/dolphin/io/_readers.py

@runtime_checkable
class StackReader(DatasetReader, Protocol):
    """An array-like interface for reading a 3D stack of input datasets.

    `StackReader` defines the abstract interface that types must conform to in order
    to be valid inputs to be read in functions like [dolphin.ps.create_ps][].
    It is a specialization of [DatasetReader][] that requires a 3D shape.
    """

    ndim: int = 3
    """int : Number of array dimensions."""

    shape: tuple[int, int, int]
    """tuple of int : Tuple of array dimensions."""

    def __len__(self) -> int:
        """Int : Number of images in the stack."""
        return self.shape[0]

ndim = 3 class-attribute instance-attribute

int : Number of array dimensions.

shape instance-attribute

tuple of int : Tuple of array dimensions.

__len__()

Int : Number of images in the stack.

Source code in src/dolphin/io/_readers.py

def __len__(self) -> int:
    """Int : Number of images in the stack."""
    return self.shape[0]

StridedBlockManager dataclass

Class to handle slicing/trimming overlapping blocks with strides.

Source code in src/dolphin/io/_blocks.py

@dataclass
class StridedBlockManager:
    """Class to handle slicing/trimming overlapping blocks with strides."""

    arr_shape: tuple[int, int]
    """(row, col) of the full-res 2D image"""
    block_shape: tuple[int, int]
    """(row, col) size of each input block to operate on at one time"""
    strides: Strides = field(default_factory=lambda: Strides(1, 1))
    """Decimation/downsampling factor in y/row and x/column direction"""
    half_window: HalfWindow = field(default_factory=lambda: HalfWindow(0, 0))
    """For multi-looking iterations, size of the (full-res) half window
    in y/row and x/column direction.
    Used to find `overlaps` between blocks and `start_offset`/`end_margin` for
    `iter_blocks`."""

    def iter_blocks(
        self,
    ) -> Iterator[
        tuple[BlockIndices, BlockIndices, BlockIndices, BlockIndices, BlockIndices]
    ]:
        """Iterate over the input/output block indices.

        Yields
        ------
        output_block : BlockIndices
            The current slices for the output raster
        out_trim : BlockIndices
            Slices to use on a processed output block to remove nodata border pixels.
            These may be relative (e.g. slice(1, -1)), not absolute like `output_block`.
        input_block : BlockIndices
            Slices used to load the full-res input data
        input_no_padding : BlockIndices
            Slices which point to the position within the full-res data without padding
        input_trim : BlockIndices
            Slices to use on full-res input block to remove padding from half-window.

        """
        out_trim = self.get_trimming_block()
        for out_block in self.iter_outputs():
            input_no_padding = dilate_block(out_block, strides=self.strides)
            input_block = pad_block(input_no_padding, margins=self.input_padding_shape)
            input_trim = get_full_res_trim(input_no_padding, input_block)
            yield (out_block, out_trim, input_block, input_no_padding, input_trim)

    def _get_out_nodata_size(self, direction: str) -> int:
        half_win = getattr(self.half_window, direction)
        stride = getattr(self.strides, direction)
        return round(half_win / stride)

    @property
    def output_shape(self) -> tuple[int, int]:
        return compute_out_shape(self.arr_shape, self.strides)

    @property
    def out_block_shape(self) -> tuple[int, int]:
        return compute_out_shape(self.block_shape, self.strides)

    @property
    def input_padding_shape(self) -> tuple[int, int]:
        """Amount of extra padding the input blocks need.

        Depends on the window size, and the strides.
        """
        return (
            self.strides.y * self._get_out_nodata_size("y"),
            self.strides.x * self._get_out_nodata_size("x"),
        )

    @property
    def output_margin(self) -> tuple[int, int]:
        """The output margins that we ignore while iterating.

        Depends on the half window and strides.

        The `half_window` is in full-res (input) coordinates (which would be the
        amount to skip with no striding), so the output margin size is smaller
        """
        return (self._get_out_nodata_size("y"), self._get_out_nodata_size("x"))

    def iter_outputs(self) -> Iterator[BlockIndices]:
        yield from iter_blocks(
            arr_shape=self.output_shape,
            block_shape=self.out_block_shape,
            overlaps=(0, 0),  # We're not overlapping in the *output* grid
            start_offsets=self.output_margin,
            end_margin=self.output_margin,
        )

    def get_trimming_block(self) -> BlockIndices:
        """Compute the slices which trim output nodata values.

        When the BlockIndex gets dilated (using `strides`) and padded (using
        `half_window`), the result will have nodata around the edges.
        The size of the nodata pixels in the full-res block is just
            (half_window['y'], half_window['x'])
        In the output (strided) coordinates, the number of nodata pixels is
        shrunk by how many strides are taken.

        Note that this is independent of which block we're on; the number of
        nodata pixels on the border is always the same.
        """
        row_nodata = self._get_out_nodata_size("y")
        col_nodata = self._get_out_nodata_size("x")
        # Extra check if we have no trimming to do: use slice(0, None)
        row_end = -row_nodata if row_nodata > 0 else None
        col_end = -col_nodata if col_nodata > 0 else None
        return BlockIndices(row_nodata, row_end, col_nodata, col_end)

arr_shape instance-attribute

(row, col) of the full-res 2D image

block_shape instance-attribute

(row, col) size of each input block to operate on at one time

half_window = field(default_factory=(lambda: HalfWindow(0, 0))) class-attribute instance-attribute

For multi-looking iterations, size of the (full-res) half window in y/row and x/column direction. Used to find overlaps between blocks and start_offset/end_margin for iter_blocks.

input_padding_shape property

Amount of extra padding the input blocks need.

Depends on the window size, and the strides.

output_margin property

The output margins that we ignore while iterating.

Depends on the half window and strides.

The half_window is in full-res (input) coordinates (which would be the amount to skip with no striding), so the output margin size is smaller

strides = field(default_factory=(lambda: Strides(1, 1))) class-attribute instance-attribute

Decimation/downsampling factor in y/row and x/column direction

get_trimming_block()

Compute the slices which trim output nodata values.

When the BlockIndex gets dilated (using strides) and padded (using half_window), the result will have nodata around the edges. The size of the nodata pixels in the full-res block is just (half_window['y'], half_window['x']) In the output (strided) coordinates, the number of nodata pixels is shrunk by how many strides are taken.

Note that this is independent of which block we're on; the number of nodata pixels on the border is always the same.

Source code in src/dolphin/io/_blocks.py

def get_trimming_block(self) -> BlockIndices:
    """Compute the slices which trim output nodata values.

    When the BlockIndex gets dilated (using `strides`) and padded (using
    `half_window`), the result will have nodata around the edges.
    The size of the nodata pixels in the full-res block is just
        (half_window['y'], half_window['x'])
    In the output (strided) coordinates, the number of nodata pixels is
    shrunk by how many strides are taken.

    Note that this is independent of which block we're on; the number of
    nodata pixels on the border is always the same.
    """
    row_nodata = self._get_out_nodata_size("y")
    col_nodata = self._get_out_nodata_size("x")
    # Extra check if we have no trimming to do: use slice(0, None)
    row_end = -row_nodata if row_nodata > 0 else None
    col_end = -col_nodata if col_nodata > 0 else None
    return BlockIndices(row_nodata, row_end, col_nodata, col_end)

iter_blocks()

Iterate over the input/output block indices.

Yields:

Name	Type	Description
output_block	BlockIndices	The current slices for the output raster
out_trim	BlockIndices	Slices to use on a processed output block to remove nodata border pixels. These may be relative (e.g. slice(1, -1)), not absolute like output_block.
input_block	BlockIndices	Slices used to load the full-res input data
input_no_padding	BlockIndices	Slices which point to the position within the full-res data without padding
input_trim	BlockIndices	Slices to use on full-res input block to remove padding from half-window.

Source code in src/dolphin/io/_blocks.py

def iter_blocks(
    self,
) -> Iterator[
    tuple[BlockIndices, BlockIndices, BlockIndices, BlockIndices, BlockIndices]
]:
    """Iterate over the input/output block indices.

    Yields
    ------
    output_block : BlockIndices
        The current slices for the output raster
    out_trim : BlockIndices
        Slices to use on a processed output block to remove nodata border pixels.
        These may be relative (e.g. slice(1, -1)), not absolute like `output_block`.
    input_block : BlockIndices
        Slices used to load the full-res input data
    input_no_padding : BlockIndices
        Slices which point to the position within the full-res data without padding
    input_trim : BlockIndices
        Slices to use on full-res input block to remove padding from half-window.

    """
    out_trim = self.get_trimming_block()
    for out_block in self.iter_outputs():
        input_no_padding = dilate_block(out_block, strides=self.strides)
        input_block = pad_block(input_no_padding, margins=self.input_padding_shape)
        input_trim = get_full_res_trim(input_no_padding, input_block)
        yield (out_block, out_trim, input_block, input_no_padding, input_trim)

VRTStack

Bases: StackReader

Class for creating a virtual stack of raster files.

Attributes:

Name	Type	Description
file_list	list[Filename]	Paths or GDAL-compatible strings (NETCDF:...) for paths to files.
outfile	(Path, optional(default=Path('slc_stack.vrt')))	Name of output file to write
dates	list[list[DateOrDatetime]]	list, where each entry is all dates matched from the corresponding file in file_list. This is used to sort the files by date. Each entry is a list because some files (compressed SLCs) may have multiple dates in the filename.
use_abs_path	(bool, optional(default=True))	Write the filepaths of the SLCs in the VRT as "relative=0"
subdataset	(str, optional)	Subdataset to use from the files in file_list, if using NetCDF files.
sort_files	(bool, optional(default=True))	Sort the files in file_list. Assumes that the naming convention will sort the files in increasing time order.
file_date_fmt	(str, optional(default='%Y%m%d'))	Format string for parsing the dates from the filenames. Passed to [opera_utils.get_dates][].

Source code in src/dolphin/io/_readers.py

class VRTStack(StackReader):
    """Class for creating a virtual stack of raster files.

    Attributes
    ----------
    file_list : list[Filename]
        Paths or GDAL-compatible strings (NETCDF:...) for paths to files.
    outfile : pathlib.Path, optional (default = Path("slc_stack.vrt"))
        Name of output file to write
    dates : list[list[DateOrDatetime]]
        list, where each entry is all dates matched from the corresponding file
        in `file_list`. This is used to sort the files by date.
        Each entry is a list because some files (compressed SLCs) may have
        multiple dates in the filename.
    use_abs_path : bool, optional (default = True)
        Write the filepaths of the SLCs in the VRT as "relative=0"
    subdataset : str, optional
        Subdataset to use from the files in `file_list`, if using NetCDF files.
    sort_files : bool, optional (default = True)
        Sort the files in `file_list`. Assumes that the naming convention
        will sort the files in increasing time order.
    file_date_fmt : str, optional (default = "%Y%m%d")
        Format string for parsing the dates from the filenames.
        Passed to [opera_utils.get_dates][].

    """

    def __init__(
        self,
        file_list: Sequence[Filename],
        outfile: Filename = "slc_stack.vrt",
        use_abs_path: bool = True,
        subdataset: Optional[str] = None,
        sort_files: bool = True,
        file_date_fmt: str = "%Y%m%d",
        write_file: bool = True,
        fail_on_overwrite: bool = False,
        skip_size_check: bool = False,
        num_threads: int = 1,
        read_masked: bool = False,
    ):
        if Path(outfile).exists() and write_file:
            if fail_on_overwrite:
                msg = (
                    f"Output file {outfile} already exists. "
                    "Please delete or specify a different output file. "
                    "To create from an existing VRT, use the `from_vrt_file` method."
                )
                raise FileExistsError(msg)
            else:
                logger.debug(f"Overwriting {outfile}")

        self._use_abs_path = use_abs_path
        files: list[Filename | S3Path]
        if any(str(f).startswith("s3://") for f in file_list):
            files = [S3Path(str(f)) for f in file_list]
        elif use_abs_path:
            # Skip resolving remote URLs (https://, http://) — Path.resolve()
            # would prepend the cwd to them.
            files = [
                p if is_remote_url(p) else utils._resolve_gdal_path(p)
                for p in file_list
            ]
        else:
            files = list(file_list)

        # Sanitize remote URLs: Path() collapses `https://` to `https:/`,
        # which GDAL's HDF5 driver cannot resolve. Restore the missing slash.
        files = [_sanitize_remote_url(f) for f in files]

        # Extract the date/datetimes from the filenames
        dates = [get_dates(f, fmt=file_date_fmt) for f in file_list]
        if sort_files:
            files, dates = sort_files_by_date(files, file_date_fmt=file_date_fmt)

        # Save the attributes
        self.file_list = files
        self.dates = dates
        self.num_threads = num_threads
        self._read_masked = read_masked

        self.outfile = Path(outfile).resolve()
        # Assumes that all files use the same subdataset (if NetCDF)
        self.subdataset = subdataset

        if not skip_size_check:
            _assert_images_same_size(self._gdal_file_strings)

        # Use the first file in the stack to get size, transform info
        ds = gdal.Open(fspath(self._gdal_file_strings[0]))
        bnd1 = ds.GetRasterBand(1)
        self.xsize = ds.RasterXSize
        self.ysize = ds.RasterYSize
        self.nodatavals = []
        for i in range(1, ds.RasterCount + 1):
            bnd = ds.GetRasterBand(i)
            self.nodatavals.append(bnd.GetNoDataValue())
        self.nodata = self.nodatavals[0]
        # Should be CFloat32
        self.gdal_dtype = gdal.GetDataTypeName(bnd1.DataType)
        # Save these for setting at the end
        self.gt = ds.GetGeoTransform()
        self.proj = ds.GetProjection()
        self.srs = ds.GetSpatialRef()
        ds = bnd1 = None

        # GDAL's HDF5 driver doesn't read NISAR's CF grid_mapping, so the
        # metadata above is identity GT + empty projection. Anything that
        # rasterizes a real-world polygon onto this VRT (bounds_mask,
        # nodata_mask, stitching crop) then silently produces an all-zero
        # result. Patch from h5py for the NISAR case before writing the VRT.
        if (
            self.subdataset
            and io._core._is_nisar_h5(self.file_list[0])
            and self.gt == (0.0, 1.0, 0.0, 0.0, 0.0, 1.0)
            and not self.proj
        ):
            from osgeo import osr

            _, _, gt, _, wkt = io._core.read_nisar_grid_metadata(
                self.file_list[0], self.subdataset
            )
            self.gt = gt
            self.proj = wkt
            self.srs = osr.SpatialReference()
            self.srs.ImportFromWkt(wkt)
        # Save the subset info

        self.xoff, self.yoff = 0, 0
        self.xsize_sub, self.ysize_sub = self.xsize, self.ysize

        if write_file:
            self._write()

    def _write(self):
        """Write out the VRT file pointing to the stack of SLCs, erroring if exists."""
        with open(self.outfile, "w") as fid:
            fid.write(
                f'<VRTDataset rasterXSize="{self.xsize_sub}"'
                f' rasterYSize="{self.ysize_sub}">\n'
            )

            for idx, filename in enumerate(self._gdal_file_strings, start=1):
                chunk_size = io.get_raster_chunk_size(filename)
                # chunks in a vrt have a min of 16, max of 2**14=16384
                # https://github.com/OSGeo/gdal/blob/2530defa1e0052827bc98696e7806037a6fec86e/frmts/vrt/vrtrasterband.cpp#L339
                if any(b < 16 for b in chunk_size) or any(
                    b > 16384 for b in chunk_size
                ):
                    chunk_str = ""
                else:
                    chunk_str = (
                        f'blockXSize="{chunk_size[0]}" blockYSize="{chunk_size[1]}"'
                    )
                outstr = f"""  <VRTRasterBand dataType="{self.gdal_dtype}" band="{idx}" {chunk_str}>
    <SimpleSource>
      <SourceFilename>{filename}</SourceFilename>
      <SourceBand>1</SourceBand>
      <SrcRect xOff="{self.xoff}" yOff="{self.yoff}" xSize="{self.xsize_sub}" ySize="{self.ysize_sub}"/>
      <DstRect xOff="0" yOff="0" xSize="{self.xsize_sub}" ySize="{self.ysize_sub}"/>
    </SimpleSource>
  </VRTRasterBand>\n"""  # noqa: E501
                fid.write(outstr)

            fid.write("</VRTDataset>")

        # Set the georeferencing metadata
        ds = gdal.Open(fspath(self.outfile), gdal.GA_Update)
        ds.SetGeoTransform(self.gt)
        ds.SetProjection(self.proj)
        ds.SetSpatialRef(self.srs)
        if self.nodata is not None:
            for i in range(ds.RasterCount):
                # ds.GetRasterBand(i + 1).SetNoDataValue(self.nodatavals[i])
                # Force to be the same nodataval for all bands
                ds.GetRasterBand(i + 1).SetNoDataValue(self.nodata)

        ds = None

    @property
    def _gdal_file_strings(self) -> list[str]:
        """Get the GDAL-compatible paths to write to the VRT.

        If we're not using .h5 or .nc, this will just be the file_list as is.
        """
        out = []
        for f in self.file_list:
            if isinstance(f, S3Path):
                out.append(f.to_gdal())
            else:
                out.append(io.format_nc_filename(f, self.subdataset))
        return out

    def __fspath__(self):
        # Allows os.fspath() to work on the object, enabling rasterio.open()
        return fspath(self.outfile)

    @classmethod
    def from_vrt_file(cls, vrt_file, new_outfile=None, **kwargs):
        """Create a new VRTStack using an existing VRT file."""
        file_list, subdataset = _parse_vrt_file(vrt_file)
        if new_outfile is None:
            # Point to the same, if none provided
            new_outfile = vrt_file

        return cls(
            file_list,
            outfile=new_outfile,
            subdataset=subdataset,
            write_file=False,
            **kwargs,
        )

    @property
    def shape(self):
        """Get the 3D shape of the stack."""
        xsize, ysize = io.get_raster_xysize(self._gdal_file_strings[0])
        return (len(self.file_list), ysize, xsize)

    def __len__(self):
        return len(self.file_list)

    def __repr__(self):
        outname = fspath(self.outfile) if self.outfile else "(not written)"
        return f"VRTStack({len(self.file_list)} bands, outfile={outname})"

    def __eq__(self, other):
        if not isinstance(other, VRTStack):
            return False
        return (
            self._gdal_file_strings == other._gdal_file_strings
            and self.outfile == other.outfile
        )

    def __hash__(self):
        return hash(self._gdal_file_strings)

    @property
    def ndim(self):
        return 3

    @property
    def dtype(self):
        return io.get_raster_dtype(self._gdal_file_strings[0])

    def __getitem__(self, index):
        if isinstance(index, int):
            if index < 0:
                index = len(self) + index
            return self.read_stack(band=index + 1)

        # TODO: raise an error if they try to skip like [::2, ::2]
        # or pass it to read_stack... but I don't think I need to support it.
        n, rows, cols = index
        if isinstance(rows, int):
            rows = slice(rows, rows + 1)
        if isinstance(cols, int):
            cols = slice(cols, cols + 1)
        if isinstance(n, int):
            if n < 0:
                n = len(self) + n
            return self.read_stack(band=n + 1, rows=rows, cols=cols)
        elif n is ...:
            n = slice(None)

        bands = list(range(1, 1 + len(self)))[n]
        if len(bands) == len(self):
            # This will use gdal's ds.ReadAsRaster, no iteration needed
            data = self.read_stack(band=None, rows=rows, cols=cols)
        else:
            # Get only the bands we need
            if self.num_threads == 1:
                data = np.stack(
                    [self.read_stack(band=i, rows=rows, cols=cols) for i in bands],
                    axis=0,
                )
            else:
                with ThreadPoolExecutor(max_workers=self.num_threads) as executor:
                    results = executor.map(
                        lambda i: self.read_stack(band=i, rows=rows, cols=cols), bands
                    )
                data = np.stack(list(results), axis=0)

        return data

    def read_stack(
        self,
        band: Optional[int] = None,
        subsample_factor: int = 1,
        rows: Optional[slice] = None,
        cols: Optional[slice] = None,
        masked: bool | None = None,
        keepdims: bool = True,
    ):
        """Read in the SLC stack."""
        if masked is None:
            masked = self._read_masked
        data = io.load_gdal(
            self.outfile,
            band=band,
            subsample_factor=subsample_factor,
            rows=rows,
            cols=cols,
            masked=masked,
        )
        # Check to get around gdal `ds.ReadAsArray()` squashing dimensions
        if len(self) == 1 and keepdims:
            # Add the front (1,) dimension which is missing for a single file
            data = data[None]
        return data

shape property

Get the 3D shape of the stack.

from_vrt_file(vrt_file, new_outfile=None, **kwargs) classmethod

Create a new VRTStack using an existing VRT file.

Source code in src/dolphin/io/_readers.py

@classmethod
def from_vrt_file(cls, vrt_file, new_outfile=None, **kwargs):
    """Create a new VRTStack using an existing VRT file."""
    file_list, subdataset = _parse_vrt_file(vrt_file)
    if new_outfile is None:
        # Point to the same, if none provided
        new_outfile = vrt_file

    return cls(
        file_list,
        outfile=new_outfile,
        subdataset=subdataset,
        write_file=False,
        **kwargs,
    )

read_stack(band=None, subsample_factor=1, rows=None, cols=None, masked=None, keepdims=True)

Read in the SLC stack.

Source code in src/dolphin/io/_readers.py

def read_stack(
    self,
    band: Optional[int] = None,
    subsample_factor: int = 1,
    rows: Optional[slice] = None,
    cols: Optional[slice] = None,
    masked: bool | None = None,
    keepdims: bool = True,
):
    """Read in the SLC stack."""
    if masked is None:
        masked = self._read_masked
    data = io.load_gdal(
        self.outfile,
        band=band,
        subsample_factor=subsample_factor,
        rows=rows,
        cols=cols,
        masked=masked,
    )
    # Check to get around gdal `ds.ReadAsArray()` squashing dimensions
    if len(self) == 1 and keepdims:
        # Add the front (1,) dimension which is missing for a single file
        data = data[None]
    return data

copy_projection(src_file, dst_file)

Copy projection/geotransform from src_file to dst_file.

Source code in src/dolphin/io/_core.py

def copy_projection(src_file: Filename, dst_file: Filename) -> None:
    """Copy projection/geotransform from `src_file` to `dst_file`."""
    ds_src = _get_gdal_ds(src_file)
    projection = ds_src.GetProjection()
    geotransform = ds_src.GetGeoTransform()
    nodata = ds_src.GetRasterBand(1).GetNoDataValue()

    if projection is None and geotransform is None:
        logger.info("No projection or geotransform found on file %s", input)
        return
    ds_dst = gdal.Open(fspath(dst_file), gdal.GA_Update)

    if geotransform is not None and geotransform != (0, 1, 0, 0, 0, 1):
        ds_dst.SetGeoTransform(geotransform)

    if projection is not None and projection != "":
        ds_dst.SetProjection(projection)

    if nodata is not None:
        ds_dst.GetRasterBand(1).SetNoDataValue(nodata)

    ds_src = ds_dst = None

format_nc_filename(filename, ds_name=None)

Format an HDF5/NetCDF filename with dataset for reading using GDAL.

If filename is already formatted, or if filename is not an HDF5/NetCDF file (based on the file extension), it is returned unchanged.

The driver prefix is chosen by filename:

• .nc → NETCDF:"file":"//ds"
• .h5 whose granule prefix is NISAR_ → HDF5:"file":"//ds" (NISAR raw HDF5; no CF metadata, and GDAL's NETCDF driver refuses it with "No such file or directory")
• every other .h5 → NETCDF:"file":"//ds" (OPERA CSLCs, COMPASS CSLCs + static_layers, etc. ship CF-1.8- compliant HDF5 and depend on the NETCDF driver to read the grid mapping; the bare HDF5 driver opens the data but reports an identity geotransform)

Parameters:

Name	Type	Description	Default
filename	str or PathLike	Filename to format.	required
ds_name	str	Dataset name to use. If not provided for a .h5 or .nc file, an error is raised.	None

Returns:

Type	Description
str	Formatted filename.

Raises:

Type	Description
ValueError	If ds_name is not provided for a .h5 or .nc file.

Source code in src/dolphin/io/_core.py

def format_nc_filename(filename: Filename, ds_name: Optional[str] = None) -> str:
    """Format an HDF5/NetCDF filename with dataset for reading using GDAL.

    If `filename` is already formatted, or if `filename` is not an HDF5/NetCDF
    file (based on the file extension), it is returned unchanged.

    The driver prefix is chosen by filename:

    - ``.nc``  → ``NETCDF:"file":"//ds"``
    - ``.h5`` whose granule prefix is ``NISAR_`` → ``HDF5:"file":"//ds"``
      (NISAR raw HDF5; no CF metadata, and GDAL's NETCDF driver refuses
      it with "No such file or directory")
    - every other ``.h5`` → ``NETCDF:"file":"//ds"``
      (OPERA CSLCs, COMPASS CSLCs + static_layers, etc. ship CF-1.8-
      compliant HDF5 and depend on the NETCDF driver to read the grid
      mapping; the bare HDF5 driver opens the data but reports an
      identity geotransform)

    Parameters
    ----------
    filename : str or PathLike
        Filename to format.
    ds_name : str, optional
        Dataset name to use. If not provided for a .h5 or .nc file, an error is raised.

    Returns
    -------
    str
        Formatted filename.

    Raises
    ------
    ValueError
        If `ds_name` is not provided for a .h5 or .nc file.

    """
    # If we've already formatted the filename, return it
    fname_clean = fspath(filename).lstrip('"').lstrip("'").rstrip('"').rstrip("'")
    if fname_clean.startswith(("NETCDF:", "HDF5:")):
        return fspath(filename)

    if not (fname_clean.endswith((".nc", ".h5"))):
        return fspath(filename)

    # Now we're definitely dealing with an HDF5/NetCDF file
    if ds_name is None:
        msg = "Must provide dataset name for HDF5/NetCDF files"
        raise ValueError(msg)

    basename = Path(fname_clean).name.upper()
    if fname_clean.endswith(".h5") and basename.startswith("NISAR_"):
        driver = "HDF5"
    else:
        driver = "NETCDF"
    return f'{driver}:"{filename}":"//{ds_name.lstrip("/")}"'

get_gtiff_options(max_error=None, compression_type='lzw', chunk_size=256, predictor=None, zlevel=1, use_16_bits=False)

Generate GTiff creation options for GDAL translate.

Parameters:

Name	Type	Description	Default
max_error	float	Maximum compression error.	None
compression_type	str	Compression type to use (default is "lzw").	'lzw'
chunk_size	int	Size of the chunks for blockxsize and blockysize (default is 256).	256
predictor	int or None	Predictor type to use (default is 3). Use None to omit the predictor.	None
zlevel	int or None	Compression level for the 'deflate' and 'zstd' compression types (default is 1). Use None to omit the zlevel.	1
use_16_bits	bool	If True, sets NBITS=16 to write float32 rasters at half precision. Default is False.	False

Returns:

Type	Description
dict[str, str] | list[str]	List of GTiff creation options formatted for GDAL (if gdal_format=True) Otherwise, a dict mapping option to value for rasterio.

Source code in src/dolphin/io/_utils.py

def get_gtiff_options(
    max_error: float | None = None,
    compression_type: str = "lzw",
    chunk_size: int = 256,
    predictor: int | None = None,
    zlevel: int | None = 1,
    use_16_bits: bool = False,
) -> dict[str, str]:
    """Generate GTiff creation options for GDAL translate.

    Parameters
    ----------
    max_error : float
        Maximum compression error.
    compression_type : str, optional
        Compression type to use (default is "lzw").
    chunk_size : int, optional
        Size of the chunks for blockxsize and blockysize (default is 256).
    predictor : int or None, optional
        Predictor type to use (default is 3). Use None to omit the predictor.
    zlevel : int or None, optional
        Compression level for the 'deflate' and 'zstd' compression types (default is 1).
        Use None to omit the zlevel.
    use_16_bits: bool
        If True, sets `NBITS=16` to write float32 rasters at half precision.
        Default is False.

    Returns
    -------
    dict[str, str] | list[str]
        List of GTiff creation options formatted for GDAL (if `gdal_format=True`)
        Otherwise, a dict mapping option to value for rasterio.

    """
    options = {
        "bigtiff": "yes",
        "tiled": "yes",
        "blockxsize": str(chunk_size),
        "blockysize": str(chunk_size),
        "compress": compression_type,
    }
    if zlevel is not None:
        options["zlevel"] = str(zlevel)
    if predictor is not None:
        options["predictor"] = str(predictor)
    if compression_type.lower().startswith("lerc") and max_error is not None:
        options["max_z_error"] = str(max_error)
    if use_16_bits:
        from ._core import _can_use_nbits16

        if _can_use_nbits16():
            options["nbits"] = "16"

    return options

get_raster_bounds(filename=None, ds=None)

Get the (left, bottom, right, top) bounds of the image.

Source code in src/dolphin/io/_core.py

def get_raster_bounds(
    filename: Optional[Filename] = None, ds: Optional[gdal.Dataset] = None
) -> Bbox:
    """Get the (left, bottom, right, top) bounds of the image."""
    if ds is None:
        if filename is None:
            msg = "Must provide either `filename` or `ds`"
            raise ValueError(msg)
        ds = _get_gdal_ds(filename)

    gt = ds.GetGeoTransform()
    xsize, ysize = ds.RasterXSize, ds.RasterYSize

    left, top = _apply_gt(gt=gt, x=0, y=0)
    right, bottom = _apply_gt(gt=gt, x=xsize, y=ysize)

    return Bbox(left, bottom, right, top)

get_raster_chunk_size(filename)

Get size the raster's chunks on disk.

This is called blockXsize, blockYsize by GDAL.

Source code in src/dolphin/io/_core.py

def get_raster_chunk_size(filename: Filename) -> list[int]:
    """Get size the raster's chunks on disk.

    This is called blockXsize, blockYsize by GDAL.
    """
    ds = _get_gdal_ds(filename)
    block_size = ds.GetRasterBand(1).GetBlockSize()
    for i in range(2, ds.RasterCount + 1):
        if block_size != ds.GetRasterBand(i).GetBlockSize():
            logger.warning(f"Warning: {filename} bands have different block shapes.")
            break
    return block_size

get_raster_crs(filename)

Get the CRS from a file.

Parameters:

Name	Type	Description	Default
filename	Filename	Path to the file to load.	required

Returns:

Type	Description
CRS	CRS.

Source code in src/dolphin/io/_core.py

def get_raster_crs(filename: Filename) -> CRS:
    """Get the CRS from a file.

    Parameters
    ----------
    filename : Filename
        Path to the file to load.

    Returns
    -------
    CRS
        CRS.

    """
    ds = _get_gdal_ds(filename)
    return CRS.from_wkt(ds.GetProjection())

get_raster_description(filename, band=1)

Get description of a raster band.

Parameters:

Name	Type	Description	Default
filename	Filename	Path to the file to load.	required
band	int	Band to get description for. Default is 1.	1

Source code in src/dolphin/io/_core.py

def get_raster_description(filename: Filename, band: int = 1):
    """Get description of a raster band.

    Parameters
    ----------
    filename : Filename
        Path to the file to load.
    band : int, optional
        Band to get description for. Default is 1.

    """
    ds = _get_gdal_ds(filename)
    bnd = ds.GetRasterBand(band)
    return bnd.GetDescription()

get_raster_driver(filename)

Get the GDAL driver ShortName from a file.

Parameters:

Name	Type	Description	Default
filename	Filename	Path to the file to load.	required

Returns:

Type	Description
str	Driver name.

Source code in src/dolphin/io/_core.py

def get_raster_driver(filename: Filename) -> str:
    """Get the GDAL driver `ShortName` from a file.

    Parameters
    ----------
    filename : Filename
        Path to the file to load.

    Returns
    -------
    str
        Driver name.

    """
    ds = _get_gdal_ds(filename)
    return ds.GetDriver().ShortName

get_raster_dtype(filename)

Get the data type from a file.

Parameters:

Name	Type	Description	Default
filename	Filename	Path to the file to load.	required

Returns:

Type	Description
dtype	Data type.

Source code in src/dolphin/io/_core.py

def get_raster_dtype(filename: Filename) -> np.dtype:
    """Get the data type from a file.

    Parameters
    ----------
    filename : Filename
        Path to the file to load.

    Returns
    -------
    np.dtype
        Data type.

    """
    ds = _get_gdal_ds(filename)
    return gdal_to_numpy_type(ds.GetRasterBand(1).DataType)

get_raster_gt(filename)

Get the geotransform from a file.

Parameters:

Name	Type	Description	Default
filename	Filename	Path to the file to load.	required

Returns:

Type	Description
List[float]	6 floats representing a GDAL Geotransform.

Source code in src/dolphin/io/_core.py

def get_raster_gt(filename: Filename) -> list[float]:
    """Get the geotransform from a file.

    Parameters
    ----------
    filename : Filename
        Path to the file to load.

    Returns
    -------
    List[float]
        6 floats representing a GDAL Geotransform.

    """
    ds = _get_gdal_ds(filename)
    return ds.GetGeoTransform()

get_raster_metadata(filename, domain='')

Get metadata from a raster file.

Parameters:

Name	Type	Description	Default
filename	Filename	Path to the file to load.	required
domain	str	Domain to get metadata for. Default is "" (all domains).	''

Returns:

Type	Description
dict	Dictionary of metadata.

Source code in src/dolphin/io/_core.py

def get_raster_metadata(filename: Filename, domain: str = ""):
    """Get metadata from a raster file.

    Parameters
    ----------
    filename : Filename
        Path to the file to load.
    domain : str, optional
        Domain to get metadata for. Default is "" (all domains).

    Returns
    -------
    dict
        Dictionary of metadata.

    """
    ds = _get_gdal_ds(filename)
    return ds.GetMetadata(domain)

get_raster_nodata(filename, band=1)

Get the nodata value from a file.

Parameters:

Name	Type	Description	Default
filename	Filename	Path to the file to load.	required
band	int	Band to get nodata value for, by default 1.	1

Returns:

Type	Description
Optional[float]	Nodata value, or None if not found.

Source code in src/dolphin/io/_core.py

def get_raster_nodata(filename: Filename, band: int = 1) -> Optional[float]:
    """Get the nodata value from a file.

    Parameters
    ----------
    filename : Filename
        Path to the file to load.
    band : int, optional
        Band to get nodata value for, by default 1.

    Returns
    -------
    Optional[float]
        Nodata value, or None if not found.

    """
    ds = _get_gdal_ds(filename)
    return ds.GetRasterBand(band).GetNoDataValue()

get_raster_units(filename, band=1)

Get units of a raster band.

Parameters:

Name	Type	Description	Default
filename	Filename	Path to the file to load.	required
band	int	Band to get units for. Default is 1.	1

Source code in src/dolphin/io/_core.py

def get_raster_units(filename: Filename, band: int = 1) -> str | None:
    """Get units of a raster band.

    Parameters
    ----------
    filename : Filename
        Path to the file to load.
    band : int
        Band to get units for.
        Default is 1.

    """
    ds = _get_gdal_ds(filename)
    bnd = ds.GetRasterBand(band)
    return bnd.GetUnitType() or None

get_raster_xysize(filename)

Get the xsize/ysize of a GDAL-readable raster.

Source code in src/dolphin/io/_core.py

def get_raster_xysize(filename: Filename) -> tuple[int, int]:
    """Get the xsize/ysize of a GDAL-readable raster."""
    ds = _get_gdal_ds(filename)
    xsize, ysize = ds.RasterXSize, ds.RasterYSize
    ds = None
    return xsize, ysize

iter_blocks(arr_shape, block_shape, overlaps=(0, 0), start_offsets=(0, 0), end_margin=(0, 0))

Create a generator to get indexes for accessing blocks of a raster.

Parameters:

Name	Type	Description	Default
arr_shape	tuple[int, int]	(num_rows, num_cols), full size of array to access	required
block_shape	tuple[int, int]	(height, width), size of blocks to load	required
overlaps	tuple[int, int]	(row_overlap, col_overlap), number of pixels to re-include from the previous block after sliding	= (0, 0)
start_offsets	tuple[int, int]	Offsets from top left to start reading from	= (0, 0)
end_margin	tuple[int, int]	Margin to avoid at the bottom/right of array	= (0, 0)

Yields:

Type	Description
BlockIndices	Iterator of BlockIndices, which can be unpacked into (slice(row_start, row_stop), slice(col_start, col_stop))

Examples:

>>> list(iter_blocks((180, 250), (100, 100)))
[BlockIndices(row_start=0, row_stop=100, col_start=0, col_stop=100), BlockIndices(row_start=0, row_stop=100, col_start=100, col_stop=200), BlockIndices(row_start=0, row_stop=100, col_start=200, col_stop=250), BlockIndices(row_start=100, row_stop=180, col_start=0, col_stop=100), BlockIndices(row_start=100, row_stop=180, col_start=100, col_stop=200), BlockIndices(row_start=100, row_stop=180, col_start=200, col_stop=250)]
>>> list(map(tuple, iter_blocks((180, 250), (100, 100), overlaps=(10, 10))))
[(slice(0, 100, None), slice(0, 100, None)), (slice(0, 100, None), slice(90, 190, None)), (slice(0, 100, None), slice(180, 250, None)), (slice(90, 180, None), slice(0, 100, None)), (slice(90, 180, None), slice(90, 190, None)), (slice(90, 180, None), slice(180, 250, None))]

Source code in src/dolphin/io/_blocks.py

def iter_blocks(
    arr_shape: tuple[int, int],
    block_shape: tuple[int, int],
    overlaps: tuple[int, int] = (0, 0),
    start_offsets: tuple[int, int] = (0, 0),
    end_margin: tuple[int, int] = (0, 0),
) -> Iterator[BlockIndices]:
    """Create a generator to get indexes for accessing blocks of a raster.

    Parameters
    ----------
    arr_shape : tuple[int, int]
        (num_rows, num_cols), full size of array to access
    block_shape : tuple[int, int]
        (height, width), size of blocks to load
    overlaps : tuple[int, int], default = (0, 0)
        (row_overlap, col_overlap), number of pixels to re-include from
        the previous block after sliding
    start_offsets : tuple[int, int], default = (0, 0)
        Offsets from top left to start reading from
    end_margin : tuple[int, int], default = (0, 0)
        Margin to avoid at the bottom/right of array

    Yields
    ------
    BlockIndices
        Iterator of BlockIndices, which can be unpacked into
        (slice(row_start, row_stop), slice(col_start, col_stop))

    Examples
    --------
        >>> list(iter_blocks((180, 250), (100, 100)))
        [BlockIndices(row_start=0, row_stop=100, col_start=0, col_stop=100), \
BlockIndices(row_start=0, row_stop=100, col_start=100, col_stop=200), \
BlockIndices(row_start=0, row_stop=100, col_start=200, col_stop=250), \
BlockIndices(row_start=100, row_stop=180, col_start=0, col_stop=100), \
BlockIndices(row_start=100, row_stop=180, col_start=100, col_stop=200), \
BlockIndices(row_start=100, row_stop=180, col_start=200, col_stop=250)]
        >>> list(map(tuple, iter_blocks((180, 250), (100, 100), overlaps=(10, 10))))
        [(slice(0, 100, None), slice(0, 100, None)), (slice(0, 100, None), \
slice(90, 190, None)), (slice(0, 100, None), slice(180, 250, None)), \
(slice(90, 180, None), slice(0, 100, None)), (slice(90, 180, None), \
slice(90, 190, None)), (slice(90, 180, None), slice(180, 250, None))]

    """
    total_rows, total_cols = arr_shape
    height, width = block_shape
    row_overlap, col_overlap = overlaps
    start_row_offset, start_col_offset = start_offsets
    last_row = total_rows - end_margin[0]
    last_col = total_cols - end_margin[1]

    if height is None:
        height = total_rows
    if width is None:
        width = total_cols

    # Check we're not moving backwards with the overlap:
    if row_overlap >= height and height != total_rows:
        msg = f"{row_overlap = } must be less than block height {height}"
        raise ValueError(msg)
    if col_overlap >= width and width != total_cols:
        msg = f"{col_overlap = } must be less than block width {width}"
        raise ValueError(msg)

    # Set up the iterating indices
    cur_row = start_row_offset
    cur_col = start_col_offset
    while cur_row < last_row:
        while cur_col < last_col:
            row_stop = min(cur_row + height, last_row)
            col_stop = min(cur_col + width, last_col)
            # yield (slice(cur_row, row_stop), slice(cur_col, col_stop))
            yield BlockIndices(cur_row, row_stop, cur_col, col_stop)

            cur_col += width
            if cur_col < last_col:  # don't bring back if already at edge
                cur_col -= col_overlap

        cur_row += height
        if cur_row < last_row:
            cur_row -= row_overlap
        cur_col = start_col_offset  # reset back to the starting offset

load_gdal(filename, *, band=None, subsample_factor=1, overview=None, rows=None, cols=None, masked=False)

Load a gdal file into a numpy array.

Parameters:

Name	Type	Description	Default
filename	str or Path	Path to the file to load.	required
band	int	Band to load. If None, load all bands as 3D array.	None
subsample_factor	int or tuple[int, int]	Subsample the data by this factor. Default is 1 (no subsampling). Uses nearest neighbor resampling.	1
overview	Optional[int]	If passed, will load an overview of the file. Raster must have existing overviews, or ValueError is raised.	None
rows	slice	Rows to load. Default is None (load all rows).	None
cols	slice	Columns to load. Default is None (load all columns).	None
masked	bool	If True, return a masked array using the raster's nodata value. Default is False.	False

Returns:

Name	Type	Description
arr	ndarray or MaskedArray	Array of shape (bands, y, x) or (y, x) if band is specified, where y = height // subsample_factor and x = width // subsample_factor.

Source code in src/dolphin/io/_core.py

def load_gdal(
    filename: Filename,
    *,
    band: Optional[int] = None,
    subsample_factor: Union[int, tuple[int, int]] = 1,
    overview: Optional[int] = None,
    rows: Optional[slice] = None,
    cols: Optional[slice] = None,
    masked: bool = False,
) -> np.ndarray | np.ma.MaskedArray:
    """Load a gdal file into a numpy array.

    Parameters
    ----------
    filename : str or Path
        Path to the file to load.
    band : int, optional
        Band to load. If None, load all bands as 3D array.
    subsample_factor : int or tuple[int, int], optional
        Subsample the data by this factor. Default is 1 (no subsampling).
        Uses nearest neighbor resampling.
    overview: int, optional
        If passed, will load an overview of the file.
        Raster must have existing overviews, or ValueError is raised.
    rows : slice, optional
        Rows to load. Default is None (load all rows).
    cols : slice, optional
        Columns to load. Default is None (load all columns).
    masked : bool, optional
        If True, return a masked array using the raster's `nodata` value.
        Default is False.

    Returns
    -------
    arr : np.ndarray or np.ma.MaskedArray
        Array of shape (bands, y, x) or (y, x) if `band` is specified,
        where y = height // subsample_factor and x = width // subsample_factor.

    """
    ds = _get_gdal_ds(filename)
    nrows, ncols = ds.RasterYSize, ds.RasterXSize

    if overview is not None:
        # We can handle the overviews most easily
        bnd = ds.GetRasterBand(band or 1)
        ovr_count = bnd.GetOverviewCount()
        if ovr_count > 0:
            idx = ovr_count + overview if overview < 0 else overview
            out = bnd.GetOverview(idx).ReadAsArray()
            bnd = ds = None
            return out
        logger.warning(f"Requested {overview = }, but none found for {filename}")

    # if rows or cols are not specified, load all rows/cols
    rows = slice(0, nrows) if rows in (None, slice(None)) else rows
    cols = slice(0, ncols) if cols in (None, slice(None)) else cols
    # Help out mypy:
    assert rows is not None
    assert cols is not None

    dt = gdal_to_numpy_type(ds.GetRasterBand(1).DataType)

    if isinstance(subsample_factor, int):
        subsample_factor = (subsample_factor, subsample_factor)

    xoff, yoff = int(cols.start), int(rows.start)
    row_stop = min(rows.stop, nrows)
    col_stop = min(cols.stop, ncols)
    xsize, ysize = int(col_stop - cols.start), int(row_stop - rows.start)
    if xsize <= 0 or ysize <= 0:
        msg = (
            f"Invalid row/col slices: {rows}, {cols} for file {filename} of size"
            f" {nrows}x{ncols}"
        )
        raise IndexError(msg)
    nrows_out, ncols_out = (
        ysize // subsample_factor[0],
        xsize // subsample_factor[1],
    )

    # Read the data, and decimate if specified
    resamp = gdal.GRA_NearestNeighbour
    if band is None:
        count = ds.RasterCount
        out = np.empty((count, nrows_out, ncols_out), dtype=dt)
        ds.ReadAsArray(xoff, yoff, xsize, ysize, buf_obj=out, resample_alg=resamp)
        if count == 1:
            out = out[0]
    else:
        out = np.empty((nrows_out, ncols_out), dtype=dt)
        bnd = ds.GetRasterBand(band)
        bnd.ReadAsArray(xoff, yoff, xsize, ysize, buf_obj=out, resample_alg=resamp)

    if not masked:
        return out
    # Get the nodata value
    nd = get_raster_nodata(filename)
    if nd is not None and np.isnan(nd):
        return np.ma.masked_invalid(out)
    else:
        return np.ma.masked_equal(out, nd)

process_blocks(readers, writer, func, block_shape=(512, 512), overlaps=(0, 0), num_threads=5)

Perform block-wise processing over blocks in readers, writing to writer.

Parameters:

Name	Type	Description	Default
readers	Sequence[StackReader]	Sequence of input readers to read data from.	required
writer	DatasetWriter	Output writer to write data to.	required
func	BlockProcessor	Function to process each block.	required
block_shape	tuple[int, int]	Shape of each block to process.	(512, 512)
overlaps	tuple[int, int]	Amount of overlap between blocks in (row, col) directions. By default (0, 0).	(0, 0)
num_threads	int	Number of threads to use, by default 5.	5

Source code in src/dolphin/io/_process.py

def process_blocks(
    readers: Sequence[StackReader],
    writer: DatasetWriter,
    func: BlockProcessor,
    block_shape: tuple[int, int] = (512, 512),
    overlaps: tuple[int, int] = (0, 0),
    num_threads: int = 5,
):
    """Perform block-wise processing over blocks in `readers`, writing to `writer`.

    Parameters
    ----------
    readers : Sequence[StackReader]
        Sequence of input readers to read data from.
    writer : DatasetWriter
        Output writer to write data to.
    func : BlockProcessor
        Function to process each block.
    block_shape : tuple[int, int], optional
        Shape of each block to process.
    overlaps : tuple[int, int], optional
        Amount of overlap between blocks in (row, col) directions.
        By default (0, 0).
    num_threads : int, optional
        Number of threads to use, by default 5.

    """
    block_manager = StridedBlockManager(
        arr_shape=readers[0].shape[-2:],
        block_shape=block_shape,
        # Here we are not using the striding mechanism
        strides=Strides(1, 1),
        # The "half window" is how much overlap we read in
        half_window=HalfWindow(*overlaps),
    )
    total_blocks = sum(1 for _ in block_manager.iter_blocks())
    pbar = tqdm(total=total_blocks)

    def write_callback(fut: Future):
        data, rows, cols = fut.result()
        writer[..., rows, cols] = data
        pbar.update()

    Executor = ThreadPoolExecutor if num_threads > 1 else DummyProcessPoolExecutor
    futures: set[Future] = set()

    # Create a helper function to perform the trimming after processing
    def _run_and_trim(
        func,
        readers,
        out_idxs: BlockIndices,
        trim_idxs: BlockIndices,
        in_idxs: BlockIndices,
    ):
        in_rows, in_cols = in_idxs
        out_rows, out_cols = out_idxs
        trim_rows, trim_cols = trim_idxs
        out_data, _, _ = func(readers=readers, rows=in_rows, cols=in_cols)
        return out_data[..., trim_rows, trim_cols], out_rows, out_cols

    with Executor(num_threads) as exc:
        for out_idxs, trim_idxs, in_idxs, _, _ in block_manager.iter_blocks():
            future = exc.submit(
                _run_and_trim, func, readers, out_idxs, trim_idxs, in_idxs
            )
            future.add_done_callback(write_callback)
            futures.add(future)

        while futures:
            done, futures = wait(futures, timeout=1, return_when=FIRST_EXCEPTION)
            for future in done:
                e = future.exception()
                if e is not None:
                    raise e

repack_raster(raster_path, output_dir=None, keep_bits=None, block_shape=(1024, 1024), **output_options)

Repack a single raster file with GDAL Translate using provided options.

Parameters:

Name	Type	Description	Default
raster_path	Path	Path to the input raster file.	required
output_dir	Path	Directory to save the repacked rasters or None for in-place repacking.	None
keep_bits	int	Number of bits to preserve in mantissa. Defaults to None. Lower numbers will truncate the mantissa more and enable more compression.	None
block_shape	int | tuple[int, int]	Size of blocks to read in at one time.	(1024, 1024)
**output_options		Keyword args passed to get_gtiff_options	{}

Returns:

Name	Type	Description
output_path	Path	Path to newly created file. If output_dir is None, this will be the same filename as raster_path

Source code in src/dolphin/io/_utils.py

def repack_raster(
    raster_path: Path,
    output_dir: Path | None = None,
    keep_bits: int | None = None,
    block_shape: int | tuple[int, int] = (1024, 1024),
    **output_options,
) -> Path:
    """Repack a single raster file with GDAL Translate using provided options.

    Parameters
    ----------
    raster_path : Path
        Path to the input raster file.
    output_dir : Path, optional
        Directory to save the repacked rasters or None for in-place repacking.
    keep_bits : int, optional
        Number of bits to preserve in mantissa. Defaults to None.
        Lower numbers will truncate the mantissa more and enable more compression.
    block_shape: int | tuple[int, int]
        Size of blocks to read in at one time.
    **output_options
        Keyword args passed to `get_gtiff_options`

    Returns
    -------
    output_path : Path
        Path to newly created file.
        If `output_dir` is None, this will be the same filename as `raster_path`

    """
    import rasterio as rio
    from rasterio.windows import Window

    from ._blocks import iter_blocks

    if isinstance(block_shape, int):
        block_shape = (block_shape, block_shape)

    if output_dir is None:
        output_file = tempfile.NamedTemporaryFile(  # noqa: SIM115
            suffix=raster_path.suffix, dir=output_dir, delete=False
        )
        output_path = Path(output_file.name)
    else:
        output_dir.mkdir(parents=True, exist_ok=True)
        output_path = output_dir / raster_path.name

    options = get_gtiff_options(**output_options)

    with rio.open(raster_path) as src:
        profile = src.profile
        profile.update(**options)
        # Work in blocks on the input raster
        blocks = iter_blocks(
            arr_shape=(src.height, src.width),
            block_shape=block_shape,
        )

        with rio.open(output_path, "w", **profile) as dst:
            for i in range(1, src.count + 1):
                for row_slice, col_slice in blocks:
                    window = Window.from_slices(rows=row_slice, cols=col_slice)
                    data = src.read(i, window=window)

                    if keep_bits is not None:
                        round_mantissa(data, keep_bits)

                    dst.write(data, i, window=window)

    if output_dir is None:
        # Overwrite the original
        shutil.move(output_path, raster_path)
        output_path = raster_path

    return output_path

repack_rasters(raster_files, output_dir=None, num_threads=4, keep_bits=None, block_shape=(1024, 1024), **output_options)

Recreate and compress a list of raster files.

Useful for rasters which were created in block and lost the full effect of compression.

Parameters:

Name	Type	Description	Default
raster_files	List[Path]	List of paths to the input raster files.	required
output_dir	Path	Directory to save the processed rasters or None for in-place processing.	None
num_threads	int	Number of threads to use (default is 4).	4
keep_bits	int	Number of bits to preserve in mantissa. Defaults to None. Lower numbers will truncate the mantissa more and enable more compression.	None
block_shape	int | tuple[int, int]	Size of blocks to read in at one time.	(1024, 1024)
**output_options		Creation options to pass to get_gtiff_options	{}

Returns:

Name	Type	Description
output_path	Path	Path to newly created file. If output_dir is None, this will be the same as raster_paths

Source code in src/dolphin/io/_utils.py

def repack_rasters(
    raster_files: list[Path],
    output_dir: Path | None = None,
    num_threads: int = 4,
    keep_bits: int | None = None,
    block_shape: int | tuple[int, int] = (1024, 1024),
    **output_options,
):
    """Recreate and compress a list of raster files.

    Useful for rasters which were created in block and lost
    the full effect of compression.

    Parameters
    ----------
    raster_files : List[Path]
        List of paths to the input raster files.
    output_dir : Path, optional
        Directory to save the processed rasters or None for in-place processing.
    num_threads : int, optional
        Number of threads to use (default is 4).
    keep_bits : int, optional
        Number of bits to preserve in mantissa. Defaults to None.
        Lower numbers will truncate the mantissa more and enable more compression.
    block_shape: int | tuple[int, int]
        Size of blocks to read in at one time.
    **output_options
        Creation options to pass to `get_gtiff_options`

    Returns
    -------
    output_path : Path
        Path to newly created file.
        If `output_dir` is None, this will be the same as `raster_paths`

    """
    from tqdm.contrib.concurrent import thread_map

    thread_map(
        lambda raster: repack_raster(
            raster,
            output_dir,
            keep_bits=keep_bits,
            block_shape=block_shape,
            **output_options,
        ),
        raster_files,
        max_workers=num_threads,
        desc="Processing Rasters",
    )

round_mantissa(z, keep_bits=10, chunk_rows=1024)

Zero out mantissa bits of elements of array in place.

Drops a specified number of bits from the floating point mantissa, leaving an array more amenable to compression.

Parameters:

Name	Type	Description	Default
z	ndarray	Real or complex array whose mantissas are to be zeroed out	required
keep_bits	int	Number of bits to preserve in mantissa. Lower numbers will truncate the mantissa more and enable more compression. Default is 10.	10
chunk_rows	int	Number of rows to process at a time to limit memory usage. Default is 1024.	1024

References

https://numcodecs.readthedocs.io/en/v0.12.1/_modules/numcodecs/bitround.html

Source code in src/dolphin/io/_utils.py

def round_mantissa(z: np.ndarray, keep_bits: int = 10, chunk_rows: int = 1024) -> None:
    """Zero out mantissa bits of elements of array in place.

    Drops a specified number of bits from the floating point mantissa,
    leaving an array more amenable to compression.

    Parameters
    ----------
    z : numpy.ndarray
        Real or complex array whose mantissas are to be zeroed out
    keep_bits : int
        Number of bits to preserve in mantissa.
        Lower numbers will truncate the mantissa more and enable
        more compression.
        Default is 10.
    chunk_rows : int
        Number of rows to process at a time to limit memory usage.
        Default is 1024.

    References
    ----------
    https://numcodecs.readthedocs.io/en/v0.12.1/_modules/numcodecs/bitround.html

    """
    max_bits = {
        "float16": 10,
        "float32": 23,
        "float64": 52,
    }
    # recurse for complex data
    if np.iscomplexobj(z):
        round_mantissa(z.real, keep_bits, chunk_rows)
        round_mantissa(z.imag, keep_bits, chunk_rows)
        return

    if not z.dtype.kind == "f" or z.dtype.itemsize > 8:
        raise TypeError("Only float arrays (16-64bit) can be bit-rounded")

    bits = max_bits[str(z.dtype)]
    # cast float to int type of same width (preserve endianness)
    a_int_dtype = np.dtype(z.dtype.str.replace("f", "i"))
    all_set = np.array(-1, dtype=a_int_dtype)
    if keep_bits == bits:
        return z
    if keep_bits > bits:
        raise ValueError("keep_bits too large for given dtype")
    maskbits = bits - keep_bits
    mask = (all_set >> maskbits) << maskbits
    half_quantum1 = (1 << (maskbits - 1)) - 1
    for i in range(0, z.shape[0], chunk_rows):
        b = z[i : i + chunk_rows].view(a_int_dtype)
        b += ((b >> maskbits) & 1) + half_quantum1
        b &= mask

set_raster_description(filename, description, band=1)

Set description on a raster band.

Parameters:

Name	Type	Description	Default
filename	Filename	Path to the file to load.	required
description	str	Description to set.	required
band	int	Band to set description for. Default is 1.	1

Source code in src/dolphin/io/_core.py

def set_raster_description(filename: Filename, description: str, band: int = 1):
    """Set description on a raster band.

    Parameters
    ----------
    filename : Filename
        Path to the file to load.
    description : str
        Description to set.
    band : int, optional
        Band to set description for. Default is 1.

    """
    ds = gdal.Open(fspath(filename), gdal.GA_Update)
    bnd = ds.GetRasterBand(band)
    bnd.SetDescription(description)
    bnd.FlushCache()
    ds = None

set_raster_metadata(filename, metadata, domain='')

Set metadata on a raster file.

Parameters:

Name	Type	Description	Default
filename	Filename	Path to the file to load.	required
metadata	dict	Dictionary of metadata to set.	required
domain	str	Domain to set metadata for. Default is "" (all domains).	''

Source code in src/dolphin/io/_core.py

def set_raster_metadata(
    filename: Filename, metadata: Mapping[str, Any], domain: str = ""
):
    """Set metadata on a raster file.

    Parameters
    ----------
    filename : Filename
        Path to the file to load.
    metadata : dict
        Dictionary of metadata to set.
    domain : str, optional
        Domain to set metadata for. Default is "" (all domains).

    """
    ds = gdal.Open(fspath(filename), gdal.GA_Update)
    # Ensure the keys/values are written as strings
    md_dict = {k: str(v) for k, v in metadata.items()}
    ds.SetMetadata(md_dict, domain)
    ds.FlushCache()
    ds = None

set_raster_nodata(filename, nodata, band=None)

Set the nodata value for a raster.

Parameters:

Name	Type	Description	Default
filename	Filename	Path to the file to load.	required
nodata	float	The nodata value to set.	required
band	int	The band to set the nodata value for, by default None (sets the nodata value for all bands).	None

Source code in src/dolphin/io/_core.py

def set_raster_nodata(filename: Filename, nodata: float, band: int | None = None):
    """Set the nodata value for a raster.

    Parameters
    ----------
    filename : Filename
        Path to the file to load.
    nodata : float
        The nodata value to set.
    band : int, optional
        The band to set the nodata value for, by default None
        (sets the nodata value for all bands).

    """
    ds = gdal.Open(fspath(filename), gdal.GA_Update)
    if band is None:
        for i in range(ds.RasterCount):
            ds.GetRasterBand(i + 1).SetNoDataValue(nodata)
    else:
        ds.GetRasterBand(band).SetNoDataValue(nodata)
    ds = None

set_raster_units(filename, units, band=None)

Set units on a raster band.

Parameters:

Name	Type	Description	Default
filename	Filename	Path to the file to load.	required
units	str	Units to set.	required
band	int	Band to set units for. Default is None, which sets for all bands.	None

Source code in src/dolphin/io/_core.py

def set_raster_units(filename: Filename, units: str, band: int | None = None) -> None:
    """Set units on a raster band.

    Parameters
    ----------
    filename : Filename
        Path to the file to load.
    units : str
        Units to set.
    band : int, optional
        Band to set units for. Default is None, which sets for all bands.

    """
    ds = gdal.Open(fspath(filename), gdal.GA_Update)
    bands = range(1, ds.RasterCount + 1) if band is None else range(band, band + 1)
    for i in bands:
        bnd = ds.GetRasterBand(i)
        bnd.SetUnitType(units)
        bnd.FlushCache()
    ds = None

write_arr(*, arr, output_name, like_filename=None, driver='GTiff', options=None, nbands=None, shape=None, dtype=None, geotransform=None, strides=None, projection=None, nodata=None, units=None, description=None)

Save an array to output_name.

If like_filename if provided, copies the projection/nodata. Options can be overridden by passing driver/nbands/dtype.

If arr is None, create an empty file with the same x/y shape as like_filename.

Parameters:

Name	Type	Description	Default
arr	ArrayLike	Array to save. If None, create an empty file.	required
output_name	str or Path	Path to save the file to.	required
like_filename	str or Path	Path to a file to copy raster shape/metadata from.	None
driver	str	GDAL driver to use. Default is "GTiff".	'GTiff'
options	list	list of options to pass to the driver. Default is DEFAULT_TIFF_OPTIONS.	None
nbands	int	Number of bands to save. Default is 1.	None
shape	tuple	(rows, cols) of desired output file. Overrides the shape of the output file, if using like_filename.	None
dtype	DTypeLike	Data type to save. Default is arr.dtype or the datatype of like_filename.	None
geotransform	list	Geotransform to save. Default is the geotransform of like_filename. See https://gdal.org/tutorials/geotransforms_tut.html .	None
strides	dict	If using like_filename, used to change the pixel size of the output file. {"x": x strides, "y": y strides}	None
projection	str or int	Projection to save. Default is the projection of like_filename. Possible values are anything parse-able by pyproj.CRS.from_user_input (including EPSG ints, WKT strings, PROJ strings, etc.)	None
nodata	float	Nodata value to save. Default is the nodata of band 1 of like_filename (if provided), or None.	None
units	str	Units of the data. Default is None. Value is stored in the metadata as "units".	None
description	str	Description of the raster bands stored in the metadata.	None

Source code in src/dolphin/io/_core.py

def write_arr(
    *,
    arr: Optional[ArrayLike],
    output_name: Filename,
    like_filename: Optional[Filename] = None,
    driver: Optional[str] = "GTiff",
    options: Optional[Sequence] = None,
    nbands: Optional[int] = None,
    shape: Optional[tuple[int, int]] = None,
    dtype: Optional[DTypeLike] = None,
    geotransform: Optional[Sequence[float]] = None,
    strides: Optional[dict[str, int]] = None,
    projection: Optional[Any] = None,
    nodata: Optional[float] = None,
    units: Optional[str] = None,
    description: Optional[str] = None,
):
    """Save an array to `output_name`.

    If `like_filename` if provided, copies the projection/nodata.
    Options can be overridden by passing `driver`/`nbands`/`dtype`.

    If arr is None, create an empty file with the same x/y shape as `like_filename`.

    Parameters
    ----------
    arr : ArrayLike, optional
        Array to save. If None, create an empty file.
    output_name : str or Path
        Path to save the file to.
    like_filename : str or Path, optional
        Path to a file to copy raster shape/metadata from.
    driver : str, optional
        GDAL driver to use. Default is "GTiff".
    options : list, optional
        list of options to pass to the driver. Default is DEFAULT_TIFF_OPTIONS.
    nbands : int, optional
        Number of bands to save. Default is 1.
    shape : tuple, optional
        (rows, cols) of desired output file.
        Overrides the shape of the output file, if using `like_filename`.
    dtype : DTypeLike, optional
        Data type to save. Default is `arr.dtype` or the datatype of like_filename.
    geotransform : list, optional
        Geotransform to save. Default is the geotransform of like_filename.
        See https://gdal.org/tutorials/geotransforms_tut.html .
    strides : dict, optional
        If using `like_filename`, used to change the pixel size of the output file.
        {"x": x strides, "y": y strides}
    projection : str or int, optional
        Projection to save. Default is the projection of like_filename.
        Possible values are anything parse-able by ``pyproj.CRS.from_user_input``
        (including EPSG ints, WKT strings, PROJ strings, etc.)
    nodata : float, optional
        Nodata value to save.
        Default is the nodata of band 1 of `like_filename` (if provided), or None.
    units : str, optional
        Units of the data. Default is None.
        Value is stored in the metadata as "units".
    description : str, optional
        Description of the raster bands stored in the metadata.

    """
    fi = FileInfo.from_user_inputs(
        arr=arr,
        output_name=output_name,
        like_filename=like_filename,
        driver=driver,
        options=options,
        nbands=nbands,
        shape=shape,
        dtype=dtype,
        geotransform=geotransform,
        strides=strides,
        projection=projection,
        nodata=nodata,
    )
    drv = gdal.GetDriverByName(fi.driver)
    ds_out = drv.Create(
        fspath(output_name),
        fi.xsize,
        fi.ysize,
        fi.nbands,
        fi.gdal_dtype,
        options=fi.options,
    )

    # Set the geo/proj information
    if fi.projection:
        # Make sure we're got a correct format for the projection
        # this still works if we're passed a WKT string
        proj = CRS.from_user_input(fi.projection).to_wkt()
        ds_out.SetProjection(proj)

    if fi.geotransform is not None:
        ds_out.SetGeoTransform(fi.geotransform)

    # Set the nodata/units/description for each band
    for i in range(fi.nbands):
        logger.debug(f"Setting nodata for band {i + 1}/{fi.nbands}")
        bnd = ds_out.GetRasterBand(i + 1)
        # Note: right now we're assuming the nodata/units/description
        if fi.nodata is not None:
            bnd.SetNoDataValue(fi.nodata)
        if units is not None:
            bnd.SetUnitType(units)
        if description is not None:
            bnd.SetDescription(description)

    # Write the actual data
    if arr is not None:
        if arr.ndim == 2:
            arr = arr[np.newaxis, ...]
        for i in range(fi.nbands):
            logger.debug(f"Writing band {i + 1}/{fi.nbands}")
            bnd = ds_out.GetRasterBand(i + 1)
            bnd.WriteArray(arr[i])

    ds_out.FlushCache()
    ds_out = None

write_block(cur_block, filename, row_start, col_start, band=None, dset=None)

Write out an ndarray to a subset of the pre-made filename.

Parameters:

Name	Type	Description	Default
cur_block	ArrayLike	2D or 3D data array	required
filename	Filename	list of output files to save to, or (if cur_block is 2D) a single file.	required
row_start	int	Row index to start writing at.	required
col_start	int	Column index to start writing at.	required
band	int	Raster band to write to within filename. If None, writes to band 1 (for 2D), or all bands if cur_block.ndim = 3.	None
dset	str	(For writing to HDF5/NetCDF files) The name of the string dataset within filename to write to.	None

Raises:

Type	Description
ValueError	If length of output_files does not match length of cur_block.

Source code in src/dolphin/io/_core.py

def write_block(
    cur_block: NDArray,
    filename: Filename,
    row_start: int,
    col_start: int,
    band: int | None = None,
    dset: str | None = None,
):
    """Write out an ndarray to a subset of the pre-made `filename`.

    Parameters
    ----------
    cur_block : ArrayLike
        2D or 3D data array
    filename : Filename
        list of output files to save to, or (if cur_block is 2D) a single file.
    row_start : int
        Row index to start writing at.
    col_start : int
        Column index to start writing at.
    band : int, optional
        Raster band to write to within `filename`.
        If None, writes to band 1 (for 2D), or all bands if `cur_block.ndim = 3`.
    dset : str
        (For writing to HDF5/NetCDF files) The name of the string dataset
        within `filename` to write to.

    Raises
    ------
    ValueError
        If length of `output_files` does not match length of `cur_block`.

    """
    if cur_block.ndim == 2 and band is None:
        # Make into 3D array shaped (1, rows, cols)
        cur_block = cur_block[np.newaxis, ...]
    # filename must be pre-made
    filename = Path(filename)
    if not filename.exists():
        msg = f"File {filename} does not exist"
        raise ValueError(msg)

    if filename.suffix in (".h5", ".hdf5", ".nc"):
        if dset is None:
            raise ValueError("Missing `dset` argument for writing to HDF5")
        _write_hdf5(cur_block, filename, row_start, col_start, dset)
    else:
        _write_gdal(cur_block, filename, row_start, col_start, band)