-
-
Notifications
You must be signed in to change notification settings - Fork 443
Expand file tree
/
Copy patharray.py
More file actions
394 lines (344 loc) · 11.5 KB
/
Copy patharray.py
File metadata and controls
394 lines (344 loc) · 11.5 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
"""
Functions to use GRASS 2D and 3D rasters with NumPy.
Usage:
>>> import grass.script as gs
>>> from grass.script import array as garray
>>>
>>> # We create a temporary region that is only valid in this python session
... gs.use_temp_region()
>>> gs.run_command("g.region", n=80, e=120, t=60, s=0, w=0, b=0, res=20, res3=20)
>>>
>>> # Lets create a raster map numpy array
... # based at the current region settings
... map2d_1 = garray.array()
>>>
>>> # Write some data
... for y in range(map2d_1.shape[0]):
... for x in range(map2d_1.shape[1]):
... map2d_1[y, x] = y + x
>>> # Lets have a look at the array
... print(map2d_1)
[[0. 1. 2. 3. 4. 5.]
[1. 2. 3. 4. 5. 6.]
[2. 3. 4. 5. 6. 7.]
[3. 4. 5. 6. 7. 8.]]
>>> # This will write the numpy array as GRASS raster map
... # with name map2d_1
... map2d_1.write(mapname="map2d_1", overwrite=True)
0
>>>
>>> # We create a new array from raster map2d_1 to modify it
... map2d_2 = garray.array(mapname="map2d_1")
>>> # Don't do map2d_2 = map2d_1 % 3
... # because: this will overwrite the internal temporary filename
... map2d_2 %= 3
>>> # Show the result
... print(map2d_2)
[[0. 1. 2. 0. 1. 2.]
[1. 2. 0. 1. 2. 0.]
[2. 0. 1. 2. 0. 1.]
[0. 1. 2. 0. 1. 2.]]
>>> # Write the result as new raster map with name map2d_2
... map2d_2.write(mapname="map2d_2", overwrite=True)
0
>>>
>>> # Here we create a 3D raster map numpy array
... # based in the current region settings
... map3d_1 = garray.array3d()
>>>
>>> # Write some data
... # Note: the 3D array has map[depth][row][column] order
... for z in range(map3d_1.shape[0]):
... for y in range(map3d_1.shape[1]):
... for x in range(map3d_1.shape[2]):
... map3d_1[z, y, x] = z + y + x
>>> # Lets have a look at the 3D array
... print(map3d_1)
[[[ 0. 1. 2. 3. 4. 5.]
[ 1. 2. 3. 4. 5. 6.]
[ 2. 3. 4. 5. 6. 7.]
[ 3. 4. 5. 6. 7. 8.]]
<BLANKLINE>
[[ 1. 2. 3. 4. 5. 6.]
[ 2. 3. 4. 5. 6. 7.]
[ 3. 4. 5. 6. 7. 8.]
[ 4. 5. 6. 7. 8. 9.]]
<BLANKLINE>
[[ 2. 3. 4. 5. 6. 7.]
[ 3. 4. 5. 6. 7. 8.]
[ 4. 5. 6. 7. 8. 9.]
[ 5. 6. 7. 8. 9. 10.]]]
>>> # This will write the numpy array as GRASS 3D raster map
... # with name map3d_1
... map3d_1.write(mapname="map3d_1", overwrite=True)
0
>>> # We create a new 3D array from 3D raster map3d_1 to modify it
... map3d_2 = garray.array3d(mapname="map3d_1")
>>> # Don't do map3d_2 = map3d_1 % 3
... # because: this will overwrite the internal temporary filename
... map3d_2 %= 3
>>> # Show the result
... print(map3d_2)
[[[0. 1. 2. 0. 1. 2.]
[1. 2. 0. 1. 2. 0.]
[2. 0. 1. 2. 0. 1.]
[0. 1. 2. 0. 1. 2.]]
<BLANKLINE>
[[1. 2. 0. 1. 2. 0.]
[2. 0. 1. 2. 0. 1.]
[0. 1. 2. 0. 1. 2.]
[1. 2. 0. 1. 2. 0.]]
<BLANKLINE>
[[2. 0. 1. 2. 0. 1.]
[0. 1. 2. 0. 1. 2.]
[1. 2. 0. 1. 2. 0.]
[2. 0. 1. 2. 0. 1.]]]
>>> # Write the result as new 3D raster map with name map3d_2
... map3d_2.write(mapname="map3d_2", overwrite=True)
0
SPDX-FileCopyrightText: 2010-2021 Glynn Clements
SPDX-FileCopyrightText: GRASS Development Team
SPDX-License-Identifier: GPL-2.0-or-later
.. sectionauthor:: Glynn Clements
"""
import numpy as np
from .utils import try_remove
from . import core as gcore
from grass.exceptions import CalledModuleError
###############################################################################
class _tempfile:
def __init__(self, env=None):
self.filename = gcore.tempfile(env=env)
def __del__(self):
try_remove(self.filename)
###############################################################################
class array(np.memmap):
# pylint: disable-next=signature-differs; W0222
def __new__(cls, mapname=None, null=None, dtype=None, env=None):
"""Define new numpy array
:param cls:
:param dtype: data type (based on map type, fallbacks to numpy.double)
:param env: environment
"""
reg = gcore.region(env=env)
r = reg["rows"]
c = reg["cols"]
shape = (r, c)
tempfile = _tempfile(env)
if mapname:
if not dtype:
try:
map_type = gcore.parse_command(
"r.info", map=mapname, format="json", env=env
)["datatype"]
if map_type == "CELL":
dtype = np.int32
elif map_type == "FCELL":
dtype = np.float32
elif map_type == "DCELL":
dtype = np.float64
except CalledModuleError:
dtype = np.double
kind = np.dtype(dtype).kind
size = np.dtype(dtype).itemsize
if kind == "f":
flags = "f"
elif kind in "biu":
if size == 8:
message = (
"64-bit integers are not supported by GRASS raster maps. "
"Use dtype=numpy.int32 or a smaller integer type."
)
raise ValueError(message)
flags = "i"
else:
raise ValueError(_("Invalid kind <%s>") % kind)
if size not in {1, 2, 4, 8}:
raise ValueError(_("Invalid size <%d>") % size)
gcore.run_command(
"r.out.bin",
flags=flags,
input=mapname,
output=tempfile.filename,
bytes=size,
null=null,
quiet=True,
overwrite=True,
env=env,
)
self = np.memmap.__new__(
cls,
filename=tempfile.filename,
dtype=dtype or np.double,
mode="r+",
shape=shape,
)
self.tempfile = tempfile
self.filename = tempfile.filename
self._env = env
return self
def write(self, mapname, title=None, null=None, overwrite=None, quiet=None):
"""Write array into raster map
:param str mapname: name for raster map
:param str title: title for raster map
:param null: null value
:param bool overwrite: True for overwriting existing raster maps
:return: 0 on success
:return: non-zero code on failure
"""
kind = self.dtype.kind
size = self.dtype.itemsize
if kind == "f":
if size == 4:
flags = "f"
elif size == 8:
flags = "d"
else:
raise ValueError(_("Invalid FP size <%d>") % size)
size = None
elif kind in "biu":
if size == 8:
message = (
"64-bit integers are not supported by GRASS raster maps. "
"Cast to a supported type before writing, e.g., "
"array.astype(numpy.int32)"
)
raise ValueError(message)
if size not in {1, 2, 4}:
raise ValueError(_("Invalid integer size <%d>") % size)
flags = None
else:
raise ValueError(_("Invalid kind <%s>") % kind)
# ensure all array content is written to the file
self.flush()
reg = gcore.region(env=self._env)
try:
gcore.run_command(
"r.in.bin",
flags=flags,
input=self.filename,
output=mapname,
title=title,
bytes=size,
anull=null,
overwrite=overwrite,
quiet=quiet,
north=reg["n"],
south=reg["s"],
east=reg["e"],
west=reg["w"],
rows=reg["rows"],
cols=reg["cols"],
env=self._env,
)
except CalledModuleError:
return 1
else:
return 0
###############################################################################
class array3d(np.memmap):
# pylint: disable-next=signature-differs; W0222
def __new__(cls, mapname=None, null=None, dtype=None, env=None):
"""Define new 3d numpy array
:param cls:
:param dtype: data type (based on map type, fallbacks to numpy.double)
:param env: environment
"""
reg = gcore.region(True, env=env)
r = reg["rows3"]
c = reg["cols3"]
d = reg["depths"]
shape = (d, r, c)
tempfile = _tempfile(env=env)
if mapname:
if not dtype:
try:
map_type = gcore.parse_command(
"r3.info", map=mapname, format="json", env=env
)["datatype"]
if map_type == "FCELL":
dtype = np.float32
elif map_type == "DCELL":
dtype = np.float64
except CalledModuleError:
dtype = np.double
kind = np.dtype(dtype).kind
size = np.dtype(dtype).itemsize
if kind == "f":
flags = None # default is double
elif kind in "biu":
flags = "i"
else:
raise ValueError(_("Invalid kind <%s>") % kind)
if size not in {1, 2, 4, 8}:
raise ValueError(_("Invalid size <%d>") % size)
gcore.run_command(
"r3.out.bin",
flags=flags,
input=mapname,
output=tempfile.filename,
bytes=size,
null=null,
quiet=True,
overwrite=True,
env=env,
)
self = np.memmap.__new__(
cls,
filename=tempfile.filename,
dtype=dtype or np.double,
mode="r+",
shape=shape,
)
self.tempfile = tempfile
self.filename = tempfile.filename
self._env = env
return self
def write(self, mapname, null=None, overwrite=None, quiet=None):
"""Write array into 3D raster map
:param str mapname: name for 3D raster map
:param null: null value
:param bool overwrite: True for overwriting existing raster maps
:return: 0 on success
:return: non-zero code on failure
"""
kind = self.dtype.kind
size = self.dtype.itemsize
flags = None
if kind == "f":
if size not in {4, 8}:
raise ValueError(_("Invalid FP size <%d>") % size)
elif kind in "biu":
if size not in {1, 2, 4, 8}:
raise ValueError(_("Invalid integer size <%d>") % size)
flags = "i"
else:
raise ValueError(_("Invalid kind <%s>") % kind)
# ensure all array content is written to the file
self.flush()
reg = gcore.region(True, env=self._env)
try:
gcore.run_command(
"r3.in.bin",
flags=flags,
input=self.filename,
output=mapname,
bytes=size,
null=null,
overwrite=overwrite,
quiet=quiet,
north=reg["n"],
south=reg["s"],
top=reg["t"],
bottom=reg["b"],
east=reg["e"],
west=reg["w"],
depths=reg["depths"],
rows=reg["rows3"],
cols=reg["cols3"],
env=self._env,
)
except CalledModuleError:
return 1
else:
return 0