On Sat, Jul 9, 2011 at 10:24 PM, Martin Brandt <martin.bra...@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 _______________________________________________ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user