[GRASS-user] building MVCs

[Glynn Clements]

Martin Brandt wrote:

> I have daily data and for each day a quality file.
> I'm trying to build 10-day MVCs (most valuable composite, taking the highest
> value) and the corresponding quality map from the daily raster images.
> 
> the MVCs work well with (example for the first one):
> r.series input="`g.mlist mapset=$maps pattern='NDVI_A'${maps}'00[1-9]'
> pattern='NDVI_A'${maps}'010' sep=,`" --o
> output=mvc_${maps}_01,max_${maps}_01 method=maximum,max_raster
> 
> but i'm struggling with the quality file. It should use the value of the
> raster which was used for the MVC. It would work with the following lines:
> 
> 
> r.mapcalc "qa_01 = if(max_'${maps}'_01 == 0, qamask_QA_A'${maps}'001, 0)"
> r.mapcalc "qa_02 = if(max_'${maps}'_01 == 1, qamask_QA_A'${maps}'002, 0)"
> r.mapcalc "qa_03 = if(max_'${maps}'_01 == 2, qamask_QA_A'${maps}'003, 0)"
> r.mapcalc "qa_04 = if(max_'${maps}'_01 == 3, qamask_QA_A'${maps}'004, 0)"
> r.mapcalc "qa_05 = if(max_'${maps}'_01 == 4, qamask_QA_A'${maps}'005, 0)"
> r.mapcalc "qa_06 = if(max_'${maps}'_01 == 5, qamask_QA_A'${maps}'006, 0)"
> r.mapcalc "qa_07 = if(max_'${maps}'_01 == 6, qamask_QA_A'${maps}'007, 0)"
> r.mapcalc "qa_08 = if(max_'${maps}'_01 == 7, qamask_QA_A'${maps}'008, 0)"
> r.mapcalc "qa_09 = if(max_'${maps}'_01 == 8, qamask_QA_A'${maps}'009, 0)"
> r.mapcalc "qa_10 = if(max_'${maps}'_01 == 9, qamask_QA_A'${maps}'010, 0)"

Note that you can generate multiple output maps from a single run of
r.mapcalc, e.g.:

	r.mapcalc <<EOF
	qa_01 = if(max_'${maps}'_01 == 0, qamask_QA_A'${maps}'001, 0)
	qa_02 = if(max_'${maps}'_01 == 1, qamask_QA_A'${maps}'002, 0)
	qa_03 = if(max_'${maps}'_01 == 2, qamask_QA_A'${maps}'003, 0)
	qa_04 = if(max_'${maps}'_01 == 3, qamask_QA_A'${maps}'004, 0)
	qa_05 = if(max_'${maps}'_01 == 4, qamask_QA_A'${maps}'005, 0)
	qa_06 = if(max_'${maps}'_01 == 5, qamask_QA_A'${maps}'006, 0)
	qa_07 = if(max_'${maps}'_01 == 6, qamask_QA_A'${maps}'007, 0)
	qa_08 = if(max_'${maps}'_01 == 7, qamask_QA_A'${maps}'008, 0)
	qa_09 = if(max_'${maps}'_01 == 8, qamask_QA_A'${maps}'009, 0)
	qa_10 = if(max_'${maps}'_01 == 9, qamask_QA_A'${maps}'010, 0)
	EOF

or:

	for i in 1 2 3 4 5 6 7 8 9 10 ; do
	    printf "qa_%02d = if(max_${maps}_01 == 0, qamask_QA_A${maps}%03d, 0)\n" $i $i
	done | r.mapcalc

Doing so will normally be more efficient that running r.mapcalc
multiple times.

> r.mapcalc "qa_mvc_'${maps}'_01 = qa_01 || qa_02 || qa_03 || qa_04 || qa_05
> || qa_06 || qa_07 || qa_08 || qa_09 || qa_10" 
> 
> 
> which look in the max_raster output which raster was used and use the
> corresponding quality file then, if the time series would be regular, but of
> course it's not. Sometimes days are missing, so max_raster does not point to
> the correct image any more.
> 
> Does anyone have an idea how this could be solved?

One option is, as Markus suggests, to create null maps (to save space,
these can be reclass maps), and insert them into the list of input
maps to r.series so that the time series is evenly spaced.

Another option is to note the relationship between the max_raster
values and input maps (i.e. store the output from g.mlist), and adjust
the r.mapcalc expression to use the appropriate map.

-- 
Glynn Clements <glynn at gclements.plus.com>