← Back to branch summary

~ubuntu-branches/ubuntu/wily/grass/wily

« back to all changes in this revision

Viewing changes to scripts/r.reclass.area/r.reclass.area.py

Tags: 7.0.0~rc1+ds1-1~exp1
http://bugs.debian.org/776349
* New upstream release candidate.
* Repack upstream tarball, remove precompiled Python objects.
* Add upstream metadata.
* Update gbp.conf and Vcs-Git URL to use the experimental branch.
* Update watch file for GRASS 7.0.
* Drop build dependencies for Tcl/Tk, add build dependencies:
  python-numpy, libnetcdf-dev, netcdf-bin, libblas-dev, liblapack-dev
* Update Vcs-Browser URL to use cgit instead of gitweb.
* Update paths to use grass70.
* Add configure options: --with-netcdf, --with-blas, --with-lapack,
  remove --with-tcltk-includes.
* Update patches for GRASS 7.
* Update copyright file, changes:
  - Update copyright years
  - Group files by license
  - Remove unused license sections
* Add patches for various typos.
* Fix desktop file with patch instead of d/rules.
* Use minimal dh rules.
* Bump Standards-Version to 3.9.6, no changes.
* Use dpkg-maintscript-helper to replace directories with symlinks.
  (closes: #776349)
* Update my email to use @debian.org address.

Show diffs side-by-side

added added

removed removed

Lines of Context:
scripts/r.reclass.area/r.reclass.area.py
 
1
#!/usr/bin/env python
 
2
 
 
3
############################################################################
 
4
#
 
5
# MODULE:       r.reclass.area
 
6
# AUTHOR(S):    NRCS
 
7
#               Converted to Python by Glynn Clements
 
8
#               Added rmarea method by Luca Delucchi
 
9
# PURPOSE:      Reclasses a raster map greater or less than user specified area size (in hectares)
 
10
# COPYRIGHT:    (C) 1999,2008,2014 by the GRASS Development Team
 
11
#
 
12
#               This program is free software under the GNU General Public
 
13
#               License (>=v2). Read the file COPYING that comes with GRASS
 
14
#               for details.
 
15
#
 
16
#############################################################################
 
17
# 10/2013: added option to use a pre-clumped input map (Eric Goddard)
 
18
# 8/2012: added fp maps support, cleanup, removed tabs AK
 
19
# 3/2007: added label support MN
 
20
# 3/2004: added parser support MN
 
21
# 11/2001 added mapset support markus
 
22
# 2/2001 fixes markus
 
23
# 2000: updated to GRASS 5
 
24
# 1998 from NRCS, slightly modified for GRASS 4.2.1
 
25
 
 
26
#%module
 
27
#% description: Reclasses a raster map greater or less than user specified area size (in hectares).
 
28
#% keyword: raster
 
29
#% keyword: statistics
 
30
#% keyword: aggregation
 
31
#%end
 
32
 
 
33
#%option G_OPT_R_INPUT
 
34
#%end
 
35
 
 
36
#%option G_OPT_R_OUTPUT
 
37
#%end
 
38
 
 
39
#%option
 
40
#% key: value
 
41
#% type: double
 
42
#% description: Value option that sets the area size limit (in hectares)
 
43
#% required: yes
 
44
#% guisection: Area
 
45
#%end
 
46
 
 
47
#%option
 
48
#% key: mode
 
49
#% type: string
 
50
#% description: Lesser or greater than specified value
 
51
#% options: lesser,greater
 
52
#% required: yes
 
53
#% guisection: Area
 
54
#%end
 
55
 
 
56
#%option
 
57
#% key: method
 
58
#% type: string
 
59
#% description: Method used for reclassification
 
60
#% options: reclass,rmarea
 
61
#% answer: reclass
 
62
#% guisection: Area
 
63
#%end
 
64
 
 
65
#%flag
 
66
#% key: c
 
67
#% description: Input map is clumped
 
68
#%end
 
69
 
 
70
#%flag
 
71
#% key: d
 
72
#% description: Clumps including diagonal neighbors
 
73
#%end
 
74
 
 
75
import sys
 
76
import os
 
77
import atexit
 
78
import grass.script as grass
 
79
 
 
80
TMPRAST = []
 
81
 
 
82
 
 
83
def reclass(inf, outf, lim, clump, diag, les):
 
84
    infile = inf
 
85
    outfile = outf
 
86
    lesser = les
 
87
    limit = lim
 
88
    clumped = clump
 
89
    diagonal = diag
 
90
 
 
91
    s = grass.read_command("g.region", flags='p')
 
92
    kv = grass.parse_key_val(s, sep=':')
 
93
    s = kv['projection'].strip().split()
 
94
    if s == '0':
 
95
        grass.fatal(_("xy-locations are not supported"))
 
96
        grass.fatal(_("Need projected data with grids in meters"))
 
97
 
 
98
    if not grass.find_file(infile)['name']:
 
99
        grass.fatal(_("Raster map <%s> not found") % infile)
 
100
 
 
101
    if clumped and diagonal:
 
102
        grass.fatal(_("flags c and d are mutually exclusive"))
 
103
 
 
104
    if clumped:
 
105
        clumpfile = infile
 
106
    else:
 
107
        clumpfile = "%s.clump.%s" % (infile.split('@')[0], outfile)
 
108
        TMPRAST.append(clumpfile)
 
109
 
 
110
        if not grass.overwrite():
 
111
            if grass.find_file(clumpfile)['name']:
 
112
                grass.fatal(_("Temporary raster map <%s> exists") % clumpfile)
 
113
        if diagonal:
 
114
            grass.message(_("Generating a clumped raster file including "
 
115
                            "diagonal neighbors..."))
 
116
            grass.run_command('r.clump', flags='d', input=infile,
 
117
                              output=clumpfile)
 
118
        else:
 
119
            grass.message(_("Generating a clumped raster file ..."))
 
120
            grass.run_command('r.clump', input=infile, output=clumpfile)
 
121
 
 
122
    if lesser:
 
123
        grass.message(_("Generating a reclass map with area size less than "
 
124
                        "or equal to %f hectares...") % limit)
 
125
    else:
 
126
        grass.message(_("Generating a reclass map with area size greater "
 
127
                        "than or equal to %f hectares...") % limit)
 
128
 
 
129
    recfile = outfile + '.recl'
 
130
    TMPRAST.append(recfile)
 
131
 
 
132
    sflags = 'aln'
 
133
    if grass.raster_info(infile)['datatype'] in ('FCELL', 'DCELL'):
 
134
        sflags += 'i'
 
135
    p1 = grass.pipe_command('r.stats', flags=sflags, input=(clumpfile, infile),
 
136
                            sep=';')
 
137
    p2 = grass.feed_command('r.reclass', input=clumpfile, output=recfile,
 
138
                            rules='-')
 
139
    rules = ''
 
140
    for line in p1.stdout:
 
141
        f = line.rstrip(os.linesep).split(';')
 
142
        if len(f) < 5:
 
143
            continue
 
144
        hectares = float(f[4]) * 0.0001
 
145
        if lesser:
 
146
            test = hectares <= limit
 
147
        else:
 
148
            test = hectares >= limit
 
149
        if test:
 
150
            rules += "%s = %s %s\n" % (f[0], f[2], f[3])
 
151
    if rules:
 
152
        p2.stdin.write(rules)
 
153
    p1.wait()
 
154
    p2.stdin.close()
 
155
    p2.wait()
 
156
    if p2.returncode != 0:
 
157
        if lesser:
 
158
            grass.fatal(_("No areas of size less than or equal to %f "
 
159
                          "hectares found.") % limit)
 
160
        else:
 
161
            grass.fatal(_("No areas of size greater than or equal to %f "
 
162
                          "hectares found.") % limit)
 
163
    grass.mapcalc("$outfile = $recfile", outfile=outfile, recfile=recfile)
 
164
 
 
165
 
 
166
def rmarea(infile, outfile, thresh, coef):
 
167
    # transform user input from hectares to map units (kept this for future)
 
168
    # thresh = thresh * 10000.0 / (float(coef)**2)
 
169
    # grass.debug("Threshold: %d, coeff linear: %s, coef squared: %d" % (thresh, coef, (float(coef)**2)), 0)
 
170
 
 
171
    # transform user input from hectares to meters because currently v.clean
 
172
    # rmarea accept only meters as threshold
 
173
    thresh = thresh * 10000.0
 
174
    vectfile = "%s_vect_%s" % (infile.split('@')[0], outfile)
 
175
    TMPRAST.append(vectfile)
 
176
    grass.run_command('r.to.vect', input=infile, output=vectfile, type='area')
 
177
    cleanfile = "%s_clean_%s" % (infile.split('@')[0], outfile)
 
178
    TMPRAST.append(cleanfile)
 
179
    grass.run_command('v.clean', input=vectfile, output=cleanfile,
 
180
                      tool='rmarea', threshold=thresh)
 
181
 
 
182
    grass.run_command('v.to.rast', input=cleanfile, output=outfile,
 
183
                      use='attr', attrcolumn='value')
 
184
 
 
185
 
 
186
def main():
 
187
    infile = options['input']
 
188
    value = options['value']
 
189
    mode = options['mode']
 
190
    outfile = options['output']
 
191
    global method
 
192
    method = options['method']
 
193
    clumped = flags['c']
 
194
    diagonal = flags['d']
 
195
 
 
196
    # check for unsupported locations
 
197
    in_proj = grass.parse_command('g.proj', flags='g')
 
198
    if in_proj['unit'].lower() == 'degree':
 
199
        grass.fatal(_("Latitude-longitude locations are not supported"))
 
200
    if in_proj['name'].lower() == 'xy_location_unprojected':
 
201
        grass.fatal(_("xy-locations are not supported"))
 
202
 
 
203
    # check lesser and greater parameters
 
204
    limit = float(value)
 
205
    if mode == 'greater' and method == 'rmarea':
 
206
        grass.fatal(_("You have to specify mode='lesser' with method='rmarea'"))
 
207
    
 
208
    if not grass.find_file(infile)['name']:
 
209
        grass.fatal(_("Raster map <%s> not found") % infile)
 
210
 
 
211
    if method == 'reclass':
 
212
        reclass(infile, outfile, limit, clumped, diagonal, mode == 'lesser')
 
213
    elif method == 'rmarea':
 
214
        rmarea(infile, outfile, limit, in_proj['meters'])
 
215
 
 
216
    grass.message(_("Generating output raster map <%s>...") % outfile)
 
217
 
 
218
 
 
219
def cleanup():
 
220
    """!Delete temporary maps"""
 
221
    TMPRAST.reverse()  # reclassed map first
 
222
    for mapp in TMPRAST:
 
223
        if method == 'rmarea':
 
224
            grass.run_command("g.remove", flags='f', type='vector', name=mapp,
 
225
                              quiet=True)
 
226
        else:
 
227
            grass.run_command("g.remove", flags='f', type='raster', name=mapp,
 
228
                              quiet=True)
 
229
 
 
230
if __name__ == "__main__":
 
231
    options, flags = grass.parser()
 
232
    atexit.register(cleanup)
 
233
    sys.exit(main())