diff --git a/CHANGES.txt b/CHANGES.txt
index 4fd3c31a..b895f20a 100644
--- a/CHANGES.txt
+++ b/CHANGES.txt
@@ -1,6 +1,10 @@
 Changes
 =======
 
+0.22.0 (?)
+----------
+- Return masked arrays in the boundless read case (#338).
+
 0.21.0 (2015-04-22)
 -------------------
 - New rio-mask command (#323).
diff --git a/rasterio/__init__.py b/rasterio/__init__.py
index 323bfc53..bad24e5a 100644
--- a/rasterio/__init__.py
+++ b/rasterio/__init__.py
@@ -23,7 +23,7 @@ from rasterio import _err, coords, enums
 
 __all__ = [
     'band', 'open', 'drivers', 'copy', 'pad']
-__version__ = "0.21.0"
+__version__ = "0.22.0"
 
 log = logging.getLogger('rasterio')
 class NullHandler(logging.Handler):
diff --git a/rasterio/_io.pyx b/rasterio/_io.pyx
index 8e1fe174..989ca52b 100644
--- a/rasterio/_io.pyx
+++ b/rasterio/_io.pyx
@@ -323,11 +323,12 @@ cdef int io_multi_cint16(
     cdef int i, j, k
     cdef np.int16_t real, imag
 
-    buf = np.empty(
+    buf = np.zeros(
             (out.shape[0], 2*out.shape[2]*out.shape[1]), 
             dtype=np.int16)
     cdef np.int16_t[:, :] buf_view = buf
 
+    
     with nogil:
         bandmap = <int *>_gdal.CPLMalloc(count*sizeof(int))
         for i in range(count):
@@ -715,8 +716,6 @@ cdef class RasterReader(_base.DatasetReader):
                 raise ValueError(
                     "'out' shape %s does not match window shape %s" %
                     (out.shape, win_shape))
-            #if masked is None:
-            #    masked = hasattr(out, 'mask')
 
         # Masking
         # -------
@@ -790,8 +789,25 @@ cdef class RasterReader(_base.DatasetReader):
                 data = np.empty(win_shape[:-2] + buffer_shape, dtype)
                 data = self._read(indexes, data, overlap, dtype)
 
+                if masked:
+                    mask = np.empty(win_shape[:-2] + buffer_shape, 'uint8')
+                    mask = ~self._read(
+                        indexes, mask, overlap, 'uint8', masks=True
+                        ).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.
@@ -808,6 +824,19 @@ cdef class RasterReader(_base.DatasetReader):
                         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:]
 
diff --git a/tests/test_read_boundless.py b/tests/test_read_boundless.py
index e6c2a1f5..dc29c36e 100644
--- a/tests/test_read_boundless.py
+++ b/tests/test_read_boundless.py
@@ -16,6 +16,7 @@ def test_read_boundless_natural_extent():
         assert abs(data[0].mean() - src.read(1).mean()) < 0.0001
         assert data.any()
 
+
 def test_read_boundless_beyond():
     with rasterio.open('tests/data/RGB.byte.tif') as src:
         data = src.read(window=((-200, -100), (-200, -100)), boundless=True)
@@ -49,3 +50,22 @@ def test_read_boundless_resample():
         assert data.shape == (3, 800, 800)
         assert data.any()
         assert data[0,798,798] == 13
+
+
+def test_read_boundless_masked_no_overlap():
+    with rasterio.open('tests/data/RGB.byte.tif') as src:
+        data = src.read(
+            window=((-200, -100), (-200, -100)), boundless=True, masked=True)
+        assert data.shape == (3, 100, 100)
+        assert data.mask.all()
+
+
+def test_read_boundless_masked_overlap():
+    with rasterio.open('tests/data/RGB.byte.tif') as src:
+        data = src.read(
+            window=((-200, 200), (-200, 200)), boundless=True, masked=True)
+        assert data.shape == (3, 400, 400)
+        assert data.mask.any()
+        assert not data.mask.all()
+        assert data.mask[0,399,399] == False
+        assert data.mask[0,0,0] == True