i.oifcalc 1.8 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364
  1. #!/bin/sh
  2. # Parameters: 1, 2, 3 - number of LANDSAT band
  3. # written by Markus Neteler 21. Jul. 1998, 5. Jan. 1999
  4. # neteler geog.uni-hannover.de | neteler itc.it
  5. # submodule of i.oif
  6. if [ $# -ne 3 ] ; then
  7. g.message -e "You must use i.oif instead of i.oifcalc (this is a submodule)"
  8. exit 1
  9. fi
  10. if [ $1 -lt 1 ] || [ $1 -gt 7 ] || [ $2 -lt 1 ] || [ $2 -gt 7 ] || [ $3 -lt 1 ] || [ $3 -gt 7 ] ; then
  11. g.message -e "Invalid channel combination ($1 $2 $3)"
  12. exit 1
  13. fi
  14. PROG=`basename $0`
  15. #### check if we have awk
  16. if [ ! -x "`which awk`" ] ; then
  17. g.message -e "awk required, please install awk or gawk first"
  18. exit 1
  19. fi
  20. # setting environment, so that awk works properly in all languages
  21. unset LC_ALL
  22. LC_NUMERIC=C
  23. export LC_NUMERIC
  24. # calculate SUM of Stddeviations:
  25. "$GISBASE"/etc/i.oif/m.cutmatrix "$temp_stddev" 1 $1 > "$temp_sum"
  26. "$GISBASE"/etc/i.oif/m.cutmatrix "$temp_stddev" 1 $2 >> "$temp_sum"
  27. "$GISBASE"/etc/i.oif/m.cutmatrix "$temp_stddev" 1 $3 >> "$temp_sum"
  28. cat "$temp_sum" | awk -v temp_file="$temp_calc" \
  29. 'BEGIN {sum = 0.0}
  30. NR == 1{}
  31. {sum += $1 ; N++}
  32. END{
  33. print sum > temp_file
  34. }'
  35. # calculate SUM of absolute(Correlation values):
  36. "$GISBASE"/etc/i.oif/m.cutmatrix "$temp_correlation" $1 $2 > "$temp_sum"
  37. "$GISBASE"/etc/i.oif/m.cutmatrix "$temp_correlation" $1 $3 >> "$temp_sum"
  38. "$GISBASE"/etc/i.oif/m.cutmatrix "$temp_correlation" $2 $3 >> "$temp_sum"
  39. cat "$temp_sum" | awk -v temp_file="$temp_calc" \
  40. 'BEGIN {sum = 0.0}
  41. NR == 1{}
  42. {sum += sqrt($1*$1) ; N++} # sqrt(x^2) is my ABS-function
  43. END{
  44. print sum >> temp_file
  45. }'
  46. # Calculate OIF index:
  47. # Divide (SUM of Stddeviations) and (SUM of Correlation)
  48. cat "$temp_calc" | awk 'BEGIN {}
  49. {count = divisor; divisor = $1;} # value shift - no better idea ;-)
  50. END{
  51. print count / divisor
  52. }'