We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
rasterio.features.rasterize crashes (segfault) when out_shape has more than two dimensions
I just spent an embarassing amount of time trying to debug a segfault in rasterize, and the cause ended up being that I passed `out_shape` as a 3-d array (e.g. `(1, 21, 18)`), when it should be 2-d (e.g. `(21, 18)`).

The error-causing shape was the same as an array I read from rasterio itself, using dataset.read() with the `window` kwarg set. That is, I had a `band = dataset.read(window=some_window)`, and passed `band.shape` to `rasterize()`.

Should `rasterize()` ever accept a shape with more than two elements? If not, maybe rasterio could fail early with a descriptive message, like it already does if you pass a wrong `dtype`, for example.


--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `rasterio/features.py`
Content:
```
1 """Functions for working with features in a raster dataset."""
2 
3 
4 import logging
5 
6 import numpy as np
7 
8 from rasterio._features import _shapes, _sieve, _rasterize, _bounds
9 from rasterio.dtypes import validate_dtype, can_cast_dtype, get_minimum_dtype
10 from rasterio.env import ensure_env
11 from rasterio.transform import IDENTITY, guard_transform
12 
13 
14 log = logging.getLogger(__name__)
15 
16 
17 @ensure_env
18 def geometry_mask(
19         geometries,
20         out_shape,
21         transform,
22         all_touched=False,
23         invert=False):
24     """Create a mask from shapes.
25 
26     By default, mask is intended for use as a
27     numpy mask, where pixels that overlap shapes are False.
28 
29     Parameters
30     ----------
31     geometries : iterable over geometries (GeoJSON-like objects)
32     out_shape : tuple or list
33         Shape of output numpy ndarray.
34     transform : Affine transformation object
35         Transformation from pixel coordinates of `image` to the
36         coordinate system of the input `shapes`. See the `transform`
37         property of dataset objects.
38     all_touched : boolean, optional
39         If True, all pixels touched by geometries will be burned in.  If
40         false, only pixels whose center is within the polygon or that
41         are selected by Bresenham's line algorithm will be burned in.
42     invert: boolean, optional
43         If True, mask will be True for pixels that overlap shapes.
44         False by default.
45 
46     Returns
47     -------
48     out : numpy ndarray of type 'bool'
49         Result
50     """
51     fill, mask_value = (0, 1) if invert else (1, 0)
52 
53     return rasterize(
54         geometries,
55         out_shape=out_shape,
56         transform=transform,
57         all_touched=all_touched,
58         fill=fill,
59         default_value=mask_value).astype('bool')
60 
61 
62 @ensure_env
63 def shapes(image, mask=None, connectivity=4, transform=IDENTITY):
64     """Yield (polygon, value for each set of adjacent pixels of the same value.
65 
66     Parameters
67     ----------
68     image : numpy ndarray or rasterio Band object
69         (RasterReader, bidx namedtuple).
70         Data type must be one of rasterio.int16, rasterio.int32,
71         rasterio.uint8, rasterio.uint16, or rasterio.float32.
72     mask : numpy ndarray or rasterio Band object, optional
73         Values of False or 0 will be excluded from feature generation
74         Must evaluate to bool (rasterio.bool_ or rasterio.uint8)
75     connectivity : int, optional
76         Use 4 or 8 pixel connectivity for grouping pixels into features
77     transform : Affine transformation, optional
78         If not provided, feature coordinates will be generated based on pixel
79         coordinates
80 
81     Yields
82     -------
83     tuple
84         A pair of (polygon, value) for each feature found in the image.
85         Polygons are GeoJSON-like dicts and the values are the associated value
86         from the image, in the data type of the image.
87         Note: due to floating point precision issues, values returned from a
88         floating point image may not exactly match the original values.
89 
90     Notes
91     -----
92     The amount of memory used by this algorithm is proportional to the number
93     and complexity of polygons produced.  This algorithm is most appropriate
94     for simple thematic data.  Data with high pixel-to-pixel variability, such
95     as imagery, may produce one polygon per pixel and consume large amounts of
96     memory.
97 
98     """
99     transform = guard_transform(transform)
100     for s, v in _shapes(image, mask, connectivity, transform.to_gdal()):
101         yield s, v
102 
103 
104 @ensure_env
105 def sieve(image, size, out=None, mask=None, connectivity=4):
106     """Replace small polygons in `image` with value of their largest neighbor.
107 
108     Polygons are found for each set of neighboring pixels of the same value.
109 
110     Parameters
111     ----------
112     image : numpy ndarray or rasterio Band object
113         (RasterReader, bidx namedtuple)
114         Must be of type rasterio.int16, rasterio.int32, rasterio.uint8,
115         rasterio.uint16, or rasterio.float32
116     size : int
117         minimum polygon size (number of pixels) to retain.
118     out : numpy ndarray, optional
119         Array of same shape and data type as `image` in which to store results.
120     mask : numpy ndarray or rasterio Band object, optional
121         Values of False or 0 will be excluded from feature generation
122         Must evaluate to bool (rasterio.bool_ or rasterio.uint8)
123     connectivity : int, optional
124         Use 4 or 8 pixel connectivity for grouping pixels into features
125 
126     Returns
127     -------
128     out : numpy ndarray
129         Result
130 
131     Notes
132     -----
133     GDAL only supports values that can be cast to 32-bit integers for this
134     operation.
135 
136     The amount of memory used by this algorithm is proportional to the number
137     and complexity of polygons found in the image.  This algorithm is most
138     appropriate for simple thematic data.  Data with high pixel-to-pixel
139     variability, such as imagery, may produce one polygon per pixel and consume
140     large amounts of memory.
141 
142     """
143 
144     if out is None:
145         out = np.zeros(image.shape, image.dtype)
146     _sieve(image, size, out, mask, connectivity)
147     return out
148 
149 
150 @ensure_env
151 def rasterize(
152         shapes,
153         out_shape=None,
154         fill=0,
155         out=None,
156         transform=IDENTITY,
157         all_touched=False,
158         default_value=1,
159         dtype=None):
160     """Return an image array with input geometries burned in.
161 
162     Parameters
163     ----------
164     shapes : iterable of (geometry, value) pairs or iterable over
165         geometries. `geometry` can either be an object that implements
166         the geo interface or GeoJSON-like object.
167     out_shape : tuple or list
168         Shape of output numpy ndarray.
169     fill : int or float, optional
170         Used as fill value for all areas not covered by input
171         geometries.
172     out : numpy ndarray, optional
173         Array of same shape and data type as `image` in which to store
174         results.
175     transform : Affine transformation object, optional
176         Transformation from pixel coordinates of `image` to the
177         coordinate system of the input `shapes`. See the `transform`
178         property of dataset objects.
179     all_touched : boolean, optional
180         If True, all pixels touched by geometries will be burned in.  If
181         false, only pixels whose center is within the polygon or that
182         are selected by Bresenham's line algorithm will be burned in.
183     default_value : int or float, optional
184         Used as value for all geometries, if not provided in `shapes`.
185     dtype : rasterio or numpy data type, optional
186         Used as data type for results, if `out` is not provided.
187 
188     Returns
189     -------
190     out : numpy ndarray
191         Results
192 
193     Notes
194     -----
195     Valid data types for `fill`, `default_value`, `out`, `dtype` and
196     shape values are rasterio.int16, rasterio.int32, rasterio.uint8,
197     rasterio.uint16, rasterio.uint32, rasterio.float32,
198     rasterio.float64.
199 
200     """
201     valid_dtypes = (
202         'int16', 'int32', 'uint8', 'uint16', 'uint32', 'float32', 'float64'
203     )
204 
205     def format_invalid_dtype(param):
206         return '{0} dtype must be one of: {1}'.format(
207             param, ', '.join(valid_dtypes)
208         )
209 
210     def format_cast_error(param, dtype):
211         return '{0} cannot be cast to specified dtype: {1}'.format(param, dtype)
212 
213     if fill != 0:
214         fill_array = np.array([fill])
215         if not validate_dtype(fill_array, valid_dtypes):
216             raise ValueError(format_invalid_dtype('fill'))
217 
218         if dtype is not None and not can_cast_dtype(fill_array, dtype):
219             raise ValueError(format_cast_error('fill', dtype))
220 
221     if default_value != 1:
222         default_value_array = np.array([default_value])
223         if not validate_dtype(default_value_array, valid_dtypes):
224             raise ValueError(format_invalid_dtype('default_value'))
225 
226         if dtype is not None and not can_cast_dtype(default_value_array, dtype):
227             raise ValueError(format_cast_error('default_vaue', dtype))
228 
229     if dtype is not None and np.dtype(dtype).name not in valid_dtypes:
230         raise ValueError(format_invalid_dtype('dtype'))
231 
232     valid_shapes = []
233     shape_values = []
234     for index, item in enumerate(shapes):
235         if isinstance(item, (tuple, list)):
236             geom, value = item
237         else:
238             geom = item
239             value = default_value
240         geom = getattr(geom, '__geo_interface__', None) or geom
241 
242         # not isinstance(geom, dict) or
243         if 'type' in geom or 'coordinates' in geom:
244             valid_shapes.append((geom, value))
245             shape_values.append(value)
246 
247         else:
248             raise ValueError(
249                 'Invalid geometry object at index {0}'.format(index)
250             )
251 
252     if not valid_shapes:
253         raise ValueError('No valid geometry objects found for rasterize')
254 
255     shape_values = np.array(shape_values)
256 
257     if not validate_dtype(shape_values, valid_dtypes):
258         raise ValueError(format_invalid_dtype('shape values'))
259 
260     if dtype is None:
261         dtype = get_minimum_dtype(np.append(shape_values, fill))
262 
263     elif not can_cast_dtype(shape_values, dtype):
264         raise ValueError(format_cast_error('shape values', dtype))
265 
266     if out is not None:
267         if np.dtype(out.dtype).name not in valid_dtypes:
268             raise ValueError(format_invalid_dtype('out'))
269 
270         if not can_cast_dtype(shape_values, out.dtype):
271             raise ValueError(format_cast_error('shape values', out.dtype.name))
272 
273     elif out_shape is not None:
274         out = np.empty(out_shape, dtype=dtype)
275         out.fill(fill)
276 
277     else:
278         raise ValueError('Either an out_shape or image must be provided')
279 
280     transform = guard_transform(transform)
281     _rasterize(valid_shapes, out, transform.to_gdal(), all_touched)
282     return out
283 
284 
285 def bounds(geometry):
286     """Return a (minx, miny, maxx, maxy) bounding box.
287 
288     From Fiona 1.4.8. Modified to return bbox from geometry if available.
289 
290     Parameters
291     ----------
292     geometry: GeoJSON-like feature, feature collection, or geometry.
293 
294     Returns
295     -------
296     tuple
297         Bounding box: (minx, miny, maxx, maxy)
298     """
299     if 'bbox' in geometry:
300         return tuple(geometry['bbox'])
301 
302     geom = geometry.get('geometry') or geometry
303     return _bounds(geom)
304
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/rasterio/features.py b/rasterio/features.py
--- a/rasterio/features.py
+++ b/rasterio/features.py
@@ -164,7 +164,7 @@
     shapes : iterable of (geometry, value) pairs or iterable over
         geometries. `geometry` can either be an object that implements
         the geo interface or GeoJSON-like object.
-    out_shape : tuple or list
+    out_shape : tuple or list with 2 integers
         Shape of output numpy ndarray.
     fill : int or float, optional
         Used as fill value for all areas not covered by input
@@ -271,6 +271,10 @@
             raise ValueError(format_cast_error('shape values', out.dtype.name))
 
     elif out_shape is not None:
+
+        if len(out_shape) != 2:
+            raise ValueError('Invalid out_shape, must be 2D')
+
         out = np.empty(out_shape, dtype=dtype)
         out.fill(fill)

{"golden_diff": "diff --git a/rasterio/features.py b/rasterio/features.py\n--- a/rasterio/features.py\n+++ b/rasterio/features.py\n@@ -164,7 +164,7 @@\n     shapes : iterable of (geometry, value) pairs or iterable over\n         geometries. `geometry` can either be an object that implements\n         the geo interface or GeoJSON-like object.\n-    out_shape : tuple or list\n+    out_shape : tuple or list with 2 integers\n         Shape of output numpy ndarray.\n     fill : int or float, optional\n         Used as fill value for all areas not covered by input\n@@ -271,6 +271,10 @@\n             raise ValueError(format_cast_error('shape values', out.dtype.name))\n \n     elif out_shape is not None:\n+\n+        if len(out_shape) != 2:\n+            raise ValueError('Invalid out_shape, must be 2D')\n+\n         out = np.empty(out_shape, dtype=dtype)\n         out.fill(fill)\n", "issue": "rasterio.features.rasterize crashes (segfault) when out_shape has more than two dimensions\nI just spent an embarassing amount of time trying to debug a segfault in rasterize, and the cause ended up being that I passed `out_shape` as a 3-d array (e.g. `(1, 21, 18)`), when it should be 2-d (e.g. `(21, 18)`).\n\nThe error-causing shape was the same as an array I read from rasterio itself, using dataset.read() with the `window` kwarg set. That is, I had a `band = dataset.read(window=some_window)`, and passed `band.shape` to `rasterize()`.\n\nShould `rasterize()` ever accept a shape with more than two elements? If not, maybe rasterio could fail early with a descriptive message, like it already does if you pass a wrong `dtype`, for example.\n\n", "before_files": [{"content": "\"\"\"Functions for working with features in a raster dataset.\"\"\"\n\n\nimport logging\n\nimport numpy as np\n\nfrom rasterio._features import _shapes, _sieve, _rasterize, _bounds\nfrom rasterio.dtypes import validate_dtype, can_cast_dtype, get_minimum_dtype\nfrom rasterio.env import ensure_env\nfrom rasterio.transform import IDENTITY, guard_transform\n\n\nlog = logging.getLogger(__name__)\n\n\n@ensure_env\ndef geometry_mask(\n        geometries,\n        out_shape,\n        transform,\n        all_touched=False,\n        invert=False):\n    \"\"\"Create a mask from shapes.\n\n    By default, mask is intended for use as a\n    numpy mask, where pixels that overlap shapes are False.\n\n    Parameters\n    ----------\n    geometries : iterable over geometries (GeoJSON-like objects)\n    out_shape : tuple or list\n        Shape of output numpy ndarray.\n    transform : Affine transformation object\n        Transformation from pixel coordinates of `image` to the\n        coordinate system of the input `shapes`. See the `transform`\n        property of dataset objects.\n    all_touched : boolean, optional\n        If True, all pixels touched by geometries will be burned in.  If\n        false, only pixels whose center is within the polygon or that\n        are selected by Bresenham's line algorithm will be burned in.\n    invert: boolean, optional\n        If True, mask will be True for pixels that overlap shapes.\n        False by default.\n\n    Returns\n    -------\n    out : numpy ndarray of type 'bool'\n        Result\n    \"\"\"\n    fill, mask_value = (0, 1) if invert else (1, 0)\n\n    return rasterize(\n        geometries,\n        out_shape=out_shape,\n        transform=transform,\n        all_touched=all_touched,\n        fill=fill,\n        default_value=mask_value).astype('bool')\n\n\n@ensure_env\ndef shapes(image, mask=None, connectivity=4, transform=IDENTITY):\n    \"\"\"Yield (polygon, value for each set of adjacent pixels of the same value.\n\n    Parameters\n    ----------\n    image : numpy ndarray or rasterio Band object\n        (RasterReader, bidx namedtuple).\n        Data type must be one of rasterio.int16, rasterio.int32,\n        rasterio.uint8, rasterio.uint16, or rasterio.float32.\n    mask : numpy ndarray or rasterio Band object, optional\n        Values of False or 0 will be excluded from feature generation\n        Must evaluate to bool (rasterio.bool_ or rasterio.uint8)\n    connectivity : int, optional\n        Use 4 or 8 pixel connectivity for grouping pixels into features\n    transform : Affine transformation, optional\n        If not provided, feature coordinates will be generated based on pixel\n        coordinates\n\n    Yields\n    -------\n    tuple\n        A pair of (polygon, value) for each feature found in the image.\n        Polygons are GeoJSON-like dicts and the values are the associated value\n        from the image, in the data type of the image.\n        Note: due to floating point precision issues, values returned from a\n        floating point image may not exactly match the original values.\n\n    Notes\n    -----\n    The amount of memory used by this algorithm is proportional to the number\n    and complexity of polygons produced.  This algorithm is most appropriate\n    for simple thematic data.  Data with high pixel-to-pixel variability, such\n    as imagery, may produce one polygon per pixel and consume large amounts of\n    memory.\n\n    \"\"\"\n    transform = guard_transform(transform)\n    for s, v in _shapes(image, mask, connectivity, transform.to_gdal()):\n        yield s, v\n\n\n@ensure_env\ndef sieve(image, size, out=None, mask=None, connectivity=4):\n    \"\"\"Replace small polygons in `image` with value of their largest neighbor.\n\n    Polygons are found for each set of neighboring pixels of the same value.\n\n    Parameters\n    ----------\n    image : numpy ndarray or rasterio Band object\n        (RasterReader, bidx namedtuple)\n        Must be of type rasterio.int16, rasterio.int32, rasterio.uint8,\n        rasterio.uint16, or rasterio.float32\n    size : int\n        minimum polygon size (number of pixels) to retain.\n    out : numpy ndarray, optional\n        Array of same shape and data type as `image` in which to store results.\n    mask : numpy ndarray or rasterio Band object, optional\n        Values of False or 0 will be excluded from feature generation\n        Must evaluate to bool (rasterio.bool_ or rasterio.uint8)\n    connectivity : int, optional\n        Use 4 or 8 pixel connectivity for grouping pixels into features\n\n    Returns\n    -------\n    out : numpy ndarray\n        Result\n\n    Notes\n    -----\n    GDAL only supports values that can be cast to 32-bit integers for this\n    operation.\n\n    The amount of memory used by this algorithm is proportional to the number\n    and complexity of polygons found in the image.  This algorithm is most\n    appropriate for simple thematic data.  Data with high pixel-to-pixel\n    variability, such as imagery, may produce one polygon per pixel and consume\n    large amounts of memory.\n\n    \"\"\"\n\n    if out is None:\n        out = np.zeros(image.shape, image.dtype)\n    _sieve(image, size, out, mask, connectivity)\n    return out\n\n\n@ensure_env\ndef rasterize(\n        shapes,\n        out_shape=None,\n        fill=0,\n        out=None,\n        transform=IDENTITY,\n        all_touched=False,\n        default_value=1,\n        dtype=None):\n    \"\"\"Return an image array with input geometries burned in.\n\n    Parameters\n    ----------\n    shapes : iterable of (geometry, value) pairs or iterable over\n        geometries. `geometry` can either be an object that implements\n        the geo interface or GeoJSON-like object.\n    out_shape : tuple or list\n        Shape of output numpy ndarray.\n    fill : int or float, optional\n        Used as fill value for all areas not covered by input\n        geometries.\n    out : numpy ndarray, optional\n        Array of same shape and data type as `image` in which to store\n        results.\n    transform : Affine transformation object, optional\n        Transformation from pixel coordinates of `image` to the\n        coordinate system of the input `shapes`. See the `transform`\n        property of dataset objects.\n    all_touched : boolean, optional\n        If True, all pixels touched by geometries will be burned in.  If\n        false, only pixels whose center is within the polygon or that\n        are selected by Bresenham's line algorithm will be burned in.\n    default_value : int or float, optional\n        Used as value for all geometries, if not provided in `shapes`.\n    dtype : rasterio or numpy data type, optional\n        Used as data type for results, if `out` is not provided.\n\n    Returns\n    -------\n    out : numpy ndarray\n        Results\n\n    Notes\n    -----\n    Valid data types for `fill`, `default_value`, `out`, `dtype` and\n    shape values are rasterio.int16, rasterio.int32, rasterio.uint8,\n    rasterio.uint16, rasterio.uint32, rasterio.float32,\n    rasterio.float64.\n\n    \"\"\"\n    valid_dtypes = (\n        'int16', 'int32', 'uint8', 'uint16', 'uint32', 'float32', 'float64'\n    )\n\n    def format_invalid_dtype(param):\n        return '{0} dtype must be one of: {1}'.format(\n            param, ', '.join(valid_dtypes)\n        )\n\n    def format_cast_error(param, dtype):\n        return '{0} cannot be cast to specified dtype: {1}'.format(param, dtype)\n\n    if fill != 0:\n        fill_array = np.array([fill])\n        if not validate_dtype(fill_array, valid_dtypes):\n            raise ValueError(format_invalid_dtype('fill'))\n\n        if dtype is not None and not can_cast_dtype(fill_array, dtype):\n            raise ValueError(format_cast_error('fill', dtype))\n\n    if default_value != 1:\n        default_value_array = np.array([default_value])\n        if not validate_dtype(default_value_array, valid_dtypes):\n            raise ValueError(format_invalid_dtype('default_value'))\n\n        if dtype is not None and not can_cast_dtype(default_value_array, dtype):\n            raise ValueError(format_cast_error('default_vaue', dtype))\n\n    if dtype is not None and np.dtype(dtype).name not in valid_dtypes:\n        raise ValueError(format_invalid_dtype('dtype'))\n\n    valid_shapes = []\n    shape_values = []\n    for index, item in enumerate(shapes):\n        if isinstance(item, (tuple, list)):\n            geom, value = item\n        else:\n            geom = item\n            value = default_value\n        geom = getattr(geom, '__geo_interface__', None) or geom\n\n        # not isinstance(geom, dict) or\n        if 'type' in geom or 'coordinates' in geom:\n            valid_shapes.append((geom, value))\n            shape_values.append(value)\n\n        else:\n            raise ValueError(\n                'Invalid geometry object at index {0}'.format(index)\n            )\n\n    if not valid_shapes:\n        raise ValueError('No valid geometry objects found for rasterize')\n\n    shape_values = np.array(shape_values)\n\n    if not validate_dtype(shape_values, valid_dtypes):\n        raise ValueError(format_invalid_dtype('shape values'))\n\n    if dtype is None:\n        dtype = get_minimum_dtype(np.append(shape_values, fill))\n\n    elif not can_cast_dtype(shape_values, dtype):\n        raise ValueError(format_cast_error('shape values', dtype))\n\n    if out is not None:\n        if np.dtype(out.dtype).name not in valid_dtypes:\n            raise ValueError(format_invalid_dtype('out'))\n\n        if not can_cast_dtype(shape_values, out.dtype):\n            raise ValueError(format_cast_error('shape values', out.dtype.name))\n\n    elif out_shape is not None:\n        out = np.empty(out_shape, dtype=dtype)\n        out.fill(fill)\n\n    else:\n        raise ValueError('Either an out_shape or image must be provided')\n\n    transform = guard_transform(transform)\n    _rasterize(valid_shapes, out, transform.to_gdal(), all_touched)\n    return out\n\n\ndef bounds(geometry):\n    \"\"\"Return a (minx, miny, maxx, maxy) bounding box.\n\n    From Fiona 1.4.8. Modified to return bbox from geometry if available.\n\n    Parameters\n    ----------\n    geometry: GeoJSON-like feature, feature collection, or geometry.\n\n    Returns\n    -------\n    tuple\n        Bounding box: (minx, miny, maxx, maxy)\n    \"\"\"\n    if 'bbox' in geometry:\n        return tuple(geometry['bbox'])\n\n    geom = geometry.get('geometry') or geometry\n    return _bounds(geom)\n", "path": "rasterio/features.py"}], "after_files": [{"content": "\"\"\"Functions for working with features in a raster dataset.\"\"\"\n\n\nimport logging\n\nimport numpy as np\n\nfrom rasterio._features import _shapes, _sieve, _rasterize, _bounds\nfrom rasterio.dtypes import validate_dtype, can_cast_dtype, get_minimum_dtype\nfrom rasterio.env import ensure_env\nfrom rasterio.transform import IDENTITY, guard_transform\n\n\nlog = logging.getLogger(__name__)\n\n\n@ensure_env\ndef geometry_mask(\n        geometries,\n        out_shape,\n        transform,\n        all_touched=False,\n        invert=False):\n    \"\"\"Create a mask from shapes.\n\n    By default, mask is intended for use as a\n    numpy mask, where pixels that overlap shapes are False.\n\n    Parameters\n    ----------\n    geometries : iterable over geometries (GeoJSON-like objects)\n    out_shape : tuple or list\n        Shape of output numpy ndarray.\n    transform : Affine transformation object\n        Transformation from pixel coordinates of `image` to the\n        coordinate system of the input `shapes`. See the `transform`\n        property of dataset objects.\n    all_touched : boolean, optional\n        If True, all pixels touched by geometries will be burned in.  If\n        false, only pixels whose center is within the polygon or that\n        are selected by Bresenham's line algorithm will be burned in.\n    invert: boolean, optional\n        If True, mask will be True for pixels that overlap shapes.\n        False by default.\n\n    Returns\n    -------\n    out : numpy ndarray of type 'bool'\n        Result\n    \"\"\"\n    fill, mask_value = (0, 1) if invert else (1, 0)\n\n    return rasterize(\n        geometries,\n        out_shape=out_shape,\n        transform=transform,\n        all_touched=all_touched,\n        fill=fill,\n        default_value=mask_value).astype('bool')\n\n\n@ensure_env\ndef shapes(image, mask=None, connectivity=4, transform=IDENTITY):\n    \"\"\"Yield (polygon, value for each set of adjacent pixels of the same value.\n\n    Parameters\n    ----------\n    image : numpy ndarray or rasterio Band object\n        (RasterReader, bidx namedtuple).\n        Data type must be one of rasterio.int16, rasterio.int32,\n        rasterio.uint8, rasterio.uint16, or rasterio.float32.\n    mask : numpy ndarray or rasterio Band object, optional\n        Values of False or 0 will be excluded from feature generation\n        Must evaluate to bool (rasterio.bool_ or rasterio.uint8)\n    connectivity : int, optional\n        Use 4 or 8 pixel connectivity for grouping pixels into features\n    transform : Affine transformation, optional\n        If not provided, feature coordinates will be generated based on pixel\n        coordinates\n\n    Yields\n    -------\n    tuple\n        A pair of (polygon, value) for each feature found in the image.\n        Polygons are GeoJSON-like dicts and the values are the associated value\n        from the image, in the data type of the image.\n        Note: due to floating point precision issues, values returned from a\n        floating point image may not exactly match the original values.\n\n    Notes\n    -----\n    The amount of memory used by this algorithm is proportional to the number\n    and complexity of polygons produced.  This algorithm is most appropriate\n    for simple thematic data.  Data with high pixel-to-pixel variability, such\n    as imagery, may produce one polygon per pixel and consume large amounts of\n    memory.\n\n    \"\"\"\n    transform = guard_transform(transform)\n    for s, v in _shapes(image, mask, connectivity, transform.to_gdal()):\n        yield s, v\n\n\n@ensure_env\ndef sieve(image, size, out=None, mask=None, connectivity=4):\n    \"\"\"Replace small polygons in `image` with value of their largest neighbor.\n\n    Polygons are found for each set of neighboring pixels of the same value.\n\n    Parameters\n    ----------\n    image : numpy ndarray or rasterio Band object\n        (RasterReader, bidx namedtuple)\n        Must be of type rasterio.int16, rasterio.int32, rasterio.uint8,\n        rasterio.uint16, or rasterio.float32\n    size : int\n        minimum polygon size (number of pixels) to retain.\n    out : numpy ndarray, optional\n        Array of same shape and data type as `image` in which to store results.\n    mask : numpy ndarray or rasterio Band object, optional\n        Values of False or 0 will be excluded from feature generation\n        Must evaluate to bool (rasterio.bool_ or rasterio.uint8)\n    connectivity : int, optional\n        Use 4 or 8 pixel connectivity for grouping pixels into features\n\n    Returns\n    -------\n    out : numpy ndarray\n        Result\n\n    Notes\n    -----\n    GDAL only supports values that can be cast to 32-bit integers for this\n    operation.\n\n    The amount of memory used by this algorithm is proportional to the number\n    and complexity of polygons found in the image.  This algorithm is most\n    appropriate for simple thematic data.  Data with high pixel-to-pixel\n    variability, such as imagery, may produce one polygon per pixel and consume\n    large amounts of memory.\n\n    \"\"\"\n\n    if out is None:\n        out = np.zeros(image.shape, image.dtype)\n    _sieve(image, size, out, mask, connectivity)\n    return out\n\n\n@ensure_env\ndef rasterize(\n        shapes,\n        out_shape=None,\n        fill=0,\n        out=None,\n        transform=IDENTITY,\n        all_touched=False,\n        default_value=1,\n        dtype=None):\n    \"\"\"Return an image array with input geometries burned in.\n\n    Parameters\n    ----------\n    shapes : iterable of (geometry, value) pairs or iterable over\n        geometries. `geometry` can either be an object that implements\n        the geo interface or GeoJSON-like object.\n    out_shape : tuple or list with 2 integers\n        Shape of output numpy ndarray.\n    fill : int or float, optional\n        Used as fill value for all areas not covered by input\n        geometries.\n    out : numpy ndarray, optional\n        Array of same shape and data type as `image` in which to store\n        results.\n    transform : Affine transformation object, optional\n        Transformation from pixel coordinates of `image` to the\n        coordinate system of the input `shapes`. See the `transform`\n        property of dataset objects.\n    all_touched : boolean, optional\n        If True, all pixels touched by geometries will be burned in.  If\n        false, only pixels whose center is within the polygon or that\n        are selected by Bresenham's line algorithm will be burned in.\n    default_value : int or float, optional\n        Used as value for all geometries, if not provided in `shapes`.\n    dtype : rasterio or numpy data type, optional\n        Used as data type for results, if `out` is not provided.\n\n    Returns\n    -------\n    out : numpy ndarray\n        Results\n\n    Notes\n    -----\n    Valid data types for `fill`, `default_value`, `out`, `dtype` and\n    shape values are rasterio.int16, rasterio.int32, rasterio.uint8,\n    rasterio.uint16, rasterio.uint32, rasterio.float32,\n    rasterio.float64.\n\n    \"\"\"\n    valid_dtypes = (\n        'int16', 'int32', 'uint8', 'uint16', 'uint32', 'float32', 'float64'\n    )\n\n    def format_invalid_dtype(param):\n        return '{0} dtype must be one of: {1}'.format(\n            param, ', '.join(valid_dtypes)\n        )\n\n    def format_cast_error(param, dtype):\n        return '{0} cannot be cast to specified dtype: {1}'.format(param, dtype)\n\n    if fill != 0:\n        fill_array = np.array([fill])\n        if not validate_dtype(fill_array, valid_dtypes):\n            raise ValueError(format_invalid_dtype('fill'))\n\n        if dtype is not None and not can_cast_dtype(fill_array, dtype):\n            raise ValueError(format_cast_error('fill', dtype))\n\n    if default_value != 1:\n        default_value_array = np.array([default_value])\n        if not validate_dtype(default_value_array, valid_dtypes):\n            raise ValueError(format_invalid_dtype('default_value'))\n\n        if dtype is not None and not can_cast_dtype(default_value_array, dtype):\n            raise ValueError(format_cast_error('default_vaue', dtype))\n\n    if dtype is not None and np.dtype(dtype).name not in valid_dtypes:\n        raise ValueError(format_invalid_dtype('dtype'))\n\n    valid_shapes = []\n    shape_values = []\n    for index, item in enumerate(shapes):\n        if isinstance(item, (tuple, list)):\n            geom, value = item\n        else:\n            geom = item\n            value = default_value\n        geom = getattr(geom, '__geo_interface__', None) or geom\n\n        # not isinstance(geom, dict) or\n        if 'type' in geom or 'coordinates' in geom:\n            valid_shapes.append((geom, value))\n            shape_values.append(value)\n\n        else:\n            raise ValueError(\n                'Invalid geometry object at index {0}'.format(index)\n            )\n\n    if not valid_shapes:\n        raise ValueError('No valid geometry objects found for rasterize')\n\n    shape_values = np.array(shape_values)\n\n    if not validate_dtype(shape_values, valid_dtypes):\n        raise ValueError(format_invalid_dtype('shape values'))\n\n    if dtype is None:\n        dtype = get_minimum_dtype(np.append(shape_values, fill))\n\n    elif not can_cast_dtype(shape_values, dtype):\n        raise ValueError(format_cast_error('shape values', dtype))\n\n    if out is not None:\n        if np.dtype(out.dtype).name not in valid_dtypes:\n            raise ValueError(format_invalid_dtype('out'))\n\n        if not can_cast_dtype(shape_values, out.dtype):\n            raise ValueError(format_cast_error('shape values', out.dtype.name))\n\n    elif out_shape is not None:\n\n        if len(out_shape) != 2:\n            raise ValueError('Invalid out_shape, must be 2D')\n\n        out = np.empty(out_shape, dtype=dtype)\n        out.fill(fill)\n\n    else:\n        raise ValueError('Either an out_shape or image must be provided')\n\n    transform = guard_transform(transform)\n    _rasterize(valid_shapes, out, transform.to_gdal(), all_touched)\n    return out\n\n\ndef bounds(geometry):\n    \"\"\"Return a (minx, miny, maxx, maxy) bounding box.\n\n    From Fiona 1.4.8. Modified to return bbox from geometry if available.\n\n    Parameters\n    ----------\n    geometry: GeoJSON-like feature, feature collection, or geometry.\n\n    Returns\n    -------\n    tuple\n        Bounding box: (minx, miny, maxx, maxy)\n    \"\"\"\n    if 'bbox' in geometry:\n        return tuple(geometry['bbox'])\n\n    geom = geometry.get('geometry') or geometry\n    return _bounds(geom)\n", "path": "rasterio/features.py"}]}