I'm not sure that you found an answer to this yet, Ahmet - if you did,
please ignore this suggestion:
1) use r.mapcalc to create an integer combined index for the combinations
of your raster layers:
tmp1 = mapa*10 + mapb
tmp2 = tmp1*10 + mapc
....
or some such way of generating an integer index differentiating the
different types. If there are raster values outside 0-9, maybe reclass
first, or use 100 instead of 10 to left shift the index where there are
10-99 discrete values.
2) In R, generate the same integer index, and make a data frame with the
integer index in one column, and the probabilities in another.
3) Use the R/GRASS interface (or r.stats -1 and scan()) to move either the
whole of your index raster, or more likely tiles made by changing the
GRASS region, to R.
4) in R, use vectorised lookup to search for the probability value given
the index value - match() is fast and powerful (I haven't checked this,
but do use match() in other code - it will probably work much better on
tiles than your big raster) and assign the probability values to a new
vector.
5) move this vector back to GRASS as FCELL or DCELL. If you are tiling,
finish off by patching the tiles back together.
6) My guess is that the tiling, import, matching, and export can all be
done in a loop within R.
Rationale: the combined index avoids all the if statements - 11K if's is
rather a lot! I did think about pushing it all out say to postgresql and
using SQL (one table of RHS raster values, one table of lookups) but there
were still seven fields on each side to index on, so it seemed at least
worth trying to combine them automatically by left shifting (probably not
the right word) and addition, which are cheap operations in r.mapcalc.
Please let me know how you tackled this problem, your solution will have
general interest, and should be in the general GRASS knowledge base. I can
expand on this off-list if it is of any use.
Roger
On Wed, 11 Jun 2003, Glynn Clements wrote:
orkun wrote:
> >>I have a huge matrix containing 7 column and 10500 rows.
> >>Each columns stand for different maps. Each rows represent "unique"
> >>combinations of catagories of maps and probability values
> >>of each row.
> >>Something like:
> >>
> >> map1 map2 map3 ... prob
> >>1 cat1 x x x x
> >>2 cat2 x x x x
> >>. .
> >>. .
> >>10500 .
> >>
> >>So I have 10500 if-clause to create a probability map.
> >>I tried r.mapcalc (as from a text file) but grass gave "stack overflow"
> >>error.
> >>
> >>
> >
> >I'm not surprised.
> >
> >
> >
> >>Is there anyone making suggestions ?
> >>
> >>
> >
> >You need to describe your problem in more detail; it isn't clear what
> >you are trying to calculate.
> >
> >
> >
> Thank you for your interest
>
> You wanted to know the problem in more detail:
>
> what I did:
> 1- using r.stat command I produced a text file for statistical calculations
> 2- using R statistical package I produced probability results of each
> map's for each categories.
> it means 7 columns of maps 10500 rows plus probability of each
> combination. Actually it is a contingency table. Any data in table
> (matrix) is resulting from unique combination of maps and their
> categories.
> 3- I want to port these results to grass to create a map (of landslide
> susceptibility map) . Firstly I thought I used r.mapcalc with piping
> from text file. But because text file contains a huge matrix that is
> made up of 10500 if statement (for same number of rows) process failed.
> I attached a part of text file after modifying for r.mapcalc. you see
> all ifs are unique.
OK, so: each line is a rule; for each pixel you scan the list of rules
to find the first match, then output the associated result. Right?
The most straightforward approach (although not particularly
efficient), would be to use an iterative approach, i.e.
n=0
last='null()'
while read a1 a2 a3 a4 a5 a6 a7 a8 ; do
result=pr2.pass$n
r.mapcalc "$result = if (f28geo5==$a1 && f28slpcat==$a2 && f28concavity==$a3 && n_facing==$a4 && road_bas==$a5 &&str_buf==$a6 && f28f4==$a7, $a8, $last)"
last=$result
done < rules.txt
where rules.txt looks like:
21 6 1 1 1 1 0 701.22
7 4 2 1 1 2 1 702.06
...
Alternatively, you could combine this with your initial approach, only
matching a limited number of rules on each pass, and using several
passes.
However, I note that there seem to be very few possible values for
most of the maps. In that case, a far more efficient approach would be
to:
1. Reclass each input map so that all of the significant categories
are consecutively numbered from zero (0,1,2,3,...); any insignificant
categories would be merged into the last category.
2. Compute a composite map, using e.g.
composite = f28geo5 + 22 * (f28slpcat + 9 * (f28concavity + 3 * ...
where the numbers are the total number of categories in the preceding
map.
3. Assign the values as labels for the composite categories.
4. Use "@map" in r.mapcalc to use the label instead of the category.
--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Breiviksveien 40, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 93 93
e-mail: Roger.Bivand@nhh.no