Hi, here's the problem: I need to implement an operation in GRASS wich takes
an UNDEFINED number of raster layers and, for each position, takes the cell
of each raster, sort those values, makes some calculations and store the
result in the same position of the output raster.
I know how to do it with just two or three layers (simply, I use a r.mapcalc
sentence that takes care of the sorting issue with "IFs"), but now we don't
know how many inputs will be and even if we knew, with many layers the
mapcalc expressions would be inmense and awfully hard to generate. So, I
think there's no choice but to seek other alternatives, such as shell or
python scripting.
Now, the questions:
-Do you think that what I want to do is possible in GRASS?
-If so, how would approach this issue?
Thank you very much!
--
View this message in context: http://osgeo-org.1560.n6.nabble.com/How-to-implement-this-kind-of-operation-tp4984790.html
Sent from the Grass - Users mailing list archive at Nabble.com.
DavidRA wrote:
Hi, here’s the problem: I need to implement an operation in GRASS wich takes
an UNDEFINED number of raster layers and, for each position, takes the cell
of each raster, sort those values, makes some calculations and store the
result in the same position of the output raster.
Hey,
R.series does all the following calculations for any number of rasters at once:
average, count, median, mode, minimum, min_raster, maximum, max_raster, stddev, range, sum, variance, diversity, slope, offset, detcoeff, quart1, quart3, perc90, quantile, skewness, kurtosis
If you need more advanced calculations, I think you will need to use scripting. Or maybe you can tweak the code of r.series to put your own formulas.
Best,
Marcello.
I guess I could use the code of r.series, modifying it a bit to obtain a new
function. The problem is that... man, that code is hard to understand. I do
have some notions of C, but here I'm finding a lot of trouble to see what it
does. For example, I don't see where it does the calculations: there is a
section of the code with a comment saying "process the data", but inside of
it, where is the average calculated? The average, or any of the other
methods.
Thanks!
--
View this message in context: http://osgeo-org.1560.n6.nabble.com/How-to-implement-this-kind-of-operation-tp4984790p4985251.html
Sent from the Grass - Users mailing list archive at Nabble.com.
On Mon, Jul 2, 2012 at 10:56 AM, DavidRA <theboss777@gmail.com> wrote:
I guess I could use the code of r.series, modifying it a bit to obtain a new
function. The problem is that... man, that code is hard to understand. I do
have some notions of C, but here I'm finding a lot of trouble to see what it
does. For example, I don't see where it does the calculations: there is a
section of the code with a comment saying "process the data", but inside of
it, where is the average calculated? The average, or any of the other
methods.
What function exactly would you need? Instead of modifying r.series
yourself, you could also request an enhancement to r.series.
Markus M
It is a version of the ordered weighted averaging, with two kind of weights:
one associated with the raster layer and another one that depends on the
order of the cell. Look, this is how I do it with r.mapcalc when there are
just 2 raster layers:
where "w"i are the weights of each layer and "zi" the weights of order. The
process is simple: for each position, if any raster has NODATA, the output
raster will have a NODATA in that position. Else, if the cell of the first
raster has a bigger value than the second, we assign z1 and z2 and calculate
the average. If it is smaller, we assign z2 and z1 and do the calculation.
What I need to do now is the same, but with many rasters. I can't use
r.mapcacl bacause the number of "IF" sentences would be huge.
--
View this message in context: http://osgeo-org.1560.n6.nabble.com/How-to-implement-this-kind-of-operation-tp4984790p4985264.html
Sent from the Grass - Users mailing list archive at Nabble.com.
DavidRA wrote:
I guess I could use the code of r.series, modifying it a bit to obtain a new
function. The problem is that… man, that code is hard to understand.
Sorry if I gave you the impression that I am a skilled programmer in C.
Actually, I forgot the majority of GRASS modules are in C (or python now).
I always use shell scripts to do my stuff, so I was actually thinking of
Shell (what would be much easier to tweak).
I do
have some notions of C, but here I’m finding a lot of trouble to see what it
does. For example, I don’t see where it does the calculations: there is a
section of the code with a comment saying “process the data”, but inside of
it, where is the average calculated? The average, or any of the other
methods.
My God, it is indeed complex! I guess, somehow by including stats.h, they
are done by a simple call to the functions (somewhere, because I can’t see either),
so he is not implementing the calculations in the code.
Well, at least that is my most humble (and ignorant) opinion.
Cheers,
Marcello.
DavidRA wrote:
I guess I could use the code of r.series, modifying it a bit to obtain a new
function. The problem is that... man, that code is hard to understand. I do
have some notions of C, but here I'm finding a lot of trouble to see what it
does. For example, I don't see where it does the calculations: there is a
section of the code with a comment saying "process the data", but inside of
it, where is the average calculated? The average, or any of the other
methods.
The aggregate functions (c_ave, c_median, etc) are all in lib/stats.
The functions are also used by r.neighbors and r.resamp.stats.
r.series has a table of aggregates defined at the top of main.c. You
would need to define your aggregate function then add it to that
table.
--
Glynn Clements <glynn@gclements.plus.com>
Ok, with this I made some progress. I still don't understand many details of
the code, but the overall operation is now much clearer to me, so I'll keep
investigating on this way.
I'd also like to explore the scripting possibility. Do you think that mi
purpose could be achieved with shell or python scripts? Right now, these are
my concerns about it:
- How to handle a undefined number of input raster layers.
- How to access the pixels of each one in order to make the calculations.
- How to create the output raster, pixel by pixel, with the results of those
calculations.
Thanks all of you for your help!
--
View this message in context: http://osgeo-org.1560.n6.nabble.com/How-to-implement-this-kind-of-operation-tp4984790p4985722.html
Sent from the Grass - Users mailing list archive at Nabble.com.
DavidRA wrote:
I'd also like to explore the scripting possibility. Do you think that mi
purpose could be achieved with shell or python scripts?
It's not really feasible to do it with a "script". You can do it with
a Python program which uses the GRASS raster API; all of the core
GRASS libraries have Python bindings, so whatever you can do in C you
can do in Python. But doing it in Python would be just as much work as
doing it in C, so you're basically talking about re-writing r.series
in Python.
Extending r.series should be fairly straightforward. Look at the
existing aggregates in lib/stats/*.c for reference (you don't need a
w_* function; those aren't used by r.series). Your function can be
added to either libstats or to r.series. Once you've created a
function, add an entry in the "menu" array at the top of
r.series/main.c.
--
Glynn Clements <glynn@gclements.plus.com>
I see, then I'll leave aside the scripting option in order to delve into the
r.series code.
Thanks a lot!
--
View this message in context: http://osgeo-org.1560.n6.nabble.com/How-to-implement-this-kind-of-operation-tp4984790p4985959.html
Sent from the Grass - Users mailing list archive at Nabble.com.