[GRASS-SVN] r52231 - grass-addons/grass7/vector/v.surf.icw
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Jun 26 04:55:53 PDT 2012
Author: hamish
Date: 2012-06-26 04:55:52 -0700 (Tue, 26 Jun 2012)
New Revision: 52231
Modified:
grass-addons/grass7/vector/v.surf.icw/v.surf.icw.py
Log:
enable full parallelization
Modified: grass-addons/grass7/vector/v.surf.icw/v.surf.icw.py
===================================================================
--- grass-addons/grass7/vector/v.surf.icw/v.surf.icw.py 2012-06-26 11:55:10 UTC (rev 52230)
+++ grass-addons/grass7/vector/v.surf.icw/v.surf.icw.py 2012-06-26 11:55:52 UTC (rev 52231)
@@ -125,7 +125,8 @@
def cleanup():
grass.verbose(_("Cleanup.."))
tmp_base = 'tmp_icw_' + str(os.getpid()) + '_'
- grass.run_command('g.mremove', flags = 'f', rast = tmp_base + '*', quiet = True)
+ grass.run_command('g.mremove', flags = 'f', rast = tmp_base + '*',
+ quiet = True)
def main():
@@ -142,6 +143,8 @@
if workers is 1 and "WORKERS" in os.environ:
workers = int(os.environ["WORKERS"])
+ if workers < 1:
+ workers = 1
pid = str(os.getpid())
tmp_base = 'tmp_icw_' + pid + '_'
@@ -156,7 +159,8 @@
grass.fatal(_("Raster map <%s> not found") % post_mask)
grass.verbose(_("v.surf.icw -- Inverse Cost Weighted Interpolation"))
- grass.verbose(_("Processing %s -> %s, column=%s, Cf=%g") % (pts_input, output, column, friction))
+ grass.verbose(_("Processing %s -> %s, column=%s, Cf=%g")
+ % (pts_input, output, column, friction))
if flags['r']:
grass.verbose(_("Using (d^n)*log(d) radial basis function."))
@@ -176,7 +180,8 @@
try:
coltype = grass.vector_columns(pts_input, layer)[column]
except KeyError:
- grass.fatal(_("Data column <%s> not found in vector points map <%s>") % (column, pts_input))
+ grass.fatal(_("Data column <%s> not found in vector points map <%s>")
+ % (column, pts_input))
if coltype['type'] not in ('INTEGER', 'DOUBLE PRECISION'):
grass.fatal(_("Data column must be numberic"))
@@ -184,7 +189,8 @@
# cleanse cost area mask to a flat =1 for my porpoises
area_mask = tmp_base + 'area'
grass.mapcalc("$result = if($cost_map, 1, null())",
- result = area_mask, cost_map = cost_map)
+ result = area_mask, cost_map = cost_map,
+ quiet = True)
## done with prep work,
########################################################################
@@ -217,7 +223,7 @@
#### generate cost maps for each site in range
- grass.message(_("Generating cost maps.."))
+ grass.message(_("Generating cost maps ..."))
# avoid do-it-yourself brain surgery
points_list_orig = list(points_list)
@@ -241,8 +247,9 @@
n -= 1
continue
else:
- grass.message(_("Site %d of %d, e=%.4f n=%.4f cat=%d data=%.8g") % (num, n,
- float(easting), float(northing), cat, float(data_value)))
+ grass.message(_("Site %d of %d, e=%.4f n=%.4f cat=%d data=%.8g")
+ % (num, n, float(easting), float(northing), cat,
+ float(data_value)))
# we know the point is in the region, but is it in a non-null area of the cost surface?
rast_val = grass.read_command('r.what', input = area_mask,
@@ -261,33 +268,56 @@
grass.fatal('Data value [%s] is non-numeric' % data_value)
cost_site_name = tmp_base + 'cost_site.' + '%05d' % num
- ret = grass.run_command('r.cost', flags = 'k', input = area_mask,
- output = cost_site_name,
- coordinate = easting + ',' + northing, quiet = True)
- if ret is not 0:
- grass.fatal(_('Problem running %s') % 'r.cost')
+ proc[num-1] = grass.start_command('r.cost', flags = 'k', input = area_mask,
+ output = cost_site_name,
+ coordinate = easting + ',' + northing,
+ quiet = True)
+ # stall to wait for the nth worker to complete,
+ if num % workers is 0:
+ proc[num-1].wait()
+ num += 1
+ # make sure everyone is finished
+ for i in range(n):
+ if proc[i].wait() is not 0:
+ grass.fatal(_('Problem running %s') % 'r.cost')
+ grass.message(_("Removing anomalies at site positions ..."))
+ proc = {}
+ for i in range(n):
+ cost_site_name = tmp_base + 'cost_site.' + '%05d' % (i+1)
#max_cost="$GIS_OPT_MAX_COST" : commented out until r.null cleansing/continue code is sorted out
#start_points=tmp_idw_cost_site_$$
# we do this so the divisor exists and the weighting is huge at the exact sample spots
# more efficient to reclass to 1?
- grass.mapcalc("$cost_n_cleansed = if($cost_n == 0, 0.1, $cost_n)",
+ proc[i] = grass.mapcalc_start("$cost_n_cleansed = if($cost_n == 0, 0.1, $cost_n)",
cost_n_cleansed = cost_site_name + '.cleansed',
- cost_n = cost_site_name)
+ cost_n = cost_site_name, quiet = True)
+ # stall to wait for the nth worker to complete,
+ if i % workers is 0:
+ #print 'stalling ...'
+ proc[i].wait()
+ # make sure everyone is finished
+ for i in range(n):
+ if proc[i].wait() is not 0:
+ grass.fatal(_('Problem running %s') % 'r.mapcalc')
+
+
+ grass.message(_("Applying radial decay ..."))
+
+ proc = {}
+ for i in range(n):
+ cost_site_name = tmp_base + 'cost_site.' + '%05d' % (i+1)
grass.run_command('g.remove', rast = cost_site_name, quiet = True)
grass.run_command('g.rename', rast = cost_site_name + '.cleansed'
+ ',' + cost_site_name, quiet = True)
-
-
-
# r.to.vect then r.patch output
# v.to.rast in=tmp_idw_cost_site_29978 out=tmp_idw_cost_val_$$ use=val val=10
@@ -303,31 +333,35 @@
grass.debug("r.mapcalc expression is: [%s]" % expr)
- one_by_cost_site_sq_n = tmp_base + '1by_cost_site_sq.' + '%05d' % num
+ one_by_cost_site_sq_n = tmp_base + '1by_cost_site_sq.' + '%05d' % (i+1)
- proc[num-1] = grass.mapcalc_start("$result = " + expr,
+ proc[i] = grass.mapcalc_start("$result = " + expr,
result = one_by_cost_site_sq_n,
cost_n = cost_site_name,
- friction = friction)
-
+ friction = friction,
+ quiet = True)
# stall to wait for the nth worker to complete,
- if num % workers is 0:
- proc[num-1].wait()
+ if i % workers is 0:
+ #print 'stalling ...'
+ proc[i].wait()
# r.patch in=1by_cost_site_sqrd.${NUM},tmp_idw_cost_val_$$ out=1by_cost_site_sqrd.${NUM} --o
# g.remove rast=cost_site.$NUM
- num += 1
-
# make sure everyone is finished
for i in range(n):
- proc[i].wait()
+ if proc[i].wait() is not 0:
+ grass.fatal(_('Problem running %s') % 'r.mapcalc')
+ grass.run_command('g.mremove', flags = 'f',
+ rast = tmp_base + 'cost_site.*', quiet = True)
+ #grass.run_command('g.list', type = 'rast', mapset = '.')
+
#######################################################
#### Step 3) find sum(cost^2)
- grass.message("")
- grass.message(_("Finding sum of squares.."))
+ grass.verbose('')
+ grass.verbose(_("Finding sum of squares ..."))
#todo: test if MASK exists already, fatal exit if it does?
if post_mask:
@@ -335,15 +369,19 @@
grass.mapcalc("MASK = $maskmap", maskmap = post_mask, overwrite = True)
- grass.message(_("Summation of cost weights.."))
+ grass.message(_("Summation of cost weights ..."))
input_maps = tmp_base + '1by_cost_site_sq.%05d' % 1
for i in range(2, n+1):
input_maps += ',%s1by_cost_site_sq.%05d' % (tmp_base, i)
+ #grass.run_command('g.list', type = 'rast', mapset = '.')
+
sum_of_1by_cost_sqs = tmp_base + 'sum_of_1by_cost_sqs'
- grass.run_command('r.series', method = 'sum', input = input_maps,
+ ret = grass.run_command('r.series', method = 'sum', input = input_maps,
output = sum_of_1by_cost_sqs)
+ if ret is not 0:
+ grass.fatal(_('Problem running %s') % 'r.series')
if post_mask:
grass.message(_("Removing post_mask <%s>"), post_mask)
@@ -352,10 +390,9 @@
#######################################################
#### Step 4) ( 1/di^2 / sum(1/d^2) ) * ai
- grass.message("")
- grass.message(_("Creating partial weights.."))
+ grass.verbose('')
+ grass.message(_("Creating partial weights ..."))
-
proc = {}
num = 1
for position in points_list:
@@ -376,8 +413,8 @@
% (num, n, cat, data_value))
# we know the point is in the region, but is it in a non-null area of the cost surface?
- rast_val = grass.read_command('r.what', input = area_mask,
- east_north = '%s,%s' % (position[0], position[1])
+ rast_val = grass.read_command('r.what', map = area_mask,
+ coordinates = '%s,%s' % (position[0], position[1])
).strip().split('|')[-1]
if rast_val == '*':
grass.message(_(" -- Skipping, point lays outside of cost_map. [Probably programmer error]"))
@@ -390,11 +427,13 @@
#"( $DATA_VALUE / $N ) * (1.0 - ( cost_sq_site.$NUM / sum_of_cost_sqs ))"
#"( cost_sq_site.$NUM / sum_of_cost_sqs ) * ( $DATA_VALUE / $N )"
- proc[num-1] = grass.mapcalc_start("$partial_n = ($data * $one_by_cost_sq) / $sum_of_1by_cost_sqs",
+ proc[num-1] = grass.mapcalc_start(
+ "$partial_n = ($data * $one_by_cost_sq) / $sum_of_1by_cost_sqs",
partial_n = partial_n,
data = data_value,
one_by_cost_sq = one_by_cost_site_sq,
- sum_of_1by_cost_sqs = sum_of_1by_cost_sqs)
+ sum_of_1by_cost_sqs = sum_of_1by_cost_sqs,
+ quiet = True)
# stall to wait for the nth worker to complete,
if num % workers is 0:
@@ -411,10 +450,16 @@
for i in range(n):
proc[i].wait()
+ # free up disk space ASAP
+ grass.run_command('g.mremove', flags = 'f',
+ rast = tmp_base + '1by_cost_site_sq.*', quiet = True)
+ #grass.run_command('g.list', type = 'rast', mapset = '.')
- grass.message("")
- grass.message(_("Calculating final values.."))
+ #######################################################
+ grass.message('')
+ grass.message(_("Calculating final values ..."))
+
input_maps = tmp_base + 'partial.%05d' % 1
for i in range(2, n+1):
input_maps += ',%spartial.%05d' % (tmp_base, i)
More information about the grass-commit
mailing list