We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
Cast nodata value before filling in rio calc
At line 192, `res = res.filled(kwargs['nodata'])`, we need to convert that string to a numpy scalar or array of `res` dtype.

--- 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/rio/calc.py`
Content:
```
1 """$ rio calc"""
2 
3 from __future__ import division
4 
5 from collections import OrderedDict
6 from distutils.version import LooseVersion
7 import math
8 
9 import click
10 import snuggs
11 
12 import rasterio
13 from rasterio.features import sieve
14 from rasterio.fill import fillnodata
15 from rasterio.windows import Window
16 from rasterio.rio import options
17 from rasterio.rio.helpers import resolve_inout
18 
19 
20 def _get_bands(inputs, sources, d, i=None):
21     """Get a rasterio.Band object from calc's inputs"""
22     idx = d if d in dict(inputs) else int(d) - 1
23     src = sources[idx]
24     return (rasterio.band(src, i) if i else
25             [rasterio.band(src, j) for j in src.indexes])
26 
27 
28 def _read_array(ix, subix=None, dtype=None):
29     """Change the type of a read array"""
30     arr = snuggs._ctx.lookup(ix, subix)
31     if dtype:
32         arr = arr.astype(dtype)
33     return arr
34 
35 
36 def _chunk_output(width, height, count, itemsize, mem_limit=1):
37     """Divide the calculation output into chunks
38 
39     This function determines the chunk size such that an array of shape
40     (chunk_size, chunk_size, count) with itemsize bytes per element
41     requires no more than mem_limit megabytes of memory.
42 
43     Output chunks are described by rasterio Windows.
44 
45     Parameters
46     ----------
47     width : int
48         Output width
49     height : int
50         Output height
51     count : int
52         Number of output bands
53     itemsize : int
54         Number of bytes per pixel
55     mem_limit : int, default
56         The maximum size in memory of a chunk array
57 
58     Returns
59     -------
60     sequence of Windows
61     """
62     max_pixels = mem_limit * 1.0e+6 / itemsize * count
63     chunk_size = int(math.floor(math.sqrt(max_pixels)))
64     ncols = int(math.ceil(width / chunk_size))
65     nrows = int(math.ceil(height / chunk_size))
66     chunk_windows = []
67 
68     for col in range(ncols):
69         col_offset = col * chunk_size
70         w = min(chunk_size, width - col_offset)
71         for row in range(nrows):
72             row_offset = row * chunk_size
73             h = min(chunk_size, height - row_offset)
74             chunk_windows.append(((row, col), Window(col_offset, row_offset, w, h)))
75 
76     return chunk_windows
77 
78 
79 @click.command(short_help="Raster data calculator.")
80 @click.argument('command')
81 @options.files_inout_arg
82 @options.output_opt
83 @options.format_opt
84 @click.option('--name', multiple=True,
85               help='Specify an input file with a unique short (alphas only) '
86                    'name for use in commands like '
87                    '"a=tests/data/RGB.byte.tif".')
88 @options.dtype_opt
89 @options.masked_opt
90 @options.overwrite_opt
91 @click.option("--mem-limit", type=int, default=64, help="Limit on memory used to perform calculations, in MB.")
92 @options.creation_options
93 @click.pass_context
94 def calc(ctx, command, files, output, driver, name, dtype, masked, overwrite, mem_limit, creation_options):
95     """A raster data calculator
96 
97     Evaluates an expression using input datasets and writes the result
98     to a new dataset.
99 
100     Command syntax is lisp-like. An expression consists of an operator
101     or function name and one or more strings, numbers, or expressions
102     enclosed in parentheses. Functions include ``read`` (gets a raster
103     array) and ``asarray`` (makes a 3-D array from 2-D arrays).
104 
105     \b
106         * (read i) evaluates to the i-th input dataset (a 3-D array).
107         * (read i j) evaluates to the j-th band of the i-th dataset (a
108           2-D array).
109         * (take foo j) evaluates to the j-th band of a dataset named foo
110           (see help on the --name option above).
111         * Standard numpy array operators (+, -, *, /) are available.
112         * When the final result is a list of arrays, a multiple band
113           output file is written.
114         * When the final result is a single array, a single band output
115           file is written.
116 
117     Example:
118 
119     \b
120          $ rio calc "(+ 2 (* 0.95 (read 1)))" tests/data/RGB.byte.tif \\
121          > /tmp/out.tif
122 
123     The command above produces a 3-band GeoTIFF with all values scaled
124     by 0.95 and incremented by 2.
125 
126     \b
127         $ rio calc "(asarray (+ 125 (read 1)) (read 1) (read 1))" \\
128         > tests/data/shade.tif /tmp/out.tif
129 
130     The command above produces a 3-band RGB GeoTIFF, with red levels
131     incremented by 125, from the single-band input.
132 
133     The maximum amount of memory used to perform caculations defaults to
134     64 MB. This number can be increased to improve speed of calculation.
135 
136     """
137     import numpy as np
138 
139     dst = None
140     sources = []
141 
142     try:
143         with ctx.obj['env']:
144             output, files = resolve_inout(files=files, output=output,
145                                           overwrite=overwrite)
146             inputs = ([tuple(n.split('=')) for n in name] +
147                       [(None, n) for n in files])
148             sources = [rasterio.open(path) for name, path in inputs]
149 
150             first = sources[0]
151             kwargs = first.profile
152             kwargs.update(**creation_options)
153             dtype = dtype or first.meta['dtype']
154             kwargs['dtype'] = dtype
155             kwargs.pop("driver", None)
156             if driver:
157                 kwargs['driver'] = driver
158 
159             # Extend snuggs.
160             snuggs.func_map['read'] = _read_array
161             snuggs.func_map['band'] = lambda d, i: _get_bands(inputs, sources, d, i)
162             snuggs.func_map['bands'] = lambda d: _get_bands(inputs, sources, d)
163             snuggs.func_map['fillnodata'] = lambda *args: fillnodata(*args)
164             snuggs.func_map['sieve'] = lambda *args: sieve(*args)
165 
166             # The windows iterator is initialized with a single sample.
167             # The actual work windows will be added in the second
168             # iteration of the loop.
169             work_windows = [(None, Window(0, 0, 16, 16))]
170 
171             for ij, window in work_windows:
172 
173                 ctxkwds = OrderedDict()
174 
175                 for i, ((name, path), src) in enumerate(zip(inputs, sources)):
176 
177                     # Using the class method instead of instance
178                     # method. Latter raises
179                     #
180                     # TypeError: astype() got an unexpected keyword
181                     # argument 'copy'
182                     #
183                     # possibly something to do with the instance being
184                     # a masked array.
185                     ctxkwds[name or '_i%d' % (i + 1)] = src.read(masked=masked, window=window)
186 
187                 res = snuggs.eval(command, **ctxkwds)
188 
189                 if (isinstance(res, np.ma.core.MaskedArray) and (
190                         tuple(LooseVersion(np.__version__).version) < (1, 9) or
191                         tuple(LooseVersion(np.__version__).version) > (1, 10))):
192                     res = res.filled(kwargs['nodata'])
193 
194                 if len(res.shape) == 3:
195                     results = np.ndarray.astype(res, dtype, copy=False)
196                 else:
197                     results = np.asanyarray(
198                         [np.ndarray.astype(res, dtype, copy=False)])
199 
200                 # The first iteration is only to get sample results and from them
201                 # compute some properties of the output dataset.
202                 if dst is None:
203                     kwargs['count'] = results.shape[0]
204                     dst = rasterio.open(output, 'w', **kwargs)
205                     work_windows.extend(_chunk_output(dst.width, dst.height, dst.count, np.dtype(dst.dtypes[0]).itemsize, mem_limit=mem_limit))
206 
207                 # In subsequent iterations we write results.
208                 else:
209                     dst.write(results, window=window)
210 
211     except snuggs.ExpressionError as err:
212         click.echo("Expression Error:")
213         click.echo("  {}".format(err.text))
214         click.echo(" {}^".format(" " * err.offset))
215         click.echo(err)
216         raise click.Abort()
217 
218     finally:
219         if dst:
220             dst.close()
221         for src in sources:
222             src.close()
223
```
--- 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/rio/calc.py b/rasterio/rio/calc.py
--- a/rasterio/rio/calc.py
+++ b/rasterio/rio/calc.py
@@ -185,17 +185,14 @@
                     ctxkwds[name or '_i%d' % (i + 1)] = src.read(masked=masked, window=window)
 
                 res = snuggs.eval(command, **ctxkwds)
-
-                if (isinstance(res, np.ma.core.MaskedArray) and (
-                        tuple(LooseVersion(np.__version__).version) < (1, 9) or
-                        tuple(LooseVersion(np.__version__).version) > (1, 10))):
-                    res = res.filled(kwargs['nodata'])
-
-                if len(res.shape) == 3:
-                    results = np.ndarray.astype(res, dtype, copy=False)
-                else:
-                    results = np.asanyarray(
-                        [np.ndarray.astype(res, dtype, copy=False)])
+                results = res.astype(dtype, copy=False)
+
+                if isinstance(results, np.ma.core.MaskedArray):
+                    results = results.filled(float(kwargs['nodata']))
+                    if len(results.shape) == 2:
+                        results = np.ma.asanyarray([results])
+                elif len(results.shape) == 2:
+                    results = np.asanyarray([results])
 
                 # The first iteration is only to get sample results and from them
                 # compute some properties of the output dataset.

{"golden_diff": "diff --git a/rasterio/rio/calc.py b/rasterio/rio/calc.py\n--- a/rasterio/rio/calc.py\n+++ b/rasterio/rio/calc.py\n@@ -185,17 +185,14 @@\n                     ctxkwds[name or '_i%d' % (i + 1)] = src.read(masked=masked, window=window)\n \n                 res = snuggs.eval(command, **ctxkwds)\n-\n-                if (isinstance(res, np.ma.core.MaskedArray) and (\n-                        tuple(LooseVersion(np.__version__).version) < (1, 9) or\n-                        tuple(LooseVersion(np.__version__).version) > (1, 10))):\n-                    res = res.filled(kwargs['nodata'])\n-\n-                if len(res.shape) == 3:\n-                    results = np.ndarray.astype(res, dtype, copy=False)\n-                else:\n-                    results = np.asanyarray(\n-                        [np.ndarray.astype(res, dtype, copy=False)])\n+                results = res.astype(dtype, copy=False)\n+\n+                if isinstance(results, np.ma.core.MaskedArray):\n+                    results = results.filled(float(kwargs['nodata']))\n+                    if len(results.shape) == 2:\n+                        results = np.ma.asanyarray([results])\n+                elif len(results.shape) == 2:\n+                    results = np.asanyarray([results])\n \n                 # The first iteration is only to get sample results and from them\n                 # compute some properties of the output dataset.\n", "issue": "Cast nodata value before filling in rio calc\nAt line 192, `res = res.filled(kwargs['nodata'])`, we need to convert that string to a numpy scalar or array of `res` dtype.\n", "before_files": [{"content": "\"\"\"$ rio calc\"\"\"\n\nfrom __future__ import division\n\nfrom collections import OrderedDict\nfrom distutils.version import LooseVersion\nimport math\n\nimport click\nimport snuggs\n\nimport rasterio\nfrom rasterio.features import sieve\nfrom rasterio.fill import fillnodata\nfrom rasterio.windows import Window\nfrom rasterio.rio import options\nfrom rasterio.rio.helpers import resolve_inout\n\n\ndef _get_bands(inputs, sources, d, i=None):\n    \"\"\"Get a rasterio.Band object from calc's inputs\"\"\"\n    idx = d if d in dict(inputs) else int(d) - 1\n    src = sources[idx]\n    return (rasterio.band(src, i) if i else\n            [rasterio.band(src, j) for j in src.indexes])\n\n\ndef _read_array(ix, subix=None, dtype=None):\n    \"\"\"Change the type of a read array\"\"\"\n    arr = snuggs._ctx.lookup(ix, subix)\n    if dtype:\n        arr = arr.astype(dtype)\n    return arr\n\n\ndef _chunk_output(width, height, count, itemsize, mem_limit=1):\n    \"\"\"Divide the calculation output into chunks\n\n    This function determines the chunk size such that an array of shape\n    (chunk_size, chunk_size, count) with itemsize bytes per element\n    requires no more than mem_limit megabytes of memory.\n\n    Output chunks are described by rasterio Windows.\n\n    Parameters\n    ----------\n    width : int\n        Output width\n    height : int\n        Output height\n    count : int\n        Number of output bands\n    itemsize : int\n        Number of bytes per pixel\n    mem_limit : int, default\n        The maximum size in memory of a chunk array\n\n    Returns\n    -------\n    sequence of Windows\n    \"\"\"\n    max_pixels = mem_limit * 1.0e+6 / itemsize * count\n    chunk_size = int(math.floor(math.sqrt(max_pixels)))\n    ncols = int(math.ceil(width / chunk_size))\n    nrows = int(math.ceil(height / chunk_size))\n    chunk_windows = []\n\n    for col in range(ncols):\n        col_offset = col * chunk_size\n        w = min(chunk_size, width - col_offset)\n        for row in range(nrows):\n            row_offset = row * chunk_size\n            h = min(chunk_size, height - row_offset)\n            chunk_windows.append(((row, col), Window(col_offset, row_offset, w, h)))\n\n    return chunk_windows\n\n\n@click.command(short_help=\"Raster data calculator.\")\n@click.argument('command')\n@options.files_inout_arg\n@options.output_opt\n@options.format_opt\n@click.option('--name', multiple=True,\n              help='Specify an input file with a unique short (alphas only) '\n                   'name for use in commands like '\n                   '\"a=tests/data/RGB.byte.tif\".')\n@options.dtype_opt\n@options.masked_opt\n@options.overwrite_opt\n@click.option(\"--mem-limit\", type=int, default=64, help=\"Limit on memory used to perform calculations, in MB.\")\n@options.creation_options\n@click.pass_context\ndef calc(ctx, command, files, output, driver, name, dtype, masked, overwrite, mem_limit, creation_options):\n    \"\"\"A raster data calculator\n\n    Evaluates an expression using input datasets and writes the result\n    to a new dataset.\n\n    Command syntax is lisp-like. An expression consists of an operator\n    or function name and one or more strings, numbers, or expressions\n    enclosed in parentheses. Functions include ``read`` (gets a raster\n    array) and ``asarray`` (makes a 3-D array from 2-D arrays).\n\n    \\b\n        * (read i) evaluates to the i-th input dataset (a 3-D array).\n        * (read i j) evaluates to the j-th band of the i-th dataset (a\n          2-D array).\n        * (take foo j) evaluates to the j-th band of a dataset named foo\n          (see help on the --name option above).\n        * Standard numpy array operators (+, -, *, /) are available.\n        * When the final result is a list of arrays, a multiple band\n          output file is written.\n        * When the final result is a single array, a single band output\n          file is written.\n\n    Example:\n\n    \\b\n         $ rio calc \"(+ 2 (* 0.95 (read 1)))\" tests/data/RGB.byte.tif \\\\\n         > /tmp/out.tif\n\n    The command above produces a 3-band GeoTIFF with all values scaled\n    by 0.95 and incremented by 2.\n\n    \\b\n        $ rio calc \"(asarray (+ 125 (read 1)) (read 1) (read 1))\" \\\\\n        > tests/data/shade.tif /tmp/out.tif\n\n    The command above produces a 3-band RGB GeoTIFF, with red levels\n    incremented by 125, from the single-band input.\n\n    The maximum amount of memory used to perform caculations defaults to\n    64 MB. This number can be increased to improve speed of calculation.\n\n    \"\"\"\n    import numpy as np\n\n    dst = None\n    sources = []\n\n    try:\n        with ctx.obj['env']:\n            output, files = resolve_inout(files=files, output=output,\n                                          overwrite=overwrite)\n            inputs = ([tuple(n.split('=')) for n in name] +\n                      [(None, n) for n in files])\n            sources = [rasterio.open(path) for name, path in inputs]\n\n            first = sources[0]\n            kwargs = first.profile\n            kwargs.update(**creation_options)\n            dtype = dtype or first.meta['dtype']\n            kwargs['dtype'] = dtype\n            kwargs.pop(\"driver\", None)\n            if driver:\n                kwargs['driver'] = driver\n\n            # Extend snuggs.\n            snuggs.func_map['read'] = _read_array\n            snuggs.func_map['band'] = lambda d, i: _get_bands(inputs, sources, d, i)\n            snuggs.func_map['bands'] = lambda d: _get_bands(inputs, sources, d)\n            snuggs.func_map['fillnodata'] = lambda *args: fillnodata(*args)\n            snuggs.func_map['sieve'] = lambda *args: sieve(*args)\n\n            # The windows iterator is initialized with a single sample.\n            # The actual work windows will be added in the second\n            # iteration of the loop.\n            work_windows = [(None, Window(0, 0, 16, 16))]\n\n            for ij, window in work_windows:\n\n                ctxkwds = OrderedDict()\n\n                for i, ((name, path), src) in enumerate(zip(inputs, sources)):\n\n                    # Using the class method instead of instance\n                    # method. Latter raises\n                    #\n                    # TypeError: astype() got an unexpected keyword\n                    # argument 'copy'\n                    #\n                    # possibly something to do with the instance being\n                    # a masked array.\n                    ctxkwds[name or '_i%d' % (i + 1)] = src.read(masked=masked, window=window)\n\n                res = snuggs.eval(command, **ctxkwds)\n\n                if (isinstance(res, np.ma.core.MaskedArray) and (\n                        tuple(LooseVersion(np.__version__).version) < (1, 9) or\n                        tuple(LooseVersion(np.__version__).version) > (1, 10))):\n                    res = res.filled(kwargs['nodata'])\n\n                if len(res.shape) == 3:\n                    results = np.ndarray.astype(res, dtype, copy=False)\n                else:\n                    results = np.asanyarray(\n                        [np.ndarray.astype(res, dtype, copy=False)])\n\n                # The first iteration is only to get sample results and from them\n                # compute some properties of the output dataset.\n                if dst is None:\n                    kwargs['count'] = results.shape[0]\n                    dst = rasterio.open(output, 'w', **kwargs)\n                    work_windows.extend(_chunk_output(dst.width, dst.height, dst.count, np.dtype(dst.dtypes[0]).itemsize, mem_limit=mem_limit))\n\n                # In subsequent iterations we write results.\n                else:\n                    dst.write(results, window=window)\n\n    except snuggs.ExpressionError as err:\n        click.echo(\"Expression Error:\")\n        click.echo(\"  {}\".format(err.text))\n        click.echo(\" {}^\".format(\" \" * err.offset))\n        click.echo(err)\n        raise click.Abort()\n\n    finally:\n        if dst:\n            dst.close()\n        for src in sources:\n            src.close()\n", "path": "rasterio/rio/calc.py"}], "after_files": [{"content": "\"\"\"$ rio calc\"\"\"\n\nfrom __future__ import division\n\nfrom collections import OrderedDict\nfrom distutils.version import LooseVersion\nimport math\n\nimport click\nimport snuggs\n\nimport rasterio\nfrom rasterio.features import sieve\nfrom rasterio.fill import fillnodata\nfrom rasterio.windows import Window\nfrom rasterio.rio import options\nfrom rasterio.rio.helpers import resolve_inout\n\n\ndef _get_bands(inputs, sources, d, i=None):\n    \"\"\"Get a rasterio.Band object from calc's inputs\"\"\"\n    idx = d if d in dict(inputs) else int(d) - 1\n    src = sources[idx]\n    return (rasterio.band(src, i) if i else\n            [rasterio.band(src, j) for j in src.indexes])\n\n\ndef _read_array(ix, subix=None, dtype=None):\n    \"\"\"Change the type of a read array\"\"\"\n    arr = snuggs._ctx.lookup(ix, subix)\n    if dtype:\n        arr = arr.astype(dtype)\n    return arr\n\n\ndef _chunk_output(width, height, count, itemsize, mem_limit=1):\n    \"\"\"Divide the calculation output into chunks\n\n    This function determines the chunk size such that an array of shape\n    (chunk_size, chunk_size, count) with itemsize bytes per element\n    requires no more than mem_limit megabytes of memory.\n\n    Output chunks are described by rasterio Windows.\n\n    Parameters\n    ----------\n    width : int\n        Output width\n    height : int\n        Output height\n    count : int\n        Number of output bands\n    itemsize : int\n        Number of bytes per pixel\n    mem_limit : int, default\n        The maximum size in memory of a chunk array\n\n    Returns\n    -------\n    sequence of Windows\n    \"\"\"\n    max_pixels = mem_limit * 1.0e+6 / itemsize * count\n    chunk_size = int(math.floor(math.sqrt(max_pixels)))\n    ncols = int(math.ceil(width / chunk_size))\n    nrows = int(math.ceil(height / chunk_size))\n    chunk_windows = []\n\n    for col in range(ncols):\n        col_offset = col * chunk_size\n        w = min(chunk_size, width - col_offset)\n        for row in range(nrows):\n            row_offset = row * chunk_size\n            h = min(chunk_size, height - row_offset)\n            chunk_windows.append(((row, col), Window(col_offset, row_offset, w, h)))\n\n    return chunk_windows\n\n\n@click.command(short_help=\"Raster data calculator.\")\n@click.argument('command')\n@options.files_inout_arg\n@options.output_opt\n@options.format_opt\n@click.option('--name', multiple=True,\n              help='Specify an input file with a unique short (alphas only) '\n                   'name for use in commands like '\n                   '\"a=tests/data/RGB.byte.tif\".')\n@options.dtype_opt\n@options.masked_opt\n@options.overwrite_opt\n@click.option(\"--mem-limit\", type=int, default=64, help=\"Limit on memory used to perform calculations, in MB.\")\n@options.creation_options\n@click.pass_context\ndef calc(ctx, command, files, output, driver, name, dtype, masked, overwrite, mem_limit, creation_options):\n    \"\"\"A raster data calculator\n\n    Evaluates an expression using input datasets and writes the result\n    to a new dataset.\n\n    Command syntax is lisp-like. An expression consists of an operator\n    or function name and one or more strings, numbers, or expressions\n    enclosed in parentheses. Functions include ``read`` (gets a raster\n    array) and ``asarray`` (makes a 3-D array from 2-D arrays).\n\n    \\b\n        * (read i) evaluates to the i-th input dataset (a 3-D array).\n        * (read i j) evaluates to the j-th band of the i-th dataset (a\n          2-D array).\n        * (take foo j) evaluates to the j-th band of a dataset named foo\n          (see help on the --name option above).\n        * Standard numpy array operators (+, -, *, /) are available.\n        * When the final result is a list of arrays, a multiple band\n          output file is written.\n        * When the final result is a single array, a single band output\n          file is written.\n\n    Example:\n\n    \\b\n         $ rio calc \"(+ 2 (* 0.95 (read 1)))\" tests/data/RGB.byte.tif \\\\\n         > /tmp/out.tif\n\n    The command above produces a 3-band GeoTIFF with all values scaled\n    by 0.95 and incremented by 2.\n\n    \\b\n        $ rio calc \"(asarray (+ 125 (read 1)) (read 1) (read 1))\" \\\\\n        > tests/data/shade.tif /tmp/out.tif\n\n    The command above produces a 3-band RGB GeoTIFF, with red levels\n    incremented by 125, from the single-band input.\n\n    The maximum amount of memory used to perform caculations defaults to\n    64 MB. This number can be increased to improve speed of calculation.\n\n    \"\"\"\n    import numpy as np\n\n    dst = None\n    sources = []\n\n    try:\n        with ctx.obj['env']:\n            output, files = resolve_inout(files=files, output=output,\n                                          overwrite=overwrite)\n            inputs = ([tuple(n.split('=')) for n in name] +\n                      [(None, n) for n in files])\n            sources = [rasterio.open(path) for name, path in inputs]\n\n            first = sources[0]\n            kwargs = first.profile\n            kwargs.update(**creation_options)\n            dtype = dtype or first.meta['dtype']\n            kwargs['dtype'] = dtype\n            kwargs.pop(\"driver\", None)\n            if driver:\n                kwargs['driver'] = driver\n\n            # Extend snuggs.\n            snuggs.func_map['read'] = _read_array\n            snuggs.func_map['band'] = lambda d, i: _get_bands(inputs, sources, d, i)\n            snuggs.func_map['bands'] = lambda d: _get_bands(inputs, sources, d)\n            snuggs.func_map['fillnodata'] = lambda *args: fillnodata(*args)\n            snuggs.func_map['sieve'] = lambda *args: sieve(*args)\n\n            # The windows iterator is initialized with a single sample.\n            # The actual work windows will be added in the second\n            # iteration of the loop.\n            work_windows = [(None, Window(0, 0, 16, 16))]\n\n            for ij, window in work_windows:\n\n                ctxkwds = OrderedDict()\n\n                for i, ((name, path), src) in enumerate(zip(inputs, sources)):\n\n                    # Using the class method instead of instance\n                    # method. Latter raises\n                    #\n                    # TypeError: astype() got an unexpected keyword\n                    # argument 'copy'\n                    #\n                    # possibly something to do with the instance being\n                    # a masked array.\n                    ctxkwds[name or '_i%d' % (i + 1)] = src.read(masked=masked, window=window)\n\n                res = snuggs.eval(command, **ctxkwds)\n                results = res.astype(dtype, copy=False)\n\n                if isinstance(results, np.ma.core.MaskedArray):\n                    results = results.filled(float(kwargs['nodata']))\n                    if len(results.shape) == 2:\n                        results = np.ma.asanyarray([results])\n                elif len(results.shape) == 2:\n                    results = np.asanyarray([results])\n\n                # The first iteration is only to get sample results and from them\n                # compute some properties of the output dataset.\n                if dst is None:\n                    kwargs['count'] = results.shape[0]\n                    dst = rasterio.open(output, 'w', **kwargs)\n                    work_windows.extend(_chunk_output(dst.width, dst.height, dst.count, np.dtype(dst.dtypes[0]).itemsize, mem_limit=mem_limit))\n\n                # In subsequent iterations we write results.\n                else:\n                    dst.write(results, window=window)\n\n    except snuggs.ExpressionError as err:\n        click.echo(\"Expression Error:\")\n        click.echo(\"  {}\".format(err.text))\n        click.echo(\" {}^\".format(\" \" * err.offset))\n        click.echo(err)\n        raise click.Abort()\n\n    finally:\n        if dst:\n            dst.close()\n        for src in sources:\n            src.close()\n", "path": "rasterio/rio/calc.py"}]}