[GRASS-dev] grass7 - python script - moving window: how to simplify/accelerate r.mapcalc?
Helmut Kudrnovsky
hellik at web.de
Sat May 31 06:58:51 PDT 2014
>The most straightforward change is to eliminate the if() functions.
>Any boolean value (e.g. the result of a comparison) is an integer,
>with 1 for true, 0 for false. So "if(a<b,1,0)" is equivalent to just
>"(a<b)". This should also provide a minor increase in efficiency.
ah thanks. a quick diff between r.mapcalc with and without the if()
function:
| Range of data: min = 0 max = 0
| Data Description:
|
| generated by r.mapcalc
|
|
|
| Comments:
|
| elevation_percentile_step2 at user1 - elevation_percentile_step2a at user1
|
yes, it works.
now a test in the python shell of the wxgui with nc sample dataset:
>>>import grass.script as grass
>>>offsets = [d
for j in xrange(1,6+1)
for i in [j,-j]
for d in [(i,0),(0,i),(i,i),(i,-i)]]
>>>offsets
[(1, 0), (0, 1), (1, 1), (1, -1), (-1, 0), (0, -1), (-1, -1), (-1, 1), (2,
0), (0, 2), (2, 2), (2, -2), (-2, 0), (0, -2), (-2, -2), (-2, 2), (3, 0),
(0, 3), (3, 3), (3, -3), (-3, 0), (0, -3), (-3, -3), (-3, 3), (4, 0), (0,
4), (4, 4), (4, -4), (-4, 0), (0, -4), (-4, -4), (-4, 4), (5, 0), (0, 5),
(5, 5), (5, -5), (-5, 0), (0, -5), (-5, -5), (-5, 5), (6, 0), (0, 6), (6,
6), (6, -6), (-6, 0), (0, -6), (-6, -6), (-6, 6)]
>>>terms = ["(myelevnc[%d,%d] < myelevnc)" % d
for d in offsets]
>>>terms
['(myelevnc[1,0] < myelevnc)', '(myelevnc[0,1] < myelevnc)', '(myelevnc[1,1]
< myelevnc)', '(myelevnc[1,-1] < myelevnc)', '(myelevnc[-1,0] < myelevnc)',
'(myelevnc[0,-1] < myelevnc)', '(myelevnc[-1,-1] < myelevnc)',
'(myelevnc[-1,1] < myelevnc)', '(myelevnc[2,0] < myelevnc)', '(myelevnc[0,2]
< myelevnc)', '(myelevnc[2,2] < myelevnc)', '(myelevnc[2,-2] < myelevnc)',
'(myelevnc[-2,0] < myelevnc)', '(myelevnc[0,-2] < myelevnc)',
'(myelevnc[-2,-2] < myelevnc)', '(myelevnc[-2,2] < myelevnc)',
'(myelevnc[3,0] < myelevnc)', '(myelevnc[0,3] < myelevnc)', '(myelevnc[3,3]
< myelevnc)', '(myelevnc[3,-3] < myelevnc)', '(myelevnc[-3,0] < myelevnc)',
'(myelevnc[0,-3] < myelevnc)', '(myelevnc[-3,-3] < myelevnc)',
'(myelevnc[-3,3] < myelevnc)', '(myelevnc[4,0] < myelevnc)', '(myelevnc[0,4]
< myelevnc)', '(myelevnc[4,4] < myelevnc)', '(myelevnc[4,-4] < myelevnc)',
'(myelevnc[-4,0] < myelevnc)', '(myelevnc[0,-4] < myelevnc)',
'(myelevnc[-4,-4] < myelevnc)', '(myelevnc[-4,4] < myelevnc)',
'(myelevnc[5,0] < myelevnc)', '(myelevnc[0,5] < myelevnc)', '(myelevnc[5,5]
< myelevnc)', '(myelevnc[5,-5] < myelevnc)', '(myelevnc[-5,0] < myelevnc)',
'(myelevnc[0,-5] < myelevnc)', '(myelevnc[-5,-5] < myelevnc)',
'(myelevnc[-5,5] < myelevnc)', '(myelevnc[6,0] < myelevnc)', '(myelevnc[0,6]
< myelevnc)', '(myelevnc[6,6] < myelevnc)', '(myelevnc[6,-6] < myelevnc)',
'(myelevnc[-6,0] < myelevnc)', '(myelevnc[0,-6] < myelevnc)',
'(myelevnc[-6,-6] < myelevnc)', '(myelevnc[-6,6] < myelevnc)']
>>>expr = "elevation_percentile = (100.0 / 48.0) * \\\n(%s)" % " +
\\\n".join(terms)
>>>expr
elevation_percentile = (100.0 / 48.0) * \
((myelevnc[1,0] < myelevnc) + \
(myelevnc[0,1] < myelevnc) + \
(myelevnc[1,1] < myelevnc) + \
(myelevnc[1,-1] < myelevnc) + \
(myelevnc[-1,0] < myelevnc) + \
(myelevnc[0,-1] < myelevnc) + \
(myelevnc[-1,-1] < myelevnc) + \
(myelevnc[-1,1] < myelevnc) + \
(myelevnc[2,0] < myelevnc) + \
(myelevnc[0,2] < myelevnc) + \
(myelevnc[2,2] < myelevnc) + \
(myelevnc[2,-2] < myelevnc) + \
(myelevnc[-2,0] < myelevnc) + \
(myelevnc[0,-2] < myelevnc) + \
(myelevnc[-2,-2] < myelevnc) + \
(myelevnc[-2,2] < myelevnc) + \
(myelevnc[3,0] < myelevnc) + \
(myelevnc[0,3] < myelevnc) + \
(myelevnc[3,3] < myelevnc) + \
(myelevnc[3,-3] < myelevnc) + \
(myelevnc[-3,0] < myelevnc) + \
(myelevnc[0,-3] < myelevnc) + \
(myelevnc[-3,-3] < myelevnc) + \
(myelevnc[-3,3] < myelevnc) + \
(myelevnc[4,0] < myelevnc) + \
(myelevnc[0,4] < myelevnc) + \
(myelevnc[4,4] < myelevnc) + \
(myelevnc[4,-4] < myelevnc) + \
(myelevnc[-4,0] < myelevnc) + \
(myelevnc[0,-4] < myelevnc) + \
(myelevnc[-4,-4] < myelevnc) + \
(myelevnc[-4,4] < myelevnc) + \
(myelevnc[5,0] < myelevnc) + \
(myelevnc[0,5] < myelevnc) + \
(myelevnc[5,5] < myelevnc) + \
(myelevnc[5,-5] < myelevnc) + \
(myelevnc[-5,0] < myelevnc) + \
(myelevnc[0,-5] < myelevnc) + \
(myelevnc[-5,-5] < myelevnc) + \
(myelevnc[-5,5] < myelevnc) + \
(myelevnc[6,0] < myelevnc) + \
(myelevnc[0,6] < myelevnc) + \
(myelevnc[6,6] < myelevnc) + \
(myelevnc[6,-6] < myelevnc) + \
(myelevnc[-6,0] < myelevnc) + \
(myelevnc[0,-6] < myelevnc) + \
(myelevnc[-6,-6] < myelevnc) + \
(myelevnc[-6,6] < myelevnc))
the expression expr is build correct, thanks.
but then while do
>>>grass.mapcalc("%s" % expr)
or
>>>grass.mapcalc( expr )
an error raise up with following message with DEBUG=3
GRASS 7.1.svn> D2/3: G_file_name(): path =
C:\grassdata/nc_spm_08_grass7/user1
D1/3: G_find_raster2(): name=\ mapset=
D2/3: G_file_name(): path = C:\grassdata/nc_spm_08_grass7/user1/SEARCH_PATH
D2/3: G_file_name(): path = C:\grassdata/nc_spm_08_grass7/PERMANENT
D2/3: G_file_name(): path = C:\grassdata/nc_spm_08_grass7/user1/cell/\
D2/3: G_file_name(): path = C:\grassdata/nc_spm_08_grass7/PERMANENT/cell/\
GRASS_INFO_WARNING(8628,1): 'cell/\' wurde in mehreren Mapsets gefunden
(auch in
<PERMANENT>).
GRASS_INFO_END(8628,1)
GRASS_INFO_WARNING(8628,2): Verwende <\@user1>
GRASS_INFO_END(8628,2)
D1/3: G_find_raster2(): name=\ mapset=user1
D2/3: G_file_name(): path = C:\grassdata/nc_spm_08_grass7/user1/cell/\
D2/3: G_file_name(): path = C:\grassdata/nc_spm_08_grass7/user1/fcell/\
D1/3: G_find_raster2(): name=\ mapset=user1
D2/3: G_file_name(): path = C:\grassdata/nc_spm_08_grass7/user1/cell/\
GRASS_INFO_ERROR(8628,3): Kann
'C:\grassdata/nc_spm_08_grass7/user1/cell_misc/\/
f_format' nicht finden.
GRASS_INFO_END(8628,3)
D2/3: G_file_name(): path = C:\grassdata/nc_spm_08_grass7/user1
GRASS_INFO_ERROR(5412,1): An error occurred while running r.mapcalc
GRASS_INFO_END(5412,1)
? ... not a very informative error message
just a another quick test
mapinput = "myelevnc"
mapinput
myelevnc
mapoutput = "myelevnc_1"
mapoutput
myelevnc_1
expression = "%s = %s + 1.0" % (mapoutput, mapinput)
expression
myelevnc_1 = myelevnc + 1.0
grass.mapcalc( expression )
grass.mapcalc( expression )
None
grass.read_command("r.info", map = mapoutput)
+----------------------------------------------------------------------------+
| Layer: myelevnc_1 Date: Sat May 31 15:20:02 2014
|
| Mapset: user1 Login of Creator: myricaria
|
| Location: nc_spm_08_grass7
|
| DataBase: C:\grassdata
|
| Title: ( myelevnc_1 )
|
| Timestamp: none
|
|----------------------------------------------------------------------------|
|
|
| Type of Map: raster Number of Categories: 0
|
| Data Type: DCELL
|
| Rows: 1350
|
| Columns: 1500
|
| Total Cells: 2025000
|
| Projection: Lambert Conformal Conic
|
| N: 228500 S: 215000 Res: 10
|
| E: 645000 W: 630000 Res: 10
|
| Range of data: min = 56.5787925720215 max = 157.329864501953
|
|
|
| Data Description:
|
| generated by r.mapcalc
|
|
|
| Comments:
|
| myelevnc + 1
|
|
|
+----------------------------------------------------------------------------+
this works.
any idea?
>But if you want any significant gain in efficiency, the obvious
>solution is to modify r.neighbors to add a new aggregate (similar to
>interspersion, but using less-than rather than not-equal). The weight=
>option can be used to specify the window shape (for an aggregate which
>lacks a weighted version, the weights are converted to a mask, where
>cells with a non-zero weight are used and those with a zero weight
>ignored).
[...]
>This type of problem comes up often enough that we really want an
>aggregate (method) which evaluates an arbitrary expression with the
>cell values and weights as inputs.
that would be nice.
thanks.
-----
best regards
Helmut
--
View this message in context: http://osgeo-org.1560.x6.nabble.com/grass7-python-script-moving-window-how-to-simplify-accelerate-r-mapcalc-tp5143161p5143415.html
Sent from the Grass - Dev mailing list archive at Nabble.com.
More information about the grass-dev
mailing list