4

I need calculated monthly mean of Snowcover data by MODIS product, and then modify this script gdal_calc.py in this section

# create alphabetic list for storing input layers
AlphaList=["A","B","C","D","E","F","G","H","I","J","K","L","M",
       "N","O","P","Q","R","S","T","U","V","W","X","Y","Z","a","b","c","d","e"]

and wrote this script:

!/bin/bash echo -n "./gdal_calc.py " > comandos_py.sh array=(A B C D E F G H I J K L M N O P Q R S T U V W X Y Z a b c d e) cont=0 for file in MCD10A1_*.tif ; do echo "map: $file" letra=${array[$cont]} echo -n "-$letra ${file} " >> comandos_py.sh cont=$[cont+1] done echo -n "--outfile=MEAN_2013_06.tif --calc=\"(A+B+C+D+E+F+G+H+I+J+K+M+N+O+P+Q+R+S+T+U+V+W+X+Y+Z+a+b+c+d)/30\" --NoDataValue=-9999 --type=Float32" >> comandos_py.sh

both running well..... but the output file is incorrect...

in another way! load all file on QGIS 2.6, open RASTER CALCULATOR insert 30 image SNOWCOVER DAILY sum and divided by 30 and the output file IS CORRECT....

can someone explain to me why this happens?

2

Because this array (input raster) includes 31 values:

array=(A B C D E F G H I J K L M  N O P Q R S T U V W X Y Z a b c d e)

and this one (output raster) only includes 29; where raster "L" and raster "e" are missing:

(A+B+C+D+E+F+G+H+I+J+K+M+N+O+P+Q+R+S+T+U+V+W+X+Y+Z+a+b+c+d)/30

However, the script would be (you have to repeat the letters):

!/bin/bash
echo -n "gdal_calc.py " > comandos_py.sh
array=(A B C D E F G H I J K L M N O P Q R S T U V W X Y Z A B C D)
cont=0
for file in MCD10A1_*.tif ; do
  echo "map: $file"
  letra=${array[$cont]}
  echo -n "-$letra ${file} " >> comandos_py.sh
  cont=$[cont+1]
done
echo -n "--outfile=MEAN_2013_06.tif --calc=\"(A+B+C+D+E+F+G+H+I+J+K+L+M+N+O+P+Q+R+S+T+U+V+W+X+Y+Z+A+B+C+D)/30\" --NoDataValue=-9999 --type=Float32" >> comandos_py.sh 

I tested it with 30 raster and it works! The first pixel of MEAN_2013_6.tif is the mean of the 30 raster; as it was manually corroborated with the gnome calculator.

enter image description here

Editing Note:

Running my script version:

zeito@debian:~/Desktop/raster$ ./script
map: MCD10A1_10.tif
map: MCD10A1_11.tif
map: MCD10A1_12.tif
map: MCD10A1_13.tif
map: MCD10A1_14.tif
map: MCD10A1_15.tif
map: MCD10A1_16.tif
map: MCD10A1_17.tif
map: MCD10A1_18.tif
map: MCD10A1_19.tif
map: MCD10A1_1.tif
map: MCD10A1_20.tif
map: MCD10A1_21.tif
map: MCD10A1_22.tif
map: MCD10A1_23.tif
map: MCD10A1_24.tif
map: MCD10A1_25.tif
map: MCD10A1_26.tif
map: MCD10A1_27.tif
map: MCD10A1_28.tif
map: MCD10A1_29.tif
map: MCD10A1_2.tif
map: MCD10A1_30.tif
map: MCD10A1_3.tif
map: MCD10A1_4.tif
map: MCD10A1_5.tif
map: MCD10A1_6.tif
map: MCD10A1_7.tif
map: MCD10A1_8.tif
map: MCD10A1_9.tif

No errors!

Changing running mode in comandos_py and running it.

chmod +x comandos_py.sh
./comandos_py.sh
0 .. 100 - Done

No errors again!

1
  • I dont understand how this works. I just tried gdal_calc.py -A a1.tif -A r2.tif --outfile=res.tif --calc="(A + A)/2" and res.tif just has the same values as r2.tif. How can we give two different rasters the same argument and expect gdal to know the difference? Jul 15 at 21:22
0

Thanks a lot for you Answer, but I think so the problem ..... look at this!!!

david@XXXX:/media/fuse/dat3/modis/mod10a1/raw/monthly/regular/stack/dayly/2013_06$ ./crea_comandos_py.sh map: MCD10A1_2013152.tif map: MCD10A1_2013153.tif map: MCD10A1_2013154.tif map: MCD10A1_2013155.tif map: MCD10A1_2013156.tif map: MCD10A1_2013157.tif map: MCD10A1_2013158.tif map: MCD10A1_2013159.tif map: MCD10A1_2013160.tif map: MCD10A1_2013161.tif map: MCD10A1_2013162.tif map: MCD10A1_2013163.tif map: MCD10A1_2013164.tif map: MCD10A1_2013165.tif map: MCD10A1_2013166.tif map: MCD10A1_2013167.tif map: MCD10A1_2013168.tif map: MCD10A1_2013169.tif map: MCD10A1_2013170.tif map: MCD10A1_2013171.tif map: MCD10A1_2013172.tif map: MCD10A1_2013173.tif map: MCD10A1_2013174.tif map: MCD10A1_2013175.tif map: MCD10A1_2013176.tif map: MCD10A1_2013177.tif map: MCD10A1_2013178.tif map: MCD10A1_2013179.tif map: MCD10A1_2013180.tif map: MCD10A1_2013181.tif david@XXXX:/media/fuse/dat3/modis/mod10a1/raw/monthly/regular/stack/dayly/2013_06$ sh comandos_py.sh Traceback (most recent call last): File "./gdal_calc.py", line 296, in main() File "./gdal_calc.py", line 275, in main eval('parser.add_option("-%s", dest="%s", help="input gdal raster file, note you can use any letter A-Z")' %(myAlpha, myAlpha)) File "", line 1, in File "/usr/lib/python2.7/optparse.py", line 1020, in add_option self._check_conflict(option) File "/usr/lib/python2.7/optparse.py", line 995, in _check_conflict option) optparse.OptionConflictError: option -A: conflicting option string(s): -A david@XXXX:/media/fuse/dat3/modis/mod10a1/raw/monthly/regular/stack/dayly/2013_06$ cat comandos_py.sh ./gdal_calc.py -A MCD10A1_2013152.tif -B MCD10A1_2013153.tif -C MCD10A1_2013154.tif -D MCD10A1_2013155.tif -E MCD10A1_2013156.tif -F MCD10A1_2013157.tif -G MCD10A1_2013158.tif -H MCD10A1_2013159.tif -I MCD10A1_2013160.tif -J MCD10A1_2013161.tif -K MCD10A1_2013162.tif -L MCD10A1_2013163.tif -M MCD10A1_2013164.tif -N MCD10A1_2013165.tif -O MCD10A1_2013166.tif -P MCD10A1_2013167.tif -Q MCD10A1_2013168.tif -R MCD10A1_2013169.tif -S MCD10A1_2013170.tif -T MCD10A1_2013171.tif -U MCD10A1_2013172.tif -V MCD10A1_2013173.tif -W MCD10A1_2013174.tif -X MCD10A1_2013175.tif -Y MCD10A1_2013176.tif -Z MCD10A1_2013177.tif -A MCD10A1_2013178.tif -B MCD10A1_2013179.tif -C MCD10A1_2013180.tif -D MCD10A1_2013181.tif --outfile=MEAN_2013_06.tif --calc="(A+B+C+D+E+F+G+H+I+J+K+L+M+N+O+P+Q+R+S+T+U+V+W+X+Y+Z+A+B+C+D)/30" --NoDataValue=-9999 --type=Float32david@gazp8:/media/fuse/dat3/modis/mod10a1/raw/monthly/regular/stack/dayly/2013_06$ ..

It so strange!!! Raster Calculator running ok enter image description here enter image description here

Thanks a lot!!

1
  • Please, edit your answer with formated code because is unclear. I edited my answer with my results (running my own version of your script ). I didn't get any errors or strange results.
    – xunilk
    May 19 '15 at 6:16
-1

xunilk's answer didnt work for me, so I had to implement a very hacky solution and make edits to gdal_calc.py directly. But the edits were simple.

Just make the variable AlphaList on line 57 to be:

AlphaList=["A","B","C","D","E","F","G","H","I","J","K","L","M",
           "N","O","P","Q","R","S","T","U","V","W","X","Y","Z",
           "a","b","c","d","e","f","g","h","i","j","k","l","m",
           "n","o","p","q","r","s","t","u","v","w","x","y","z"]

And then the script will be able to take tiffs marked with lowercase and accept up to 52 variables.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.