[GRASS-user] building MVCs

GRASS GRASS user
1

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?

-----
Department of Geography and
Regional Research
UZA II, Althanstr. 14
1090 Wien, AUSTRIA
http://geooekologie.univie.ac.at
--
View this message in context: http://osgeo-org.1803224.n2.nabble.com/building-MVCs-tp6566371p6566371.html
Sent from the Grass - Users mailing list archive at Nabble.com.

2

On Sat, Jul 9, 2011 at 10:24 PM, Martin Brandt
<martin.brandt@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

3

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@gclements.plus.com>

4

brilliant ideas, thanks a lot!

-----
Department of Geography and
Regional Research
UZA II, Althanstr. 14
1090 Wien, AUSTRIA
http://geooekologie.univie.ac.at
--
View this message in context: http://osgeo-org.1803224.n2.nabble.com/building-MVCs-tp6566371p6567459.html
Sent from the Grass - Users mailing list archive at Nabble.com.