#!/bin/csh -f
########### Coodinate Asign
# confuse9の計算部分を修正した
#
#

set area_name="Paddy" 
set x0=1110
set y0=4920
set PIXEL=720
set LINE=620
@ x1 = $x0 + $PIXEL - 1
@ y1 = $y0 + $LINE - 1

########### TM image directory
#
#
#
set TM1=~/IWATE616/iwa_616y1    #
set TM2=~/IWATE616/iwa_616y2    #
set TM3=~/IWATE616/iwa_616y3    #
set TM4=~/IWATE616/iwa_616y4    #
set TM5=~/IWATE616/iwa_616y5    #
set TM6=~/IWATE616/iwa_616y7    #

########### DEM(Degital Elevetion Model) directory 
set DEM=~/DEM/iwa_utm              #
set SLP=~/DEM/iwa_slp              #
set ASP=~/DEM/iwa_asp              #

########### Vegetation Map directory
set VEG=~/VEG/iwaveg71

########### date of TM scene
#File    Time(GST)  North    East  Sun-angle Sun-zenith
#840426A 0:41:47    40:20   141:30   54       130.2
#840426B 0:42:11    38:54   141:03   54.7     128.2
#840604D 0:48:55    40:20   139:53   60.8     119.4
#840816A 0:44:25    40:20   141:17   53.6     128.3
#850616A 0:44:43    40:20   141:24   61.3     118
#850616B 0:45:07    38:54   140:58   61.7     115.2
#851022A 0:43:30    40:20   141:23   35       152.1
#851022B 0:43:54    38:54   140:56   35.4     150.5
#870521A 0:38:36    40:20   141:27   59.2     123.3
#870521B 0:39:00    38:53   140:59   59.8     120.8
#900427A 0:34:48    40:20   141:24   53.4     129
#900427B 0:35:12    38:54   140:57   54.1     127.1
#900614A 0:34:18    40:12   141:15   90       358
#931028A 0:36:45    40:19   141:26   32.5     151.5
#931028B 0:37:08    38:53   140:59   33.6     140.6
#940422A 0:34:35    40:20   141:22   51.4     129.1
#940524A 0:34:23    38:54   140:54   59.3     118.2
#950612A 0:20:58    40:20   141:22   57.9     112.4
#960529A 0:26:29    40:20   141:30   57.7     115.9
#960529B 0:26:53    38:54   141:03   58.1     113.4
#970613X 1:47:54    40:01   141:07   70       147
#970613Y 1:48:02    39:32   140:57   70       146
#
set mon=6
set date=16
set time=0.745


############ $area_name.doc log save
if ( -e $area_name.doc ) then      #
   rm $area_name.doc               #
endif                          #

############ 
echo Image Center Position Information > $area_name.doc  #タイトルを記録保存
echo "" >> $area_name.doc
echo '------------------------------------' >> $area_name.doc  #単なる区切り線
date >> $area_name.doc             #処理開始時間の保存
hostname >> $area_name.doc         #処理しているマシン名の保存
pwd >>  $area_name.doc             #現在位置（ディレクトリ）の記録
whoami >> $area_name.doc           #誰が作業したものか記録
echo '------------------------------------' >> $area_name.doc  #単なる区切り線

########### subset oparation

   if ( -e veg.img ) then
      rm veg.img               #もしveg.imgがあるなら、古いデータとみなし削除する
   endif
   if ( -e pre.bsq ) then
      rm pre.bsq               #もしpre.bsqがあるなら、古いデータとみなし削除する
   endif
   if ( -e pre.img ) then
      rm pre.img               #もしpre.imgがあるなら削除（通常は最初から無い）
   endif

########### Vegitation Image data load
subset $VEG.img veg.img $y0 $y1 $x0 $x1
paint $PIXEL $LINE veg.img ~/VEG/iwaveg7.tbl veg.ppm
cjpeg veg.ppm > veg.jpg
rm veg.ppm


########### TM data load
subset $TM1.img pre.img $y0 $y1 $x0 $x1   #バンド1画像切り出し
cat pre.img > pre.bsq ; rm pre.img        #切り出した画像でpre.bsq作成
subset $TM2.img pre.img $y0 $y1 $x0 $x1   #バンド2画像切り出し
cat pre.img >> pre.bsq ; rm pre.img       #切り出した画像でpre.bsq作成
subset $TM3.img pre.img $y0 $y1 $x0 $x1   #バンド3画像切り出し
cat pre.img >> pre.bsq ; rm pre.img       #切り出した画像でpre.bsq作成
subset $TM4.img pre.img $y0 $y1 $x0 $x1   #バンド4画像切り出し
cat pre.img >> pre.bsq ; rm pre.img       #切り出した画像でpre.bsq作成
subset $TM5.img pre.img $y0 $y1 $x0 $x1   #バンド5画像切り出し
cat pre.img >> pre.bsq ; rm pre.img       #切り出した画像でpre.bsq作成
subset $TM6.img pre.img $y0 $y1 $x0 $x1   #バンド7画像切り出し
cat pre.img >> pre.bsq ; rm pre.img       #切り出した画像でpre.bsq作成
############
#pre.bsqのdocファイル作成
rm pre.doc
echo 'title       : #'$area_name > pre.doc
echo 'data type   : byte' >> pre.doc
echo 'rows        : '$LINE >> pre.doc
echo 'columns     : '$PIXEL >> pre.doc
echo 'bands       : 5' >> pre.doc
echo 'class       : 0' >> pre.doc

############
xvp pre.bsq $PIXEL $LINE 2 1 0 3.5 4.5 1.8 pre321.ppm
xvp pre.bsq $PIXEL $LINE 3 4 5 2.0 2.0 2.0 pre456.ppm
cjpeg pre321.ppm > pre321.jpg
cjpeg pre456.ppm > pre456.jpg
rm pre321.ppm pre456.ppm


#############Iwate(4350x6800)l
#############Left/UP(x0,y0)=467400,4483500 Right/Down(x1,y1)=597900,4279500
#############140:37:30E,40:30:00N, 142:07:30E,38:40:00N
#画像の中心座標のUTM座標値計算　UTM;Universal Traverse Mercator Map 横メルカトル地図座標

@ cenx = ( $x0 + $x1 ) / 2      #中心の座標ｘ  
@ ceny = ( $y0 + $y1 ) / 2      #中心の座標ｙ
@ utmx = 467400 + 30 * $cenx    #ＵＴＭ座標値による値の算出
@ utmy = 4483500 - 30 *  $ceny  #ＵＴＭ座標値による値の算出

#############Data log od Center Scene
#画像の中心位置における情報の記録
echo Area       : $area_name >> $area_name.doc
echo x0,x1      : $x0,$x1 >> $area_name.doc
echo y0,y1      : $y0,$y1 >> $area_name.doc
echo Pixel,Line : $cenx,$ceny >> $area_name.doc
echo UTMX,Y     : $utmx,$utmy >> $area_name.doc
echo -n LON,LAT    :  >> $area_name.doc           #LON;経度、LAT;緯度
utm2lla $utmx $utmy >> $area_name.doc             #UTMから緯度経度へ変換する命令
echo '' >> $area_name.doc
echo '------------------------------------' >> $area_name.doc

############# 6S data set up
#6S（シックスエス）は、大気放射伝達シミュレータ。
#大気圏の条件（エアロゾルや標高値、緯度経度値、海上陸上など）
#でどのように地面の放射輝度が変化するか調査するプログラム
#
#
if ( -e pre.org ) then
   rm pre.org             #6Sに引き渡す計算条件をpre.orgに格納する
endif

echo   '7             (Landsat TM)' > pre.org                  #ランドサット画像の放射輝度計算
echo  -n $mon $date $time ' ' >> pre.org                       #月日と緯度経度を入れる
utm2lla $utmx $utmy >> pre.org                                 #緯度経度はutm2llaで計算（改行なし）
echo '  (for Solar Zenith Angle)' >> pre.org                   #
echo '2               (Mid Latitude Summer)' >> pre.org        #夏の季節
echo '3               (AEROSOLS MODEL)' >> pre.org             #海洋性エアロゾル
echo '50              (Direct Visibility Input km)' >> pre.org #見通し距離５０ｋｍ
echo '0.00            (TARGET ALTITUDE IN KM)' >> pre.org      #対象地域の標高値
echo '-1000           (SATELLITE CASE)' >> pre.org             #衛星の条件（不明）
echo '29              (Landsat TM Band 1)' >> pre.org          #ＴＭバンド画像１での結果を求める
echo '0               (HOMOGENEOUS CASE)' >> pre.org           #均質条件（不明）
echo '0               (NO DIRECTIONAL EFFECTS)' >> pre.org     #指向性効果（不明）
echo '1               (VEGETATION)' >> pre.org                 #地上被覆物を１（不明）
echo '-9              (ATMOSPHERIC CORRECTION FOR BRDF)' >> pre.org    #双方向反射BRDFの補正（不明）

############# 6s Calculate
if ( -e pre.ang ) then
   rm pre.ang              #6sの計算結果はpre.angに保存する
endif

#6sによるＴＭ画像撮影時の太陽高度角、方位角の計算実行。結果をpre.angへ。
#
cat pre.org  | sixs | grep 'solar zenith' > pre.ang

############# Sun angle,zenith value
#計算結果には不要な項も含まれており、削除や修正を行う部分
#awkは、BASICのような計算実行するUNIXコマンド
#suna;太陽入射角。6sの結果を９０度から引いた値になる。
#sunz;太陽方位角(zenith)。これは真北を０度とした値。
#
set suna = `awk '{ print 90-$5 }' < pre.ang`
set sunz = `awk '{ print $10 }' < pre.ang`

#計算結果を$area_name.docにも記録しておく
echo -n 'Sun angle,zenith :'  >> $area_name.doc
echo $suna, $sunz >> $area_name.doc
echo '------------------------------------' >> $area_name.doc


############# DEM data copy
#DEM(Digital Elevation Model)標高データ
#SLP(Slope)斜度、地形の傾きを示すデータ
#ASP(Aspect)方向、地形が真北なら０のデータ
#これらは、大気地形効果補正、太陽光の陰影除去のために使う
#
zsubset $DEM.img dem.img $y0 $y1 $x0 $x1
zsubset $SLP.img slp.img $y0 $y1 $x0 $x1
zsubset $ASP.img asp.img $y0 $y1 $x0 $x1
#gvc2 dem.img $PIXEL $LINE dem.pgm          #標高図をｘｖで確認観察する場合

############# Incident Data get
#Incident,入射。太陽光がどういう角度で入射しているかを示すデータ
#この処理の結果出力されるのが、「Ｃｏｓβ」値である。
#この計算のために、地形の傾き、方向、太陽の角度、方位情報が必要
#
zincident slp asp inc $suna $sunz
#gvc21 inc.img $PIXEL $LINE inc.pgm    #Ｃｏｓβの図（陰影の図）を確認したいとき


echo '------------------------------------' >> $area_name.doc
echo 'veg.img ' >> $area_name.doc
echo '  No.  Frequency' >> $area_name.doc

freq veg | awk 'NR>8 {print $0}' >> $area_name.doc

############ Correct operation
#実際の修正コサイン法による大気地形効果補正
#これで出力されるpost.bsqが、補正された結果であり、pre.bsqに比べ
#影の部分が平坦になったり、標高差による輝度の違いも少なくなっている
#

#繰り返し計算・大気地形効果補正による分類精度の変化


set cn1=0
while ( $cn1 < 6 )
set cn2=0
while ( $cn2 < 6 )
set cn3=0
while ( $cn3 < 6 )

set file_name="paddy$cn1$cn2$cn3"

@ k1 = 56 + $cn1 * 2
@ k2 = 18 + $cn2 * 2
@ k3 = 12 + $cn3 * 2

if ( -e pre.crt ) then
   rm pre.crt
endif

#echo '59.32  23.94  17.66  51.73  28.87  8.10'

echo -n $k1 " " $k2 " " $k3 " "> pre.crt
echo '51.73  28.87  8.10' >> pre.crt
echo '-0.003 -0.001  -0.001  0  0  0' >> pre.crt
echo '0.001 0.001 0.001 0.001 0.001 0.001' >> pre.crt
echo '0.6841 0.7755 0.8334 0.8860 0.9526 0.9680' >> pre.crt
echo '0.1392 0.0840 0.0897 0.0220 0.0016 0.0009' >> pre.crt

echo -n $k1 " " $k2 " " $k3 " " > $file_name.doc
echo ' 51.73  28.87  8.10' >> $file_name.doc

#echo 'Modified Cosine Coefficient Table' >> $file_name.doc
#cat pre.crt >> $file_name.doc
#echo '------------------------------------' >> $file_name.doc


cat pre.crt | zcorrect2 pre dem inc post

#xvp post.bsq $PIXEL $LINE 2 1 0 3.5 4.5 1.8 post321.ppm    #大気地形効果補正画像を確認したいとき
#xvp post.bsq $PIXEL $LINE 3 4 5 2.0 2.0 2.0 post456.ppm
#cjpeg post321.ppm > post321.jpg
#cjpeg post456.ppm > post456.jpg
#rm post321.ppm post456.ppm


############ Classification
# Classification,分類作業
#大気地形効果補正したＴＭ画像を使うと、分類結果がよいとされるので、その準備の
#ために、post.bsqを、分類用ファイル名a.bsqに変更する。
#もし、補正前の画像と比較したいときは、ここのpost部分をpreにすればよい。
#
mv post.bsq a.bsq    #ファイル名の変更。大気地形効果補正画像のBSQファイルを分類用に使う
mv post.doc a.doc    #付属ファイルも名前変更

############ 分類処理開始
#分類作業のパラメータ類を記録開始
#echo 'Classification Parameter' >> $file_name.doc
#echo '------------------------------------' >> $file_name.doc

############ ピラミッドリンキング
#分類結果が「ごま塩状」の画素にならないように、予め画像にフィルタをかけておく。
#
pcomp > out

############ ISODATA法によるクラスタ分類
#教師なし分類法の代表例であるISODATA法で、特徴ある画素グループを機械的にグループ化する。
#つまり、可視光線領域で輝度が高い（コンクリート等）とか、赤外で輝度が高い（葉緑素）、
#バンド５で暗くなる（水）などの特徴を、コンピュータの判断により勝手にグループ化していく。
#そのグループのことを、「クラスタ」といい、これが分類結果になる。
#
#isodataのパラメータ
#c部分;ピラミッドリンキングでフィルタがかかった画像を入力する、という意味
#g部分;code book fileの名前で、この後に付く番号は、分類処理が何回で収束したか、その回数値である。
#  code bookには、分類した各クラスの代表的な輝度値の一覧が記録される。これを見て、後でカテゴリー対応に利用する。
#16の部分;standard number of code bookで、クラスタ数を１６個に指定。実際には１６にならず、増減がある。
#10の部分;各クラスタが、値の上でどれだけ離れてもらいたいか、その距離。
#50の部分;各クラスタが最小でも何ピクセルのまとまりになってもらいたいか、その数。
#100の部分;minimum sample size for split
#18の部分;minimum distance between clusters for merging
#24の部分;maximum s.d. for spliting。s.dは標準偏差
#実際の処理は次のようにする
#
#isodata c g 16 10 50 100 18 24
isodata c g 30 10 50 100 18 24

#上の分類コマンドそのものを記録しておく
#たとえば16を20に変えたりすると、まったく別の結果が出てくる
#
echo 'isodata c g 30 10 50 100 18 24' >> $file_name.doc
echo "" >> $file_name.doc
echo '------------------------------------' >> $file_name.doc

#g53などの後の数値を無視して、全部g.codというファイル名にする
#
if (-e g.cod) then
 rm g.cod
endif

cp g*.cod g.cod


#分類作業の記録
#echo 'Unsupervised Clastering Result (g.rst) ' >> $file_name.doc
#echo '------------------------------------' >> $file_name.doc

#まず、code bookの内容を$file_name.docに書き込む
cat g.cod >> $file_name.doc


############# クラスタ分析にもとずいてクラスタによる画像作成
#まず、cy.imgを作成する
mlbg0x c g cy 100

#echo '------------------------------------' >> $file_name.doc
#echo 'Merging Parameter' >> $file_name.doc
#echo '------------------------------------' >> $file_name.doc
#echo 'mlbg0x c g cy 100' >> $file_name.doc

############# クラスタ分析にもとずいてクラスタによる画像作成・第２段階
#cy.imgに基づいてby.imgを作成する
prcont cy.img by.img b
############# クラスタ分析にもとずいてクラスタによる画像作成・第３段階
#by.imgに基づいてay.imgを作成する。ay.imgが、クラスタ分類による結果の画像
prcont by.img ay.img a
############# 分類結果であるay.imgを、モノクロで参照するための可視化
#gvcでもよいが、こういうコマンドもある。
#pgmform ay.img ay 1

############# 後の作業で使用しない画像データ等の削除
#
#rm a.bsq a.doc asp.doc asp.img dem.doc dem.img
#rm inc.doc inc.img pre.ang pre.bsq pre.crt pre.doc 
#rm pre.org slp.doc slp.img

rm a.lnk alap.bsq alap.cov alap.doc out
rm b.bsq b.doc b.lnk blap.bsq blap.cov blap.doc
rm by.doc by.img c.bsq c.doc cy.doc cy.img

cp g.cod g.rst
#rm *.cod

#echo '------------------------------------' >> $file_name.doc
#echo "" >> $file_name.doc
#echo -n 'Classification completed Time  :' >> $file_name.doc
#date >> $file_name.doc                #分類画像入手までの処理時間記録
#echo "" >> $file_name.doc
#echo '------------------------------------' >> $file_name.doc

############## カテゴリー対応作業
#echo ""  >> $file_name.doc
#echo "ay.img"   >> $file_name.doc
#echo ""  >> $file_name.doc

############## 植生図の凡例表示
#echo "" >> $file_name.doc
#echo '------------------------------------' >> $file_name.doc
#
 set num=`freq ay | tail -n 1 | awk '{print $1}'`

 set cattotal = 7

############## カテゴリー作業の開始
if (-e ay.tp1 ) then
 rm ay.tp1
endif
@ numx = $num + 2

compa $PIXEL $LINE 7 veg.img ay.img g.rst $numx 0 > ay.tp1

echo "compa" >> $file_name.doc
echo '------------------------------------' >> $file_name.doc
compa $PIXEL $LINE 7 veg.img ay.img g.rst $numx 1 >> $file_name.doc


############ カテゴリー対応結果に基づくリクラス処理
#リクラスのための対応表を自動作成する

if (-e ay.cls ) then
 rm ay.cls 
endif

echo -n "* new number of class : " > ay.cls
echo $cattotal >> ay.cls

cat ay.tp1 >> ay.cls

rm ay.tp1

cat ay.doc | awk 'NR<=5 {print $0}' > ay.tp1
echo -n 'class       : ' >> ay.tp1
@ num1 = $num + 1
echo $num1 >> ay.tp1
cp ay.tp1 ay.doc

########### リクラスの実行 ayr.imgが、分類の最終結果になる

echo '------------------------------------' >> $file_name.doc
echo "" >> $file_name.doc
echo '------------------------------------' >> $file_name.doc
cat ay.cls >>  $file_name.doc

reclass ay ay ayr

rm ay.tp1


########### カテゴリー対応による色割り当て結果
paint $PIXEL $LINE ayr.img ~/VEG/iwaveg7.tbl ayr.ppm
cjpeg ayr.ppm  > ayr$cn1$cn2$cn3.jpg
rm ayr.ppm

#"paddy$cn1$cn2"


########### Confuse 評価の結果を$file_name.docに記録
#echo '------------------------------------' >> $file_name.doc
#echo 'veg.img ' >> $file_name.doc
#echo '  No.  Frequency' >> $file_name.doc

#freq veg | awk 'NR>8 {print $0}' >> $file_name.doc

#echo '------------------------------------' >> $file_name.doc
#echo 'ay.img ' >> $file_name.doc
#echo '  No.  Frequency' >> $file_name.doc

#freq ay | awk 'NR>8 {print $0}' >> $file_name.doc

#echo '------------------------------------' >> $file_name.doc
#echo 'ayr.img ' >> $file_name.doc
#echo '  No.  Frequency' >> $file_name.doc

#freq ayr | awk 'NR>8 {print $0}' >> $file_name.doc

#echo '------------------------------------' >> $file_name.doc
#echo 'CONFUSE ' >> $file_name.doc
#echo '------------------------------------' >> $file_name.doc

echo "0:Unclass 1:Grass 2:Conifer 3:Broad 4:Paddy 5:Town 6:Water" >> $file_name.doc
#echo '------------------------------------' >> $file_name.doc
#confuse9 veg ayr >> $file_name.doc

confuse9 veg ayr | erchk 7 >> $file_name.doc

#exit

@ cn3 = $cn3 + 1
end
@ cn2 = $cn2 + 1
end
@ cn1 = $cn1 + 1
end

#echo '------------------------------------' >> $file_name.doc
#echo -n 'Proccessing completed time =  ' >> $file_name.doc
#date >> $file_name.doc                #全作業終了時間

# End of C-Shell Script
# Sato Kiyotada , Ichinoseki National College, 1999.1