<commit_before>"""Rasterio input/output."""

from __future__ import absolute_import

include "directives.pxi"
include "gdal.pxi"

import logging
import os
import os.path
import sys
import uuid
import warnings

import numpy as np

from rasterio._base import tastes_like_gdal
from rasterio._env import driver_count, GDALEnv
from rasterio._err import (
    GDALError, CPLE_OpenFailedError, CPLE_IllegalArgError)
from rasterio.crs import CRS
from rasterio.compat import text_type, string_types
from rasterio import dtypes
from rasterio.enums import ColorInterp, MaskFlags, Resampling
from rasterio.errors import CRSError, DriverRegistrationError
from rasterio.errors import RasterioIOError
from rasterio.errors import NodataShadowWarning
from rasterio.sample import sample_gen
from rasterio.transform import Affine
from rasterio.vfs import parse_path, vsi_path
from rasterio import windows

from libc.stdio cimport FILE
cimport numpy as np

from rasterio._base cimport _osr_from_crs, get_driver_name, DatasetBase
from rasterio._err cimport exc_wrap_int, exc_wrap_pointer, exc_wrap_vsilfile
from rasterio._shim cimport (
    delete_nodata_value, io_band, io_multi_band, io_multi_mask)


log = logging.getLogger(__name__)


cdef bint in_dtype_range(value, dtype):
    """Returns True if value is in the range of dtype, else False."""
    infos = {
        'c': np.finfo,
        'f': np.finfo,
        'i': np.iinfo,
        'u': np.iinfo,
        # Cython 0.22 returns dtype.kind as an int and will not cast to a char
        99: np.finfo,
        102: np.finfo,
        105: np.iinfo,
        117: np.iinfo}

    key = np.dtype(dtype).kind
    if np.isnan(value):
        return key in ('c', 'f', 99, 102)

    rng = infos[key](dtype)
    return rng.min <= value <= rng.max


cdef int io_auto(data, GDALRasterBandH band, bint write, int resampling=0):
    """Convenience function to handle IO with a GDAL band.

    :param data: a numpy ndarray
    :param band: an instance of GDALGetRasterBand
    :param write: 1 (True) uses write mode (writes data into band),
                  0 (False) uses read mode (reads band into data)
    :return: the return value from the data-type specific IO function
    """
    cdef int ndims = len(data.shape)
    cdef float height = data.shape[-2]
    cdef float width = data.shape[-1]

    if ndims == 2:
        return io_band(band, write, 0.0, 0.0, width, height, data,
                       resampling=resampling)
    elif ndims == 3:
        indexes = np.arange(1, data.shape[0] + 1)
        return io_multi_band(band, write, 0.0, 0.0, width, height, data,
                             indexes, resampling=resampling)
    else:
        raise ValueError("Specified data must have 2 or 3 dimensions")


cdef class DatasetReaderBase(DatasetBase):

    def read(self, indexes=None, out=None, window=None, masked=False,
            out_shape=None, boundless=False, resampling=Resampling.nearest,
            fill_value=None):
        """Read raster bands as a multidimensional array

        Parameters
        ----------
        indexes : list of ints or a single int, optional
            If `indexes` is a list, the result is a 3D array, but is
            a 2D array if it is a band index number.

        out : numpy ndarray, optional
            As with Numpy ufuncs, this is an optional reference to an
            output array into which data will be placed. If the height
            and width of `out` differ from that of the specified 
            window (see below), the raster image will be decimated or
            replicated using the specified resampling method (also see
            below).

            *Note*: the method's return value may be a view on this
            array. In other words, `out` is likely to be an
            incomplete representation of the method's results.

            This parameter cannot be combined with `out_shape`.

        out_shape : tuple, optional
            A tuple describing the shape of a new output array. See
            `out` (above) for notes on image decimation and
            replication.

            Cannot combined with `out`.

        window : a pair (tuple) of pairs of ints, optional
            The optional `window` argument is a 2 item tuple. The first
            item is a tuple containing the indexes of the rows at which
            the window starts and stops and the second is a tuple
            containing the indexes of the columns at which the window
            starts and stops. For example, ((0, 2), (0, 2)) defines
            a 2x2 window at the upper left of the raster dataset.

        masked : bool, optional
            If `masked` is `True` the return value will be a masked
            array. Otherwise (the default) the return value will be a
            regular array. Masks will be exactly the inverse of the
            GDAL RFC 15 conforming arrays returned by read_masks().

        boundless : bool, optional (default `False`)
            If `True`, windows that extend beyond the dataset's extent
            are permitted and partially or completely filled arrays will
            be returned as appropriate.

        fill_value : scalar
            Fill value applied in the `boundless=True` case only.

        Returns
        -------
        Numpy ndarray or a view on a Numpy ndarray

        Note: as with Numpy ufuncs, an object is returned even if you
        use the optional `out` argument and the return value shall be
        preferentially used by callers.
        """

        cdef GDALRasterBandH band = NULL

        return2d = False
        if indexes is None:
            indexes = self.indexes
        elif isinstance(indexes, int):
            indexes = [indexes]
            return2d = True
            if out is not None and out.ndim == 2:
                out.shape = (1,) + out.shape
        if not indexes:
            raise ValueError("No indexes to read")

        check_dtypes = set()
        nodatavals = []
        # Check each index before processing 3D array
        for bidx in indexes:
            if bidx not in self.indexes:
                raise IndexError("band index out of range")
            idx = self.indexes.index(bidx)

            dtype = self.dtypes[idx]
            check_dtypes.add(dtype)

            ndv = self._nodatavals[idx]

            log.debug("Output nodata value read from file: %r", ndv)

            # Change given nodatavals to the closest value that
            # can be represented by this band's data type to
            # match GDAL's strategy.
            if fill_value:
                ndv = fill_value
                log.debug("Output nodata value set from fill value")

            elif ndv is not None:
                if np.dtype(dtype).kind in ('i', 'u'):
                    info = np.iinfo(dtype)
                    dt_min, dt_max = info.min, info.max
                elif np.dtype(dtype).kind in ('f', 'c'):
                    info = np.finfo(dtype)
                    dt_min, dt_max = info.min, info.max
                else:
                    dt_min, dt_max = False, True
                if ndv < dt_min:
                    ndv = dt_min
                elif ndv > dt_max:
                    ndv = dt_max

            nodatavals.append(ndv)

        log.debug("Output nodata values: %r", nodatavals)

        # Mixed dtype reads are not supported at this time.
        if len(check_dtypes) > 1:
            raise ValueError("more than one 'dtype' found")
        elif len(check_dtypes) == 0:
            dtype = self.dtypes[0]
        else:
            dtype = check_dtypes.pop()

        # Get the natural shape of the read window, boundless or not.
        win_shape = (len(indexes),)
        if window:
            if boundless:
                win_shape += (
                        window[0][1]-window[0][0], window[1][1]-window[1][0])
            else:
                window = windows.crop(
                    windows.evaluate(window, self.height, self.width),
                    self.height, self.width
                )
                (r_start, r_stop), (c_start, c_stop) = window
                win_shape += (r_stop - r_start, c_stop - c_start)
        else:
            win_shape += self.shape

        if out is not None and out_shape is not None:
            raise ValueError("out and out_shape are exclusive")

        # `out` takes precedence over `out_shape`.
        elif out is not None:
            if out.dtype != dtype:
                raise ValueError(
                    "the array's dtype '%s' does not match "
                    "the file's dtype '%s'" % (out.dtype, dtype))
            if out.shape[0] != win_shape[0]:
                raise ValueError(
                    "'out' shape %s does not match window shape %s" %
                    (out.shape, win_shape))

        else:
            if out_shape is not None:
                if len(out_shape) == 2:
                    out_shape = (1,) + out_shape
            else:
                out_shape = win_shape

            # We're filling in both the bounded and boundless cases.
            # TODO: profile and see if we should avoid this in the
            # bounded case.

            if boundless:
                out = np.zeros(tuple(int(round(x)) for x in out_shape), dtype=dtype)
            else:
                out = np.empty(tuple(int(round(x)) for x in out_shape), dtype=dtype)

            for i, (ndv, arr) in enumerate(zip(
                    nodatavals, out if len(out.shape) == 3 else [out])):

                if ndv is not None:
                    arr.fill(ndv)

        # Masking
        # -------
        #
        # If masked is True, we check the GDAL mask flags using
        # GDALGetMaskFlags. If GMF_ALL_VALID for all bands, we do not
        # call read_masks(), but pass `mask=False` to the masked array
        # constructor. Else, we read the GDAL mask bands using
        # read_masks(), invert them and use them in constructing masked
        # arrays.

        if masked:
            enums = self.mask_flag_enums
            all_valid = all([MaskFlags.all_valid in flags for flags in enums])
            log.debug("all_valid: %s", all_valid)
            log.debug("mask_flags: %r", enums)

        # We can jump straight to _read() in some cases. We can ignore
        # the boundless flag if there's no given window.
        if not boundless or not window:

            log.debug("Jump straight to _read()")
            log.debug("Window: %r", window)

            out = self._read(indexes, out, window, dtype,
                             resampling=resampling)

            if masked:
                if all_valid:
                    mask = np.ma.nomask
                else:
                    mask = np.zeros(out.shape, 'uint8')
                    mask = ~self._read(
                        indexes, mask, window, 'uint8', masks=True,
                        resampling=resampling).astype('bool')

                kwds = {'mask': mask}
                # Set a fill value only if the read bands share a
                # single nodata value.
                if len(set(nodatavals)) == 1:
                    if nodatavals[0] is not None:
                        kwds['fill_value'] = nodatavals[0]
                out = np.ma.array(out, **kwds)

        else:
            # Compute the overlap between the dataset and the boundless window.
            overlap = ((
                max(min(window[0][0], self.height), 0),
                max(min(window[0][1], self.height), 0)), (
                max(min(window[1][0], self.width), 0),
                max(min(window[1][1], self.width), 0)))

            log.debug("Overlap: %r", overlap)

            if overlap != ((0, 0), (0, 0)):
                # Prepare a buffer.
                window_h, window_w = win_shape[-2:]
                overlap_h = overlap[0][1] - overlap[0][0]
                overlap_w = overlap[1][1] - overlap[1][0]
                scaling_h = float(out.shape[-2:][0])/window_h
                scaling_w = float(out.shape[-2:][1])/window_w
                buffer_shape = (
                        int(round(overlap_h*scaling_h)),
                        int(round(overlap_w*scaling_w)))
                data = np.zeros(win_shape[:-2] + buffer_shape, dtype)
                data = self._read(indexes, data, overlap, dtype,
                                  resampling=resampling)

                if masked:
                    mask = np.zeros(win_shape[:-2] + buffer_shape, 'uint8')
                    mask = ~self._read(
                        indexes, mask, overlap, 'uint8', masks=True,
                        resampling=resampling).astype('bool')
                    kwds = {'mask': mask}
                    if len(set(nodatavals)) == 1:
                        if nodatavals[0] is not None:
                            kwds['fill_value'] = nodatavals[0]
                    data = np.ma.array(data, **kwds)

            else:
                data = None
                if masked:
                    kwds = {'mask': True}
                    if len(set(nodatavals)) == 1:
                        if nodatavals[0] is not None:
                            kwds['fill_value'] = nodatavals[0]
                    out = np.ma.array(out, **kwds)

            if data is not None:
                # Determine where to put the data in the output window.
                data_h, data_w = buffer_shape
                roff = 0
                coff = 0
                if window[0][0] < 0:
                    roff = int(-window[0][0] * scaling_h)
                if window[1][0] < 0:
                    coff = int(-window[1][0] * scaling_w)

                for dst, src in zip(
                        out if len(out.shape) == 3 else [out],
                        data if len(data.shape) == 3 else [data]):
                    dst[roff:roff+data_h, coff:coff+data_w] = src

                if masked:
                    if not hasattr(out, 'mask'):
                        kwds = {'mask': True}
                        if len(set(nodatavals)) == 1:
                            if nodatavals[0] is not None:
                                kwds['fill_value'] = nodatavals[0]
                        out = np.ma.array(out, **kwds)

                    for dst, src in zip(
                            out.mask if len(out.shape) == 3 else [out.mask],
                            data.mask if len(data.shape) == 3 else [data.mask]):
                        dst[roff:roff+data_h, coff:coff+data_w] = src

        if return2d:
            out.shape = out.shape[1:]

        return out


    def read_masks(self, indexes=None, out=None, out_shape=None, window=None,
                   boundless=False, resampling=Resampling.nearest):
        """Read raster band masks as a multidimensional array

        Parameters
        ----------
        indexes : list of ints or a single int, optional
            If `indexes` is a list, the result is a 3D array, but is
            a 2D array if it is a band index number.

        out : numpy ndarray, optional
            As with Numpy ufuncs, this is an optional reference to an
            output array with the same dimensions and shape into which
            data will be placed.

            *Note*: the method's return value may be a view on this
            array. In other words, `out` is likely to be an
            incomplete representation of the method's results.

            Cannot combine with `out_shape`.

        out_shape : tuple, optional
            A tuple describing the output array's shape.  Allows for decimated
            reads without constructing an output Numpy array.

            Cannot combined with `out`.

        window : a pair (tuple) of pairs of ints, optional
            The optional `window` argument is a 2 item tuple. The first
            item is a tuple containing the indexes of the rows at which
            the window starts and stops and the second is a tuple
            containing the indexes of the columns at which the window
            starts and stops. For example, ((0, 2), (0, 2)) defines
            a 2x2 window at the upper left of the raster dataset.

        boundless : bool, optional (default `False`)
            If `True`, windows that extend beyond the dataset's extent
            are permitted and partially or completely filled arrays will
            be returned as appropriate.

        Returns
        -------
        Numpy ndarray or a view on a Numpy ndarray

        Note: as with Numpy ufuncs, an object is returned even if you
        use the optional `out` argument and the return value shall be
        preferentially used by callers.
        """

        return2d = False
        if indexes is None:
            indexes = self.indexes
        elif isinstance(indexes, int):
            indexes = [indexes]
            return2d = True
            if out is not None and out.ndim == 2:
                out.shape = (1,) + out.shape
        if not indexes:
            raise ValueError("No indexes to read")

        # Get the natural shape of the read window, boundless or not.
        win_shape = (len(indexes),)
        if window:
            if boundless:
                win_shape += (
                        window[0][1]-window[0][0], window[1][1]-window[1][0])
            else:
                w = windows.evaluate(window, self.height, self.width)
                minr = min(max(w[0][0], 0), self.height)
                maxr = max(0, min(w[0][1], self.height))
                minc = min(max(w[1][0], 0), self.width)
                maxc = max(0, min(w[1][1], self.width))
                win_shape += (maxr - minr, maxc - minc)
                window = ((minr, maxr), (minc, maxc))
        else:
            win_shape += self.shape

        dtype = 'uint8'

        if out is not None and out_shape is not None:
            raise ValueError("out and out_shape are exclusive")
        elif out_shape is not None:
            if len(out_shape) == 2:
                out_shape = (1,) + out_shape
            out = np.zeros(out_shape, 'uint8')

        if out is not None:
            if out.dtype != np.dtype(dtype):
                raise ValueError(
                    "the out array's dtype '%s' does not match '%s'"
                    % (out.dtype, dtype))
            if out.shape[0] != win_shape[0]:
                raise ValueError(
                    "'out' shape %s does not match window shape %s" %
                    (out.shape, win_shape))
        else:
            out = np.zeros(win_shape, 'uint8')


        # We can jump straight to _read() in some cases. We can ignore
        # the boundless flag if there's no given window.
        if not boundless or not window:
            out = self._read(indexes, out, window, dtype, masks=True,
                             resampling=resampling)

        else:
            # Compute the overlap between the dataset and the boundless window.
            overlap = ((
                max(min(window[0][0], self.height), 0),
                max(min(window[0][1], self.height), 0)), (
                max(min(window[1][0], self.width), 0),
                max(min(window[1][1], self.width), 0)))

            if overlap != ((0, 0), (0, 0)):
                # Prepare a buffer.
                window_h, window_w = win_shape[-2:]
                overlap_h = overlap[0][1] - overlap[0][0]
                overlap_w = overlap[1][1] - overlap[1][0]
                scaling_h = float(out.shape[-2:][0])/window_h
                scaling_w = float(out.shape[-2:][1])/window_w
                buffer_shape = (int(overlap_h*scaling_h), int(overlap_w*scaling_w))
                data = np.zeros(win_shape[:-2] + buffer_shape, 'uint8')
                data = self._read(indexes, data, overlap, dtype, masks=True,
                                  resampling=resampling)
            else:
                data = None

            if data is not None:
                # Determine where to put the data in the output window.
                data_h, data_w = data.shape[-2:]
                roff = 0
                coff = 0
                if window[0][0] < 0:
                    roff = int(window_h*scaling_h) - data_h
                if window[1][0] < 0:
                    coff = int(window_w*scaling_w) - data_w
                for dst, src in zip(
                        out if len(out.shape) == 3 else [out],
                        data if len(data.shape) == 3 else [data]):
                    dst[roff:roff+data_h, coff:coff+data_w] = src

        if return2d:
            out.shape = out.shape[1:]

        return out


    def _read(self, indexes, out, window, dtype, masks=False,
              resampling=Resampling.nearest):
        """Read raster bands as a multidimensional array

        If `indexes` is a list, the result is a 3D array, but
        is a 2D array if it is a band index number.

        Optional `out` argument is a reference to an output array with the
        same dimensions and shape.

        See `read_band` for usage of the optional `window` argument.

        The return type will be either a regular NumPy array, or a masked
        NumPy array depending on the `masked` argument. The return type is
        forced if either `True` or `False`, but will be chosen if `None`.
        For `masked=None` (default), the array will be the same type as
        `out` (if used), or will be masked if any of the nodatavals are
        not `None`.
        """
        cdef int aix, bidx, indexes_count
        cdef double height, width, xoff, yoff
        cdef int retval = 0
        cdef GDALDatasetH dataset = NULL

        dataset = self.handle()

        # Prepare the IO window.
        if window:
            window = windows.evaluate(window, self.height, self.width)

            log.debug("Eval'd window: %r", window)

            # yoff = <int>window[0][0]
            # xoff = <int>window[1][0]
            # height = <int>window[0][1] - yoff
            # width = <int>window[1][1] - xoff
            yoff = window[0][0]
            xoff = window[1][0]
            height = window[0][1] - yoff
            width = window[1][1] - xoff
        else:
            xoff = yoff = <int>0
            width = <int>self.width
            height = <int>self.height

        log.debug(
            "IO window xoff=%s yoff=%s width=%s height=%s",
            xoff, yoff, width, height)

        # Call io_multi* functions with C type args so that they
        # can release the GIL.
        indexes_arr = np.array(indexes, dtype=int)
        indexes_count = <int>indexes_arr.shape[0]

        if masks:
            # Warn if nodata attribute is shadowing an alpha band.
            if self.count == 4 and self.colorinterp(4) == ColorInterp.alpha:
                for flags in self.mask_flag_enums:
                    if MaskFlags.nodata in flags:
                        warnings.warn(NodataShadowWarning())

            retval = io_multi_mask(
                            self._hds, 0, xoff, yoff, width, height,
                            out, indexes_arr, resampling=resampling)

        else:
            retval = io_multi_band(self._hds, 0, xoff, yoff, width, height,
                                  out, indexes_arr, resampling=resampling)
<commit_msg>Round windows to get array shape<commit_after>"""Rasterio input/output."""

from __future__ import absolute_import

include "directives.pxi"
include "gdal.pxi"

import logging
import os
import os.path
import sys
import uuid
import warnings

import numpy as np

from rasterio._base import tastes_like_gdal
from rasterio._env import driver_count, GDALEnv
from rasterio._err import (
    GDALError, CPLE_OpenFailedError, CPLE_IllegalArgError)
from rasterio.crs import CRS
from rasterio.compat import text_type, string_types
from rasterio import dtypes
from rasterio.enums import ColorInterp, MaskFlags, Resampling
from rasterio.errors import CRSError, DriverRegistrationError
from rasterio.errors import RasterioIOError
from rasterio.errors import NodataShadowWarning
from rasterio.sample import sample_gen
from rasterio.transform import Affine
from rasterio.vfs import parse_path, vsi_path
from rasterio import windows

from libc.stdio cimport FILE
from libc.math cimport ceil
cimport numpy as np

from rasterio._base cimport _osr_from_crs, get_driver_name, DatasetBase
from rasterio._err cimport exc_wrap_int, exc_wrap_pointer, exc_wrap_vsilfile
from rasterio._shim cimport (
    delete_nodata_value, io_band, io_multi_band, io_multi_mask)


log = logging.getLogger(__name__)


cdef bint in_dtype_range(value, dtype):
    """Returns True if value is in the range of dtype, else False."""
    infos = {
        'c': np.finfo,
        'f': np.finfo,
        'i': np.iinfo,
        'u': np.iinfo,
        # Cython 0.22 returns dtype.kind as an int and will not cast to a char
        99: np.finfo,
        102: np.finfo,
        105: np.iinfo,
        117: np.iinfo}

    key = np.dtype(dtype).kind
    if np.isnan(value):
        return key in ('c', 'f', 99, 102)

    rng = infos[key](dtype)
    return rng.min <= value <= rng.max


cdef int io_auto(data, GDALRasterBandH band, bint write, int resampling=0):
    """Convenience function to handle IO with a GDAL band.

    :param data: a numpy ndarray
    :param band: an instance of GDALGetRasterBand
    :param write: 1 (True) uses write mode (writes data into band),
                  0 (False) uses read mode (reads band into data)
    :return: the return value from the data-type specific IO function
    """
    cdef int ndims = len(data.shape)
    cdef float height = data.shape[-2]
    cdef float width = data.shape[-1]

    if ndims == 2:
        return io_band(band, write, 0.0, 0.0, width, height, data,
                       resampling=resampling)
    elif ndims == 3:
        indexes = np.arange(1, data.shape[0] + 1)
        return io_multi_band(band, write, 0.0, 0.0, width, height, data,
                             indexes, resampling=resampling)
    else:
        raise ValueError("Specified data must have 2 or 3 dimensions")


cdef class DatasetReaderBase(DatasetBase):

    def read(self, indexes=None, out=None, window=None, masked=False,
            out_shape=None, boundless=False, resampling=Resampling.nearest,
            fill_value=None):
        """Read raster bands as a multidimensional array

        Parameters
        ----------
        indexes : list of ints or a single int, optional
            If `indexes` is a list, the result is a 3D array, but is
            a 2D array if it is a band index number.

        out : numpy ndarray, optional
            As with Numpy ufuncs, this is an optional reference to an
            output array into which data will be placed. If the height
            and width of `out` differ from that of the specified 
            window (see below), the raster image will be decimated or
            replicated using the specified resampling method (also see
            below).

            *Note*: the method's return value may be a view on this
            array. In other words, `out` is likely to be an
            incomplete representation of the method's results.

            This parameter cannot be combined with `out_shape`.

        out_shape : tuple, optional
            A tuple describing the shape of a new output array. See
            `out` (above) for notes on image decimation and
            replication.

            Cannot combined with `out`.

        window : a pair (tuple) of pairs of ints, optional
            The optional `window` argument is a 2 item tuple. The first
            item is a tuple containing the indexes of the rows at which
            the window starts and stops and the second is a tuple
            containing the indexes of the columns at which the window
            starts and stops. For example, ((0, 2), (0, 2)) defines
            a 2x2 window at the upper left of the raster dataset.

        masked : bool, optional
            If `masked` is `True` the return value will be a masked
            array. Otherwise (the default) the return value will be a
            regular array. Masks will be exactly the inverse of the
            GDAL RFC 15 conforming arrays returned by read_masks().

        boundless : bool, optional (default `False`)
            If `True`, windows that extend beyond the dataset's extent
            are permitted and partially or completely filled arrays will
            be returned as appropriate.

        fill_value : scalar
            Fill value applied in the `boundless=True` case only.

        Returns
        -------
        Numpy ndarray or a view on a Numpy ndarray

        Note: as with Numpy ufuncs, an object is returned even if you
        use the optional `out` argument and the return value shall be
        preferentially used by callers.
        """

        cdef GDALRasterBandH band = NULL

        return2d = False
        if indexes is None:
            indexes = self.indexes
        elif isinstance(indexes, int):
            indexes = [indexes]
            return2d = True

            if out is not None and out.ndim == 2:
                out.shape = (1,) + out.shape

        if not indexes:
            raise ValueError("No indexes to read")

        check_dtypes = set()
        nodatavals = []
        # Check each index before processing 3D array
        for bidx in indexes:
            if bidx not in self.indexes:
                raise IndexError("band index out of range")
            idx = self.indexes.index(bidx)

            dtype = self.dtypes[idx]
            check_dtypes.add(dtype)

            ndv = self._nodatavals[idx]

            log.debug("Output nodata value read from file: %r", ndv)

            # Change given nodatavals to the closest value that
            # can be represented by this band's data type to
            # match GDAL's strategy.
            if fill_value:
                ndv = fill_value
                log.debug("Output nodata value set from fill value")

            elif ndv is not None:
                if np.dtype(dtype).kind in ('i', 'u'):
                    info = np.iinfo(dtype)
                    dt_min, dt_max = info.min, info.max
                elif np.dtype(dtype).kind in ('f', 'c'):
                    info = np.finfo(dtype)
                    dt_min, dt_max = info.min, info.max
                else:
                    dt_min, dt_max = False, True
                if ndv < dt_min:
                    ndv = dt_min
                elif ndv > dt_max:
                    ndv = dt_max

            nodatavals.append(ndv)

        log.debug("Output nodata values: %r", nodatavals)

        # Mixed dtype reads are not supported at this time.
        if len(check_dtypes) > 1:
            raise ValueError("more than one 'dtype' found")
        elif len(check_dtypes) == 0:
            dtype = self.dtypes[0]
        else:
            dtype = check_dtypes.pop()

        # Get the natural shape of the read window, boundless or not.
        # The window can have float values. In this case, we round up
        # when computing the shape.

        # Stub the win_shape.
        win_shape = (len(indexes),)

        if window:

            if boundless:
                win_shape += (
                    int(round(window[0][1] - window[0][0], 6)),
                    int(round(window[1][1] - window[1][0], 6)))

            else:
                window = windows.crop(
                    windows.evaluate(window, self.height, self.width),
                    self.height, self.width)
                (r_start, r_stop), (c_start, c_stop) = window

                log.debug("Cropped window: %r", window)

                win_shape += (
                    int(round(r_stop - r_start, 6)),
                    int(round(c_stop - c_start, 6)))

        else:
            win_shape += self.shape

        if out is not None and out_shape is not None:
            raise ValueError("out and out_shape are exclusive")

        # `out` takes precedence over `out_shape`.
        elif out is not None:
            if out.dtype != dtype:
                raise ValueError(
                    "the array's dtype '%s' does not match "
                    "the file's dtype '%s'" % (out.dtype, dtype))
            if out.shape[0] != win_shape[0]:
                raise ValueError(
                    "'out' shape %s does not match window shape %s" %
                    (out.shape, win_shape))

        else:
            if out_shape is not None:
                if len(out_shape) == 2:
                    out_shape = (1,) + out_shape
            else:
                out_shape = win_shape

            # We're filling in both the bounded and boundless cases.
            # TODO: profile and see if we should avoid this in the
            # bounded case.

            if boundless:
                out = np.zeros(tuple(int(round(x)) for x in out_shape), dtype=dtype)
            else:
                out = np.empty(tuple(int(round(x)) for x in out_shape), dtype=dtype)

            for i, (ndv, arr) in enumerate(zip(
                    nodatavals, out if len(out.shape) == 3 else [out])):

                if ndv is not None:
                    arr.fill(ndv)

        # Masking
        # -------
        #
        # If masked is True, we check the GDAL mask flags using
        # GDALGetMaskFlags. If GMF_ALL_VALID for all bands, we do not
        # call read_masks(), but pass `mask=False` to the masked array
        # constructor. Else, we read the GDAL mask bands using
        # read_masks(), invert them and use them in constructing masked
        # arrays.

        if masked:
            enums = self.mask_flag_enums
            all_valid = all([MaskFlags.all_valid in flags for flags in enums])
            log.debug("all_valid: %s", all_valid)
            log.debug("mask_flags: %r", enums)

        # We can jump straight to _read() in some cases. We can ignore
        # the boundless flag if there's no given window.
        if not boundless or not window:

            log.debug("Jump straight to _read()")
            log.debug("Window: %r", window)

            out = self._read(indexes, out, window, dtype,
                             resampling=resampling)

            if masked:
                if all_valid:
                    mask = np.ma.nomask
                else:
                    mask = np.zeros(out.shape, 'uint8')
                    mask = ~self._read(
                        indexes, mask, window, 'uint8', masks=True,
                        resampling=resampling).astype('bool')

                kwds = {'mask': mask}
                # Set a fill value only if the read bands share a
                # single nodata value.
                if len(set(nodatavals)) == 1:
                    if nodatavals[0] is not None:
                        kwds['fill_value'] = nodatavals[0]
                out = np.ma.array(out, **kwds)

        else:
            # Compute the overlap between the dataset and the boundless window.
            overlap = ((
                max(min(window[0][0], self.height), 0),
                max(min(window[0][1], self.height), 0)), (
                max(min(window[1][0], self.width), 0),
                max(min(window[1][1], self.width), 0)))

            log.debug("Overlap: %r", overlap)

            if overlap != ((0, 0), (0, 0)):
                # Prepare a buffer.
                window_h, window_w = win_shape[-2:]
                overlap_h = overlap[0][1] - overlap[0][0]
                overlap_w = overlap[1][1] - overlap[1][0]
                scaling_h = float(out.shape[-2:][0])/window_h
                scaling_w = float(out.shape[-2:][1])/window_w
                buffer_shape = (
                        int(round(overlap_h*scaling_h)),
                        int(round(overlap_w*scaling_w)))
                data = np.zeros(win_shape[:-2] + buffer_shape, dtype)
                data = self._read(indexes, data, overlap, dtype,
                                  resampling=resampling)

                if masked:
                    mask = np.zeros(win_shape[:-2] + buffer_shape, 'uint8')
                    mask = ~self._read(
                        indexes, mask, overlap, 'uint8', masks=True,
                        resampling=resampling).astype('bool')
                    kwds = {'mask': mask}
                    if len(set(nodatavals)) == 1:
                        if nodatavals[0] is not None:
                            kwds['fill_value'] = nodatavals[0]
                    data = np.ma.array(data, **kwds)

            else:
                data = None
                if masked:
                    kwds = {'mask': True}
                    if len(set(nodatavals)) == 1:
                        if nodatavals[0] is not None:
                            kwds['fill_value'] = nodatavals[0]
                    out = np.ma.array(out, **kwds)

            if data is not None:
                # Determine where to put the data in the output window.
                data_h, data_w = buffer_shape
                roff = 0
                coff = 0
                if window[0][0] < 0:
                    roff = int(-window[0][0] * scaling_h)
                if window[1][0] < 0:
                    coff = int(-window[1][0] * scaling_w)

                for dst, src in zip(
                        out if len(out.shape) == 3 else [out],
                        data if len(data.shape) == 3 else [data]):
                    dst[roff:roff+data_h, coff:coff+data_w] = src

                if masked:
                    if not hasattr(out, 'mask'):
                        kwds = {'mask': True}
                        if len(set(nodatavals)) == 1:
                            if nodatavals[0] is not None:
                                kwds['fill_value'] = nodatavals[0]
                        out = np.ma.array(out, **kwds)

                    for dst, src in zip(
                            out.mask if len(out.shape) == 3 else [out.mask],
                            data.mask if len(data.shape) == 3 else [data.mask]):
                        dst[roff:roff+data_h, coff:coff+data_w] = src

        if return2d:
            out.shape = out.shape[1:]

        return out


    def read_masks(self, indexes=None, out=None, out_shape=None, window=None,
                   boundless=False, resampling=Resampling.nearest):
        """Read raster band masks as a multidimensional array

        Parameters
        ----------
        indexes : list of ints or a single int, optional
            If `indexes` is a list, the result is a 3D array, but is
            a 2D array if it is a band index number.

        out : numpy ndarray, optional
            As with Numpy ufuncs, this is an optional reference to an
            output array with the same dimensions and shape into which
            data will be placed.

            *Note*: the method's return value may be a view on this
            array. In other words, `out` is likely to be an
            incomplete representation of the method's results.

            Cannot combine with `out_shape`.

        out_shape : tuple, optional
            A tuple describing the output array's shape.  Allows for decimated
            reads without constructing an output Numpy array.

            Cannot combined with `out`.

        window : a pair (tuple) of pairs of ints, optional
            The optional `window` argument is a 2 item tuple. The first
            item is a tuple containing the indexes of the rows at which
            the window starts and stops and the second is a tuple
            containing the indexes of the columns at which the window
            starts and stops. For example, ((0, 2), (0, 2)) defines
            a 2x2 window at the upper left of the raster dataset.

        boundless : bool, optional (default `False`)
            If `True`, windows that extend beyond the dataset's extent
            are permitted and partially or completely filled arrays will
            be returned as appropriate.

        Returns
        -------
        Numpy ndarray or a view on a Numpy ndarray

        Note: as with Numpy ufuncs, an object is returned even if you
        use the optional `out` argument and the return value shall be
        preferentially used by callers.
        """

        return2d = False
        if indexes is None:
            indexes = self.indexes
        elif isinstance(indexes, int):
            indexes = [indexes]
            return2d = True
            if out is not None and out.ndim == 2:
                out.shape = (1,) + out.shape
        if not indexes:
            raise ValueError("No indexes to read")

        # Get the natural shape of the read window, boundless or not.
        win_shape = (len(indexes),)
        if window:
            if boundless:
                win_shape += (
                        window[0][1]-window[0][0], window[1][1]-window[1][0])
            else:
                w = windows.evaluate(window, self.height, self.width)
                minr = min(max(w[0][0], 0), self.height)
                maxr = max(0, min(w[0][1], self.height))
                minc = min(max(w[1][0], 0), self.width)
                maxc = max(0, min(w[1][1], self.width))
                win_shape += (maxr - minr, maxc - minc)
                window = ((minr, maxr), (minc, maxc))
        else:
            win_shape += self.shape

        dtype = 'uint8'

        if out is not None and out_shape is not None:
            raise ValueError("out and out_shape are exclusive")
        elif out_shape is not None:
            if len(out_shape) == 2:
                out_shape = (1,) + out_shape
            out = np.zeros(out_shape, 'uint8')

        if out is not None:
            if out.dtype != np.dtype(dtype):
                raise ValueError(
                    "the out array's dtype '%s' does not match '%s'"
                    % (out.dtype, dtype))
            if out.shape[0] != win_shape[0]:
                raise ValueError(
                    "'out' shape %s does not match window shape %s" %
                    (out.shape, win_shape))
        else:
            out = np.zeros(win_shape, 'uint8')


        # We can jump straight to _read() in some cases. We can ignore
        # the boundless flag if there's no given window.
        if not boundless or not window:
            out = self._read(indexes, out, window, dtype, masks=True,
                             resampling=resampling)

        else:
            # Compute the overlap between the dataset and the boundless window.
            overlap = ((
                max(min(window[0][0], self.height), 0),
                max(min(window[0][1], self.height), 0)), (
                max(min(window[1][0], self.width), 0),
                max(min(window[1][1], self.width), 0)))

            if overlap != ((0, 0), (0, 0)):
                # Prepare a buffer.
                window_h, window_w = win_shape[-2:]
                overlap_h = overlap[0][1] - overlap[0][0]
                overlap_w = overlap[1][1] - overlap[1][0]
                scaling_h = float(out.shape[-2:][0])/window_h
                scaling_w = float(out.shape[-2:][1])/window_w
                buffer_shape = (int(overlap_h*scaling_h), int(overlap_w*scaling_w))
                data = np.zeros(win_shape[:-2] + buffer_shape, 'uint8')
                data = self._read(indexes, data, overlap, dtype, masks=True,
                                  resampling=resampling)
            else:
                data = None

            if data is not None:
                # Determine where to put the data in the output window.
                data_h, data_w = data.shape[-2:]
                roff = 0
                coff = 0
                if window[0][0] < 0:
                    roff = int(window_h*scaling_h) - data_h
                if window[1][0] < 0:
                    coff = int(window_w*scaling_w) - data_w
                for dst, src in zip(
                        out if len(out.shape) == 3 else [out],
                        data if len(data.shape) == 3 else [data]):
                    dst[roff:roff+data_h, coff:coff+data_w] = src

        if return2d:
            out.shape = out.shape[1:]

        return out


    def _read(self, indexes, out, window, dtype, masks=False,
              resampling=Resampling.nearest):
        """Read raster bands as a multidimensional array

        If `indexes` is a list, the result is a 3D array, but
        is a 2D array if it is a band index number.

        Optional `out` argument is a reference to an output array with the
        same dimensions and shape.

        See `read_band` for usage of the optional `window` argument.

        The return type will be either a regular NumPy array, or a masked
        NumPy array depending on the `masked` argument. The return type is
        forced if either `True` or `False`, but will be chosen if `None`.
        For `masked=None` (default), the array will be the same type as
        `out` (if used), or will be masked if any of the nodatavals are
        not `None`.
        """
        cdef int aix, bidx, indexes_count
        cdef double height, width, xoff, yoff
        cdef int retval = 0
        cdef GDALDatasetH dataset = NULL

        if out is None:
            raise ValueError("An output array is required.")

        dataset = self.handle()

        # Turning the read window into GDAL offsets and lengths is
        # the job of _read().
        if window:
            window = windows.evaluate(window, self.height, self.width)

            log.debug("Eval'd window: %r", window)

            yoff = window[0][0]
            xoff = window[1][0]
            height = window[0][1] - yoff
            width = window[1][1] - xoff

        else:
            xoff = yoff = <int>0
            width = <int>self.width
            height = <int>self.height

        log.debug(
            "IO window xoff=%s yoff=%s width=%s height=%s",
            xoff, yoff, width, height)

        # Call io_multi* functions with C type args so that they
        # can release the GIL.
        indexes_arr = np.array(indexes, dtype=int)
        indexes_count = <int>indexes_arr.shape[0]

        if masks:
            # Warn if nodata attribute is shadowing an alpha band.
            if self.count == 4 and self.colorinterp(4) == ColorInterp.alpha:
                for flags in self.mask_flag_enums:
                    if MaskFlags.nodata in flags:
                        warnings.warn(NodataShadowWarning())

            retval = io_multi_mask(
                            self._hds, 0, xoff, yoff, width, height,
                            out, indexes_arr, resampling=resampling)

        else:
            retval = io_multi_band(self._hds, 0, xoff, yoff, width, height,
                                  out, indexes_arr, resampling=resampling)

{'fraction_non_alphanumeric': 0.08925602217011298, 'fraction_numerical': 0.010360264335962481, 'mean_word_length': 1.603274139844617, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5064203', 'n_tokens_mistral': 13621, 'n_tokens_neox': 12842, 'n_words': 4520}