#!/usr/bin/env python ############################################################################ # # MODULE: r.subdayprecip.design # # AUTHOR(S): Martin Landa # # PURPOSE: Computes subday design precipitation totals. # # COPYRIGHT: (C) 2015 Martin Landa and GRASS development team # # This program is free software under the GNU General # Public License (>=v2). Read the file COPYING that # comes with GRASS for details. # ############################################################################# #%module #% description: Computes subday design precipitation totals. #% keyword: raster #% keyword: hydrology #% keyword: precipitation #%end #%option G_OPT_V_MAP #% label: Vector map of location under analysis #%end #%option G_OPT_R_INPUTS #% key: return_period #% description: Rainfall raster maps of required return period #% options: N2,N5,N10,N20,N50,N100 #%end #%option #% key: rainlength #% description: Design rainfall length in minutes #% type: integer #% options: 0-1439 #% required: yes #%end #%option #% key: area_size #% description: Maximum area size to be processed (in km2, -1 for no limit) #% type: double #% answer: 20 #%end import os import sys import grass.script as grass from grass.pygrass.modules import Module from grass.exceptions import CalledModuleError def coeff(name, rl): a = c = None if name == 'N2': if rl < 40: a = 0.166 c = 0.701 elif rl < 120: a = 0.237 c = 0.803 elif rl < 1440: a = 0.235 c = 0.801 elif name == 'N5': if rl < 40: a = 0.171 c = 0.688 elif rl <120: a = 0.265 c = 0.803 elif rl < 1440: a = 0.324 c = 0.845 elif name == 'N10': if rl < 40: a = 0.163 c = 0.656 elif rl <120: a = 0.280 c = 0.803 elif rl < 1440: a = 0.380 c = 0.867 elif name == 'N20': if rl < 40: a = 0.169 c = 0.648 elif rl > 40 and rl < 120: a = 0.300 c = 0.803 elif rl < 1440: a = 0.463 c = 0.894 elif name == 'N50': if rl < 40: a = 0.174 c = 0.638 elif rl < 120: a = 0.323 c = 0.803 elif rl < 1440: a = 0.580 c = 0.925 elif name == 'N100': if rl < 40: a = 0.173 c = 0.625 elif rl < 120: a = 0.335 c = 0.803 elif rl < 1440: a = 0.642 c = 0.939 return a, c def main(): # check if the map is in the current mapset mapset = grass.find_file(opt['map'], element='vector')['mapset'] if not mapset or mapset != grass.gisenv()['MAPSET']: grass.fatal(_("Vector map <{}> not found in the current mapset").format(opt['map'])) # get list of existing columns try: columns = grass.vector_columns(opt['map']).keys() except CalledModuleError as e: return 1 allowed_rasters = ('N2', 'N5', 'N10', 'N20', 'N50', 'N100') # test input feature type vinfo = grass.vector_info_topo(opt['map']) if vinfo['areas'] < 1 and vinfo['points'] < 1: grass.fatal(_("No points or areas found in input vector map <{}>").format(opt['map'])) # check area size limit check_area_size = float(opt['area_size']) > 0 if check_area_size: area_col_name = 'area_{}'.format(os.getpid()) Module('v.db.addcolumn', map=opt['map'], columns='{} double precision'.format(area_col_name)) Module('v.to.db', map=opt['map'], option='area', units='kilometers', columns=area_col_name, quiet=True) areas = Module('v.db.select', flags='c', map=opt['map'], columns=area_col_name, where='{} > {}'.format(area_col_name, opt['area_size']), stdout_=grass.PIPE) large_areas = len(areas.outputs.stdout.splitlines()) if large_areas > 0: grass.warning('{} areas larger than size limit will be skipped from computation'.format(large_areas)) # extract multi values to points for rast in opt['return_period'].split(','): # check valid rasters name = grass.find_file(rast, element='cell')['name'] if not name: grass.warning('Raster map <{}> not found. ' 'Skipped.'.format(rast)) continue if name not in allowed_rasters: grass.warning('Raster map <{}> skipped. ' 'Allowed: {}'.format(rast, allowed_rasters)) continue # perform zonal statistics grass.message('Processing <{}>...'.format(rast)) table = '{}_table'.format(name) if vinfo['areas'] > 0: Module('v.rast.stats', flags='c', map=opt['map'], raster=rast, column_prefix=name, method='average', quiet=True) # handle NULL values (areas smaller than raster resolution) null_values = Module('v.db.select', map=opt['map'], columns='cat', flags='c', where="{}_average is NULL".format(name), stdout_=grass.PIPE) cats = null_values.outputs.stdout.splitlines() if len(cats) > 0: grass.warning(_("Input vector map <{}> contains very small areas (smaller than " "raster resolution). These areas will be proceeded by querying " "single raster cell.").format(opt['map'])) Module('v.what.rast', map=opt['map'], raster=rast, type='centroid', column='{}_average'.format(name), where="{}_average is NULL".format(name), quiet=True) else: # -> points Module('v.what.rast', map=opt['map'], raster=rast, column='{}_average'.format(name), quiet=True) # add column to the attribute table if not exists rl = float(opt['rainlength']) field_name='H_{}T{}'.format(name, opt['rainlength']) if field_name not in columns: Module('v.db.addcolumn', map=opt['map'], columns='{} double precision'.format(field_name)) # determine coefficient for calculation a, c = coeff(rast, rl) if a is None or c is None: grass.fatal("Unable to calculate coefficients") # calculate output values, update attribute table coef = a * rl ** (1 - c) expression = '{}_average * {}'.format(name, coef) Module('v.db.update', map=opt['map'], column=field_name, query_column=expression) if check_area_size: Module('v.db.update', map=opt['map'], column=field_name, value='-1', where='{} > {}'.format(area_col_name, opt['area_size'])) # remove unused column Module('v.db.dropcolumn', map=opt['map'], columns='{}_average'.format(name)) if check_area_size: # remove unused column Module('v.db.dropcolumn', map=opt['map'], columns=area_col_name) return 0 if __name__ == "__main__": opt, flg = grass.parser() sys.exit(main())