diff --git a/docs/windowed-rw.rst b/docs/windowed-rw.rst index 5da0b802..d8fa1864 100644 --- a/docs/windowed-rw.rst +++ b/docs/windowed-rw.rst @@ -123,3 +123,82 @@ And the result: :width: 500 :height: 300 +Advanced windows +---------------- + +Since windows are like slices, you can also use negative numbers in rasterio +windows. + +.. code-block:: python + + ((-4, None), (-4, None)) + +specifies a 4 x 4 rectangular subset with upper left at 4 rows to the left of +and 4 columns above the lower right corner of the dataset and extending to the +lower right corner of the dataset. + +Below is an example of reading a raster subset and then writing it into a +larger subset that is defined relative to the lower right corner of the +destination dataset. + +.. code-block:: python + + read_window = (350, 410), (350, 450) + + with rasterio.open('rasterio/tests/data/RGB.byte.tif') as src: + r = src.read_band(1, window=read_window) + g = src.read_band(2, window=read_window) + b = src.read_band(3, window=read_window) + + write_window = (-240, None), (-400, None) + + with rasterio.open( + '/tmp/example2.tif', 'w', + driver='GTiff', width=500, height=300, count=3, + dtype=r.dtype) as dst: + dst.write_band(1, r, window=write_window) + dst.write_band(2, g, window=write_window) + dst.write_band(3, b, window=write_window) + +This example also demonstrates decimation. + +.. image:: http://farm3.staticflickr.com/2827/11378772013_c8ab540f21_o_d.png + :width: 500 + :height: 300 + +Blocks +------ + +Raster datasets are generally composed of multiple blocks of data and +windowed reads and writes are most efficient when the windows match the +dataset's own block structure. When a file is opened to read, the shape +of blocks for any band can be had from the block_shapes property. + +.. code-block:: pycon + + >>> with rasterio.open('rasterio/tests/data/RGB.byte.tif') as src: + ... for i, shape in enumerate(src.block_shapes, 1): + ... print(i, shape) + ... + (1, (3, 791)) + (2, (3, 791)) + (3, (3, 791)) + + +The block windows themselves can be had from the block_windows function. + +.. code-block:: pycon + + >>> with rasterio.open('rasterio/tests/data/RGB.byte.tif') as src: + ... for ji, window in src.block_windows(1): + ... print(ji, window) + ... + ((0, 0), ((0, 3), (0, 791))) + ((1, 0), ((3, 6), (0, 791))) + ... + +This function returns an iterator that yields a pair of values. The second is +a window tuple that can be used in calls to read_band or write_band. The first +is the pair of row and column indexes of this block within all blocks of the +dataset. + diff --git a/rasterio/_io.pyx b/rasterio/_io.pyx index 8e6d7c7c..08e4c748 100644 --- a/rasterio/_io.pyx +++ b/rasterio/_io.pyx @@ -228,7 +228,6 @@ cdef class RasterReader: cdef char *proj_c = NULL if self._hds is NULL: raise ValueError("Null dataset") - #cdef const char *wkt = _gdal.GDALGetProjectionRef(self._hds) cdef void *osr = _gdal.OSRNewSpatialReference( _gdal.GDALGetProjectionRef(self._hds)) log.debug("Got coordinate system") @@ -456,8 +455,6 @@ cdef class RasterReader: raise ValueError("can't read closed raster file") if out is not None and out.dtype != self.dtypes[i]: raise ValueError("band and output array dtypes do not match") - #if window and out is not None and out.shape != window_shape(window): - # raise ValueError("output and window dimensions do not match") cdef void *hband = _gdal.GDALGetRasterBand(self._hds, bidx) if hband is NULL: @@ -465,7 +462,10 @@ cdef class RasterReader: dtype = self.dtypes[i] if out is None: - out_shape = window and window_shape(window) or self.shape + out_shape = ( + window + and window_shape(window, self.height, self.width) + or self.shape) out = np.zeros(out_shape, dtype) if window: window = eval_window(window, self.height, self.width) @@ -670,8 +670,6 @@ cdef class RasterUpdater(RasterReader): raise ValueError("can't read closed raster file") if src is not None and src.dtype != self.dtypes[i]: raise ValueError("band and srcput array dtypes do not match") - #if window and src is not None and src.shape != window_shape(window): - # raise ValueError("source and window dimensions do not match") cdef void *hband = _gdal.GDALGetRasterBand(self._hds, bidx) if hband is NULL: