[GRASS-user] building MVCs

[Markus Metz]

On Sat, Jul 9, 2011 at 10:24 PM, Martin Brandt
<martin.brandt at univie.ac.at> wrote:
> dear all,
>
> 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)"
>
> 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?
>
You need to create maps containing only NULL values for the missing
days, then the max_raster output of r.series should point to the
correct map number:

firstday=1
# lastday=366 for leap years
lastday=365
for map_no in `seq $firstday $lastday  | awk '{printf "%03d\n", $1}' ` ; do
    map='NDVI_A${maps}${map_no}'
    eval `g.findfile file=$map element=cell mapset=.`
    if [[ $name = "" ]] ; then
      r.mapcalc "$map = null()"
    fi
done

Now as before:

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

HTH,

Markus M