| 1 | #!/usr/bin/env python
|
|---|
| 2 |
|
|---|
| 3 | ############################################################################
|
|---|
| 4 | #
|
|---|
| 5 | # MODULE: v.rast.stats
|
|---|
| 6 | # AUTHOR(S): Markus Neteler
|
|---|
| 7 | # converted to Python by Glynn Clements
|
|---|
| 8 | # speed up by Markus Metz
|
|---|
| 9 | # add column choose by Luca Delucchi
|
|---|
| 10 | # PURPOSE: Calculates univariate statistics from a GRASS raster map
|
|---|
| 11 | # only for areas covered by vector objects on a per-category base
|
|---|
| 12 | # COPYRIGHT: (C) 2005-2016 by the GRASS Development Team
|
|---|
| 13 | #
|
|---|
| 14 | # This program is free software under the GNU General Public
|
|---|
| 15 | # License (>=v2). Read the file COPYING that comes with GRASS
|
|---|
| 16 | # for details.
|
|---|
| 17 | #
|
|---|
| 18 | #############################################################################
|
|---|
| 19 |
|
|---|
| 20 | #%module
|
|---|
| 21 | #% description: Calculates univariate statistics from a raster map based on a vector map and uploads statistics to new attribute columns.
|
|---|
| 22 | #% keyword: vector
|
|---|
| 23 | #% keyword: statistics
|
|---|
| 24 | #% keyword: raster
|
|---|
| 25 | #% keyword: univariate statistics
|
|---|
| 26 | #% keyword: zonal statistics
|
|---|
| 27 | #% keyword: sampling
|
|---|
| 28 | #% keyword: querying
|
|---|
| 29 | #%end
|
|---|
| 30 | #%flag
|
|---|
| 31 | #% key: c
|
|---|
| 32 | #% description: Continue if upload column(s) already exist
|
|---|
| 33 | #%end
|
|---|
| 34 | #%flag
|
|---|
| 35 | #% key: d
|
|---|
| 36 | #% label: Create densified lines (default: thin lines)
|
|---|
| 37 | #% description: All cells touched by the line will be set, not only those on the render path
|
|---|
| 38 | #%end
|
|---|
| 39 | #%option G_OPT_V_MAP
|
|---|
| 40 | #%end
|
|---|
| 41 | #%option G_OPT_V_FIELD
|
|---|
| 42 | #%end
|
|---|
| 43 | #%option G_OPT_DB_WHERE
|
|---|
| 44 | #%end
|
|---|
| 45 | #%option G_OPT_R_INPUTS
|
|---|
| 46 | #% key: raster
|
|---|
| 47 | #% description: Name of input raster map to calculate statistics from
|
|---|
| 48 | #%end
|
|---|
| 49 | #%option
|
|---|
| 50 | #% key: column_prefix
|
|---|
| 51 | #% type: string
|
|---|
| 52 | #% description: Column prefix for new attribute columns
|
|---|
| 53 | #% required : yes
|
|---|
| 54 | #% multiple: yes
|
|---|
| 55 | #%end
|
|---|
| 56 | #%option
|
|---|
| 57 | #% key: method
|
|---|
| 58 | #% type: string
|
|---|
| 59 | #% description: The methods to use
|
|---|
| 60 | #% required: no
|
|---|
| 61 | #% multiple: yes
|
|---|
| 62 | #% options: number,null_cells,minimum,maximum,range,average,stddev,variance,coeff_var,sum,first_quartile,median,third_quartile,percentile
|
|---|
| 63 | #% answer: number,null_cells,minimum,maximum,range,average,stddev,variance,coeff_var,sum,first_quartile,median,third_quartile,percentile
|
|---|
| 64 | #%end
|
|---|
| 65 | #%option
|
|---|
| 66 | #% key: percentile
|
|---|
| 67 | #% type: integer
|
|---|
| 68 | #% description: Percentile to calculate
|
|---|
| 69 | #% options: 0-100
|
|---|
| 70 | #% answer: 90
|
|---|
| 71 | #% required : no
|
|---|
| 72 | #%end
|
|---|
| 73 |
|
|---|
| 74 | import sys
|
|---|
| 75 | import os
|
|---|
| 76 | import atexit
|
|---|
| 77 | import grass.script as grass
|
|---|
| 78 | from grass.script.utils import decode
|
|---|
| 79 | from grass.exceptions import CalledModuleError
|
|---|
| 80 |
|
|---|
| 81 |
|
|---|
| 82 | def cleanup():
|
|---|
| 83 | if rastertmp:
|
|---|
| 84 | grass.run_command('g.remove', flags='f', type='raster',
|
|---|
| 85 | name=rastertmp, quiet=True)
|
|---|
| 86 | # for f in [tmp, tmpname, sqltmp]:
|
|---|
| 87 | # grass.try_remove(f)
|
|---|
| 88 |
|
|---|
| 89 |
|
|---|
| 90 | def main():
|
|---|
| 91 | global tmp, sqltmp, tmpname, nuldev, vector, rastertmp
|
|---|
| 92 | rastertmp = False
|
|---|
| 93 | # setup temporary files
|
|---|
| 94 | tmp = grass.tempfile()
|
|---|
| 95 | sqltmp = tmp + ".sql"
|
|---|
| 96 | # we need a random name
|
|---|
| 97 | tmpname = grass.basename(tmp)
|
|---|
| 98 |
|
|---|
| 99 | nuldev = open(os.devnull, 'w')
|
|---|
| 100 |
|
|---|
| 101 | rasters = options['raster'].split(',')
|
|---|
| 102 | colprefixes = options['column_prefix'].split(',')
|
|---|
| 103 | vector = options['map']
|
|---|
| 104 | layer = options['layer']
|
|---|
| 105 | where = options['where']
|
|---|
| 106 | percentile = options['percentile']
|
|---|
| 107 | basecols = options['method'].split(',')
|
|---|
| 108 |
|
|---|
| 109 | ### setup enviro vars ###
|
|---|
| 110 | env = grass.gisenv()
|
|---|
| 111 | mapset = env['MAPSET']
|
|---|
| 112 |
|
|---|
| 113 | vs = vector.split('@')
|
|---|
| 114 | if len(vs) > 1:
|
|---|
| 115 | vect_mapset = vs[1]
|
|---|
| 116 | else:
|
|---|
| 117 | vect_mapset = mapset
|
|---|
| 118 |
|
|---|
| 119 | # does map exist in CURRENT mapset?
|
|---|
| 120 | if vect_mapset != mapset or not grass.find_file(vector, 'vector', mapset)['file']:
|
|---|
| 121 | grass.fatal(_("Vector map <%s> not found in current mapset") % vector)
|
|---|
| 122 |
|
|---|
| 123 | # check if DBF driver used, in this case cut to 10 chars col names:
|
|---|
| 124 | try:
|
|---|
| 125 | fi = grass.vector_db(map=vector)[int(layer)]
|
|---|
| 126 | except KeyError:
|
|---|
| 127 | grass.fatal(
|
|---|
| 128 | _('There is no table connected to this map. Run v.db.connect or v.db.addtable first.'))
|
|---|
| 129 | # we need this for non-DBF driver:
|
|---|
| 130 | dbfdriver = fi['driver'] == 'dbf'
|
|---|
| 131 |
|
|---|
| 132 | # colprefix for every raster map?
|
|---|
| 133 | if len(colprefixes) != len(rasters):
|
|---|
| 134 | grass.fatal(_("Number of raster maps ({0}) different from \
|
|---|
| 135 | number of column prefixes ({1})". format(len(rasters),
|
|---|
| 136 | len(colprefixes))))
|
|---|
| 137 |
|
|---|
| 138 | vector = vs[0]
|
|---|
| 139 |
|
|---|
| 140 | rastertmp = "%s_%s" % (vector, tmpname)
|
|---|
| 141 |
|
|---|
| 142 | for raster in rasters:
|
|---|
| 143 | # check the input raster map
|
|---|
| 144 | if not grass.find_file(raster, 'cell')['file']:
|
|---|
| 145 | grass.fatal(_("Raster map <%s> not found") % raster)
|
|---|
| 146 |
|
|---|
| 147 | # save current settings:
|
|---|
| 148 | grass.use_temp_region()
|
|---|
| 149 |
|
|---|
| 150 | # Temporarily aligning region resolution to $RASTER resolution
|
|---|
| 151 | # keep boundary settings
|
|---|
| 152 | grass.run_command('g.region', align=rasters[0])
|
|---|
| 153 |
|
|---|
| 154 | # prepare base raster for zonal statistics
|
|---|
| 155 | try:
|
|---|
| 156 | nlines = grass.vector_info_topo(vector)['lines']
|
|---|
| 157 | kwargs = {}
|
|---|
| 158 | if where:
|
|---|
| 159 | kwargs['where'] = where
|
|---|
| 160 | # Create densified lines rather than thin lines
|
|---|
| 161 | if flags['d'] and nlines > 0:
|
|---|
| 162 | kwargs['flags'] = 'd'
|
|---|
| 163 |
|
|---|
| 164 | grass.run_command('v.to.rast', input=vector, layer=layer, output=rastertmp,
|
|---|
| 165 | use='cat', quiet=True, **kwargs)
|
|---|
| 166 | except CalledModuleError:
|
|---|
| 167 | grass.fatal(_("An error occurred while converting vector to raster"))
|
|---|
| 168 |
|
|---|
| 169 | # dump cats to file to avoid "too many argument" problem:
|
|---|
| 170 | p = grass.pipe_command('r.category', map=rastertmp, sep=';', quiet=True)
|
|---|
| 171 | cats = []
|
|---|
| 172 |
|
|---|
| 173 | for line in p.stdout:
|
|---|
| 174 | line = decode(line)
|
|---|
| 175 | cats.append(line.rstrip('\r\n').split(';')[0])
|
|---|
| 176 | p.wait()
|
|---|
| 177 |
|
|---|
| 178 | number = len(cats)
|
|---|
| 179 | if number < 1:
|
|---|
| 180 | grass.fatal(_("No categories found in raster map"))
|
|---|
| 181 |
|
|---|
| 182 | # Check if all categories got converted
|
|---|
| 183 | # Report categories from vector map
|
|---|
| 184 | vect_cats = grass.read_command('v.category', input=vector, option='report',
|
|---|
| 185 | flags='g').rstrip('\n').split('\n')
|
|---|
| 186 |
|
|---|
| 187 | # get number of all categories in selected layer
|
|---|
| 188 | for vcl in vect_cats:
|
|---|
| 189 | if vcl.split(' ')[0] == layer and vcl.split(' ')[1] == 'all':
|
|---|
| 190 | vect_cats_n = int(vcl.split(' ')[2])
|
|---|
| 191 |
|
|---|
| 192 | if vect_cats_n != number:
|
|---|
| 193 | grass.warning(_("Not all vector categories converted to raster. \
|
|---|
| 194 | Converted {0} of {1}.".format(number, vect_cats_n)))
|
|---|
| 195 |
|
|---|
| 196 | # check if DBF driver used, in this case cut to 10 chars col names:
|
|---|
| 197 | try:
|
|---|
| 198 | fi = grass.vector_db(map=vector)[int(layer)]
|
|---|
| 199 | except KeyError:
|
|---|
| 200 | grass.fatal(
|
|---|
| 201 | _('There is no table connected to this map. Run v.db.connect or v.db.addtable first.'))
|
|---|
| 202 | # we need this for non-DBF driver:
|
|---|
| 203 | dbfdriver = fi['driver'] == 'dbf'
|
|---|
| 204 |
|
|---|
| 205 | # Find out which table is linked to the vector map on the given layer
|
|---|
| 206 | if not fi['table']:
|
|---|
| 207 | grass.fatal(
|
|---|
| 208 | _('There is no table connected to this map. Run v.db.connect or v.db.addtable first.'))
|
|---|
| 209 |
|
|---|
| 210 | # replaced by user choiche
|
|---|
| 211 | #basecols = ['n', 'min', 'max', 'range', 'mean', 'stddev', 'variance', 'cf_var', 'sum']
|
|---|
| 212 |
|
|---|
| 213 | for i in range(len(rasters)):
|
|---|
| 214 | raster = rasters[i]
|
|---|
| 215 | colprefix = colprefixes[i]
|
|---|
| 216 | # we need at least three chars to distinguish [mea]n from [med]ian
|
|---|
| 217 | # so colprefix can't be longer than 6 chars with DBF driver
|
|---|
| 218 | if dbfdriver:
|
|---|
| 219 | colprefix = colprefix[:6]
|
|---|
| 220 | variables_dbf = {}
|
|---|
| 221 |
|
|---|
| 222 | # by default perccol variable is used only for "variables" variable
|
|---|
| 223 | perccol = "percentile"
|
|---|
| 224 | perc = None
|
|---|
| 225 | for b in basecols:
|
|---|
| 226 | if b.startswith('p'):
|
|---|
| 227 | perc = b
|
|---|
| 228 | if perc:
|
|---|
| 229 | # namespace is limited in DBF but the % value is important
|
|---|
| 230 | if dbfdriver:
|
|---|
| 231 | perccol = "per" + percentile
|
|---|
| 232 | else:
|
|---|
| 233 | perccol = "percentile_" + percentile
|
|---|
| 234 | percindex = basecols.index(perc)
|
|---|
| 235 | basecols[percindex] = perccol
|
|---|
| 236 |
|
|---|
| 237 | # dictionary with name of methods and position in "r.univar -gt" output
|
|---|
| 238 | variables = {'number': 2, 'null_cells': 2, 'minimum': 4, 'maximum': 5, 'range': 6,
|
|---|
| 239 | 'average': 7, 'stddev': 9, 'variance': 10, 'coeff_var': 11,
|
|---|
| 240 | 'sum': 12, 'first_quartile': 14, 'median': 15,
|
|---|
| 241 | 'third_quartile': 16, perccol: 17}
|
|---|
| 242 | # this list is used to set the 'e' flag for r.univar
|
|---|
| 243 | extracols = ['first_quartile', 'median', 'third_quartile', perccol]
|
|---|
| 244 | addcols = []
|
|---|
| 245 | colnames = []
|
|---|
| 246 | extstat = ""
|
|---|
| 247 | for i in basecols:
|
|---|
| 248 | # this check the complete name of out input that should be truncated
|
|---|
| 249 | for k in variables.keys():
|
|---|
| 250 | if i in k:
|
|---|
| 251 | i = k
|
|---|
| 252 | break
|
|---|
| 253 | if i in extracols:
|
|---|
| 254 | extstat = 'e'
|
|---|
| 255 | # check if column already present
|
|---|
| 256 | currcolumn = ("%s_%s" % (colprefix, i))
|
|---|
| 257 | if dbfdriver:
|
|---|
| 258 | currcolumn = currcolumn[:10]
|
|---|
| 259 | variables_dbf[currcolumn.replace("%s_" % colprefix, '')] = i
|
|---|
| 260 |
|
|---|
| 261 | colnames.append(currcolumn)
|
|---|
| 262 | if currcolumn in grass.vector_columns(vector, layer).keys():
|
|---|
| 263 | if not flags['c']:
|
|---|
| 264 | grass.fatal((_("Cannot create column <%s> (already present). ") % currcolumn) +
|
|---|
| 265 | _("Use -c flag to update values in this column."))
|
|---|
| 266 | else:
|
|---|
| 267 | if i == "n":
|
|---|
| 268 | coltype = "INTEGER"
|
|---|
| 269 | else:
|
|---|
| 270 | coltype = "DOUBLE PRECISION"
|
|---|
| 271 | addcols.append(currcolumn + ' ' + coltype)
|
|---|
| 272 |
|
|---|
| 273 | if addcols:
|
|---|
| 274 | grass.verbose(_("Adding columns '%s'") % addcols)
|
|---|
| 275 | try:
|
|---|
| 276 | grass.run_command('v.db.addcolumn', map=vector, columns=addcols,
|
|---|
| 277 | layer=layer)
|
|---|
| 278 | except CalledModuleError:
|
|---|
| 279 | grass.fatal(_("Adding columns failed. Exiting."))
|
|---|
| 280 |
|
|---|
| 281 | # calculate statistics:
|
|---|
| 282 | grass.message(_("Processing input data (%d categories)...") % number)
|
|---|
| 283 |
|
|---|
| 284 | # get rid of any earlier attempts
|
|---|
| 285 | grass.try_remove(sqltmp)
|
|---|
| 286 |
|
|---|
| 287 | f = open(sqltmp, 'w')
|
|---|
| 288 |
|
|---|
| 289 | # do the stats
|
|---|
| 290 | p = grass.pipe_command('r.univar', flags='t' + extstat, map=raster,
|
|---|
| 291 | zones=rastertmp, percentile=percentile, sep=';')
|
|---|
| 292 |
|
|---|
| 293 | first_line = 1
|
|---|
| 294 |
|
|---|
| 295 | f.write("{0}\n".format(grass.db_begin_transaction(fi['driver'])))
|
|---|
| 296 | for line in p.stdout:
|
|---|
| 297 | if first_line:
|
|---|
| 298 | first_line = 0
|
|---|
| 299 | continue
|
|---|
| 300 |
|
|---|
| 301 | vars = decode(line).rstrip('\r\n').split(';')
|
|---|
| 302 |
|
|---|
| 303 | f.write("UPDATE %s SET" % fi['table'])
|
|---|
| 304 | first_var = 1
|
|---|
| 305 | for colname in colnames:
|
|---|
| 306 | variable = colname.replace("%s_" % colprefix, '', 1)
|
|---|
| 307 | if dbfdriver:
|
|---|
| 308 | variable = variables_dbf[variable]
|
|---|
| 309 | i = variables[variable]
|
|---|
| 310 | value = vars[i]
|
|---|
| 311 | # convert nan, +nan, -nan, inf, +inf, -inf, Infinity, +Infinity,
|
|---|
| 312 | # -Infinity to NULL
|
|---|
| 313 | if value.lower().endswith('nan') or 'inf' in value.lower():
|
|---|
| 314 | value = 'NULL'
|
|---|
| 315 | if not first_var:
|
|---|
| 316 | f.write(" , ")
|
|---|
| 317 | else:
|
|---|
| 318 | first_var = 0
|
|---|
| 319 | f.write(" %s=%s" % (colname, value))
|
|---|
| 320 |
|
|---|
| 321 | f.write(" WHERE %s=%s;\n" % (fi['key'], vars[0]))
|
|---|
| 322 | f.write("{0}\n".format(grass.db_commit_transaction(fi['driver'])))
|
|---|
| 323 | p.wait()
|
|---|
| 324 | f.close()
|
|---|
| 325 |
|
|---|
| 326 | grass.message(_("Updating the database ..."))
|
|---|
| 327 | exitcode = 0
|
|---|
| 328 | try:
|
|---|
| 329 | grass.run_command('db.execute', input=sqltmp,
|
|---|
| 330 | database=fi['database'], driver=fi['driver'])
|
|---|
| 331 | grass.verbose((_("Statistics calculated from raster map <{raster}>"
|
|---|
| 332 | " and uploaded to attribute table"
|
|---|
| 333 | " of vector map <{vector}>."
|
|---|
| 334 | ).format(raster=raster, vector=vector)))
|
|---|
| 335 | except CalledModuleError:
|
|---|
| 336 | grass.warning(
|
|---|
| 337 | _("Failed to upload statistics to attribute table of vector map <%s>.") %
|
|---|
| 338 | vector)
|
|---|
| 339 | exitcode = 1
|
|---|
| 340 |
|
|---|
| 341 | sys.exit(exitcode)
|
|---|
| 342 |
|
|---|
| 343 | if __name__ == "__main__":
|
|---|
| 344 | options, flags = grass.parser()
|
|---|
| 345 | atexit.register(cleanup)
|
|---|
| 346 | main()
|
|---|
