forked from OSGeo/grass-addons
-
Notifications
You must be signed in to change notification settings - Fork 0
/
r.raoq.area.py
94 lines (82 loc) · 2.19 KB
/
r.raoq.area.py
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
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Feb 10 09:47:42 2022
@author: lucadelu
"""
############################################################################
#
# MODULE: r.raoq.area
# AUTHOR(S): Luca Delucchi, Michele Torresani
# PURPOSE: It calculates the Rao's Q value for the considered area
#
# COPYRIGHT: (C) 2022 by Luca Delucchi
#
# This program is free software under the GNU General Public
# License (>=v2). Read the file COPYING that comes with GRASS
# for details.
#
#############################################################################
#%module
#% description: Calculate the Rao's Q value for the considered area
#% keyword: raster
#%end
#%option
#% key: input
#% type: string
#% gisprompt: old,cell,raster
#% key_desc: name
#% description: Name of input raster map
#% required: yes
#%end
#%option
#% key: nprocs
#% type: integer
#% required: no
#% multiple: no
#% description: Number of processes which will be used for parallel calculation
#% answer: 1
#%end
#%flag
#% key: g
#% description: Print output value in shell script style
#%end
import sys
import multiprocessing as mp
import numpy as np
import grass.script as grass
array = None
vals = None
def forloop(x):
vals.extend([np.abs(x - y) for y in array.flat])
def main():
global array
global vals
from grass.pygrass.raster import RasterRow
nprocs = int(options["nprocs"])
if nprocs < 1:
grass.fatal("nprocs value should >= 1")
map_in = RasterRow(options["input"])
map_in.open("r")
array = np.array(map_in)
# array[array==-2147483648] = np.NaN
number = np.count_nonzero(array)
number2 = pow(number, 2)
if nprocs == 1:
vals = [np.abs(x - y) for x in array.flat for y in array.flat]
elif nprocs > 1:
manager = mp.Manager()
vals = manager.list()
pool = mp.Pool(nprocs)
listarray = list(array.flat)
pool.map(forloop, listarray)
pool.close()
# pool.join()
out = sum(vals) / number2
if flags["g"]:
print(f"raoq={out}")
else:
print(f"The Rao's Q value is {out}")
if __name__ == "__main__":
options, flags = grass.parser()
sys.exit(main())