Enable windowed read and write.

This is done by extending the existing read_band() and write_band()
methods. The windows that can be read and written most efficiently
are the 7nes that come from the new block_windows property.

Closes #6.

CHANGES.txt

@@ -6,6 +6,8 @@ Next
 - Drop six dependency (#9)
 - Add crs_wkt attribute (#10).
 - Add bounds attribute and ul() method (#11).
+- Add block_windows property (#7).
+- Enable windowed reads and writes (#6).
 
 0.2 (2013-11-24)
 ----------------

rasterio/_io.pyx

@@ -28,73 +28,87 @@ cdef void register():
     registered = 1
 
 cdef int io_ubyte(
-        void *hband, 
+        void *hband,
         int mode,
+        int xoff,
+        int yoff,
         int width, 
         int height, 
         np.ndarray[DTYPE_UBYTE_t, ndim=2, mode='c'] buffer):
     return _gdal.GDALRasterIO(
-        hband, mode, 0, 0, width, height,
+        hband, mode, xoff, yoff, width, height,
         &buffer[0, 0], buffer.shape[1], buffer.shape[0], 1, 0, 0)
 
 cdef int io_uint16(
         void *hband, 
         int mode,
+        int xoff,
+        int yoff,
         int width, 
         int height, 
         np.ndarray[DTYPE_UINT16_t, ndim=2, mode='c'] buffer):
     return _gdal.GDALRasterIO(
-        hband, mode, 0, 0, width, height,
+        hband, mode, xoff, yoff, width, height,
         &buffer[0, 0], buffer.shape[1], buffer.shape[0], 2, 0, 0)
 
 cdef int io_int16(
         void *hband, 
         int mode,
+        int xoff,
+        int yoff,
         int width, 
         int height, 
         np.ndarray[DTYPE_INT16_t, ndim=2, mode='c'] buffer):
     return _gdal.GDALRasterIO(
-        hband, mode, 0, 0, width, height,
+        hband, mode, xoff, yoff, width, height,
         &buffer[0, 0], buffer.shape[1], buffer.shape[0], 3, 0, 0)
 
 cdef int io_uint32(
         void *hband, 
         int mode,
+        int xoff,
+        int yoff,
         int width, 
         int height, 
         np.ndarray[DTYPE_UINT32_t, ndim=2, mode='c'] buffer):
     return _gdal.GDALRasterIO(
-        hband, mode, 0, 0, width, height,
+        hband, mode, xoff, yoff, width, height,
         &buffer[0, 0], buffer.shape[1], buffer.shape[0], 4, 0, 0)
 
 cdef int io_int32(
         void *hband, 
         int mode,
+        int xoff,
+        int yoff,
         int width, 
         int height, 
         np.ndarray[DTYPE_INT32_t, ndim=2, mode='c'] buffer):
     return _gdal.GDALRasterIO(
-        hband, mode, 0, 0, width, height,
+        hband, mode, xoff, yoff, width, height,
         &buffer[0, 0], buffer.shape[1], buffer.shape[0], 5, 0, 0)
 
 cdef int io_float32(
         void *hband, 
         int mode,
+        int xoff,
+        int yoff,
         int width, 
         int height, 
         np.ndarray[DTYPE_FLOAT32_t, ndim=2, mode='c'] buffer):
     return _gdal.GDALRasterIO(
-        hband, mode, 0, 0, width, height,
+        hband, mode, xoff, yoff, width, height,
         &buffer[0, 0], buffer.shape[1], buffer.shape[0], 6, 0, 0)
 
 cdef int io_float64(
         void *hband,
         int mode,
+        int xoff,
+        int yoff,
         int width, 
         int height, 
         np.ndarray[DTYPE_FLOAT64_t, ndim=2, mode='c'] buffer):
     return _gdal.GDALRasterIO(
-        hband, mode, 0, 0, width, height,
+        hband, mode, xoff, yoff, width, height,
         &buffer[0, 0], buffer.shape[1], buffer.shape[0], 7, 0, 0)
 
 
@@ -356,11 +370,14 @@ cdef class RasterReader:
         return self._transform
     transform = property(get_transform)
 
-    def read_band(self, bidx, out=None):
+    def read_band(self, bidx, out=None, window=None):
         """Read the `bidx` band into an `out` array if provided, otherwise
         return a new array.
-        
+
         Band indexes begin with 1: read_band(1) returns the first band.
+
+        The optional `window` argument takes 4 item tuples specifying
+        (xoffset, yoffset, width, height) of a raster subset.
         """
         if bidx not in self.indexes:
             raise IndexError("band index out of range")
@@ -369,26 +386,44 @@ 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")
-        dtype = self.dtypes[i]
-        if out is None:
-            out = np.zeros(self.shape, dtype)
+        if window and out is not None and out.shape[::-1] != window[2:]:
+            raise ValueError("output and window dimensions do not match")
+
         cdef void *hband = _gdal.GDALGetRasterBand(self._hds, bidx)
         if hband is NULL:
             raise ValueError("NULL band")
+
+        dtype = self.dtypes[i]
+        if out is None:
+            out_shape = window and window[2:][::-1] or self.shape
+            out = np.zeros(out_shape, dtype)
+        if window:
+            xoff, yoff, width, height = window
+        else:
+            xoff = yoff = 0
+            width = self.width
+            height = self.height
         if dtype == dtypes.ubyte:
-            retval = io_ubyte(hband, 0, self.width, self.height, out)
+            retval = io_ubyte(
+                hband, 0, xoff, yoff, width, height, out)
         elif dtype == dtypes.uint16:
-            retval = io_uint16(hband, 0, self.width, self.height, out)
+            retval = io_uint16(
+                hband, 0, xoff, yoff, width, height, out)
         elif dtype == dtypes.int16:
-            retval = io_int16(hband, 0, self.width, self.height, out)
+            retval = io_int16(
+                hband, 0, xoff, yoff, width, height, out)
         elif dtype == dtypes.uint32:
-            retval = io_uint32(hband, 0, self.width, self.height, out)
+            retval = io_uint32(
+                hband, 0, xoff, yoff, width, height, out)
         elif dtype == dtypes.int32:
-            retval = io_int32(hband, 0, self.width, self.height, out)
+            retval = io_int32(
+                hband, 0, xoff, yoff, width, height, out)
         elif dtype == dtypes.float32:
-            retval = io_float32(hband, 0, self.width, self.height, out)
+            retval = io_float32(
+                hband, 0, xoff, yoff, width, height, out)
         elif dtype == dtypes.float64:
-            retval = io_float64(hband, 0, self.width, self.height, out)
+            retval = io_float64(
+                hband, 0, xoff, yoff, width, height, out)
         else:
             raise ValueError("Invalid dtype")
         # TODO: handle errors (by retval).
@@ -543,10 +578,13 @@ cdef class RasterUpdater(RasterReader):
 
     transform = property(get_transform, write_transform)
 
-    def write_band(self, bidx, src):
+    def write_band(self, bidx, src, window=None):
         """Write the src array into the `bidx` band.
         
         Band indexes begin with 1: read_band(1) returns the first band.
+
+        The optional `window` argument takes 4 item tuples specifying
+        (xoffset, yoffset, width, height) of a raster subset to write into.
         """
         if bidx not in self.indexes:
             raise IndexError("band index out of range")
@@ -555,24 +593,41 @@ 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")
-        dtype = self.dtypes[i]
+        if window and src is not None and src.shape[::-1] != window[2:]:
+            raise ValueError("source and window dimensions do not match")
+
         cdef void *hband = _gdal.GDALGetRasterBand(self._hds, bidx)
         if hband is NULL:
             raise ValueError("NULL band")
+
+        if window:
+            xoff, yoff, width, height = window
+        else:
+            xoff = yoff = 0
+            width = self.width
+            height = self.height
+        dtype = self.dtypes[i]
         if dtype == dtypes.ubyte:
-            retval = io_ubyte(hband, 1, self.width, self.height, src)
+            retval = io_ubyte(
+                hband, 1, xoff, yoff, width, height, src)
         elif dtype == dtypes.uint16:
-            retval = io_uint16(hband, 1, self.width, self.height, src)
+            retval = io_uint16(
+                hband, 1, xoff, yoff, width, height, src)
         elif dtype == dtypes.int16:
-            retval = io_int16(hband, 1, self.width, self.height, src)
+            retval = io_int16(
+                hband, 1, xoff, yoff, width, height, src)
         elif dtype == dtypes.uint32:
-            retval = io_uint32(hband, 1, self.width, self.height, src)
+            retval = io_uint32(
+                hband, 1, xoff, yoff, width, height, src)
         elif dtype == dtypes.int32:
-            retval = io_int32(hband, 1, self.width, self.height, src)
+            retval = io_int32(
+                hband, 1, xoff, yoff, width, height, src)
         elif dtype == dtypes.float32:
-            retval = io_float32(hband, 1, self.width, self.height, src)
+            retval = io_float32(
+                hband, 1, xoff, yoff, width, height, src)
         elif dtype == dtypes.float64:
-            retval = io_float64(hband, 1, self.width, self.height, src)
+            retval = io_float64(
+                hband, 1, xoff, yoff, width, height, src)
         else:
             raise ValueError("Invalid dtype")
         # TODO: handle errors (by retval).

rasterio/tests/test_blocks.py

@@ -1,7 +1,17 @@
+import logging
+import os.path
 import unittest
+import shutil
+import subprocess
+import sys
+import tempfile
+
+import numpy
 
 import rasterio
 
+logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
+
 class RasterBlocksTest(unittest.TestCase):
     def test_blocks(self):
         with rasterio.open('rasterio/tests/data/RGB.byte.tif') as s:
@@ -14,4 +24,37 @@ def test_blocks(self):
             self.assertEqual(second, (0, 3, 791, 3))
             last = list(windows)[~0]
             self.assertEqual(last, (0, 717, 791, 1))
+    def test_block_coverage(self):
+        with rasterio.open('rasterio/tests/data/RGB.byte.tif') as s:
+            self.assertEqual(
+                s.width*s.height,
+                sum(b[2]*b[3] for b in s.block_windows(1)))
+
+class WindowReadTest(unittest.TestCase):
+    def test_read_window(self):
+        with rasterio.open('rasterio/tests/data/RGB.byte.tif') as s:
+            windows = s.block_windows(1)
+            first_window = next(windows)
+            first_block = s.read_band(1, window=first_window)
+            self.assertEqual(first_block.dtype, rasterio.ubyte)
+            self.assertEqual(first_block.shape[::-1], first_window[2:])
+
+class WindowWriteTest(unittest.TestCase):
+    def setUp(self):
+        self.tempdir = tempfile.mkdtemp()
+    def tearDown(self):
+        shutil.rmtree(self.tempdir)
+    def test_write_window(self):
+        name = os.path.join(self.tempdir, "test_write_window.tif")
+        a = numpy.ones((50, 50), dtype=rasterio.ubyte) * 127
+        with rasterio.open(
+                name, 'w', 
+                driver='GTiff', width=100, height=100, count=1, 
+                dtype=a.dtype) as s:
+            s.write_band(1, a, window=(30, 10, 50, 50))
+        info = subprocess.check_output(["gdalinfo", "-stats", name])
+        self.assert_(
+            "Minimum=0.000, Maximum=127.000, "
+            "Mean=31.750, StdDev=54.993" in info.decode('utf-8'),
+            info)