#!/bin/csh -f
#.cshrcを経由しないで立ちあげる方法(cshrcとは環境設定ファイルのこと）
#
#ＴＭ画像による土地被覆分類の処理用　シェルスクリプト    bysugupon 
#1999年01月06日更新
#左端に#がつくと、その行は実行されない。目的により、一部分の実行を行い
#その場合、#で飛び越えさせる。if then endifでもよい。
#
########### TM image directory
#ＴＭ画像ファイル名を絶対パスで表示し、img拡張子を省略した。
#バンド６は使わなかった。（理由は遠赤外で分解能が200mと荒いため。）
#

set TM1=/usr1/DATA/IMAGE/GIS/IWATE850616/iwa_616y1    #バンド1のファイル名
set TM2=/usr1/DATA/IMAGE/GIS/IWATE850616/iwa_616y2    #バンド2のファイル名
set TM3=/usr1/DATA/IMAGE/GIS/IWATE850616/iwa_616y3    #バンド3のファイル名
set TM4=/usr1/DATA/IMAGE/GIS/IWATE850616/iwa_616y4    #バンド4のファイル名
set TM5=/usr1/DATA/IMAGE/GIS/IWATE850616/iwa_616y5    #バンド5のファイル名
set TM6=/usr1/DATA/IMAGE/GIS/IWATE850616/iwa_616y7    #バンド7のファイル名



########### DEM(Degital Elevetion Model) directory 
#ＤＥＭ（デジタル標高データ）のファイル名。絶対パス表示、img拡張子省略
#utmが、30mx30mメッシュ間隔の精度0.1mの標高データ（２バイトデータ）
#slpが、30mx30mメッシュ間隔の傾斜角データ（２バイト・２の補数データ）
#aspが、30mx30mメッシュ間隔の方位角データ（２バイトデータ）
#この３つのデータで、大気地形効果補正処理を行う
set DEM=/usr1/DATA/IMAGE/GIS/DEM/iwa_utm              #標高データファイル名
set SLP=/usr1/DATA/IMAGE/GIS/DEM/iwa_slp              #傾斜角データファイル名
set ASP=/usr1/DATA/IMAGE/GIS/DEM/iwa_asp              #方位角データファイル名

########### Vegetation Map directory
#植生図のファイル名。絶対パス表示、img拡張子省略
#３は第３次植生調査（１９８７年）、４は第４次植生調査
#植生項目は全部で111項目あり、VEG/iwate3/item.docに書いてある

set VEG=/usr1/DATA/IMAGE/GIS/VEG/iwaveg7     #３次と４次調査合成植生図ファイル名を７種にマージ

########### 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             #上記850616A、1985年６月の６を記入
set date=16           #１６日
set time=0.745        #850616A 0:44:43で計算。44/60+43/3600=.745
#
########### Coodinate Asign
#調査場所の座標値x,yを記入する
#これらの記録は、最後に$area_name.docに保存される
#
set area_name="ichinoseki"     #地名を記入
set x0=622                     #4350x6800ＴＭ画像で、左上の座標値xを調べる
set y0=5422                    #左上のy座標値
set PIXEL=720                  #横幅サイズ
set LINE=620                   #縦幅サイズ
@ x1 = $x0 + $PIXEL - 1        #右下座標値xの計算
@ y1 = $y0 + $LINE - 1         #右下座標値yの計算

############ $area_name.doc log save
if ( -e $area_name.doc ) then      #現在のディレクトリにReadme.docがあるかどうか検査
  rm $area_name.doc                #もし存在するのなら古いReadme.docなので抹消
endif                          #if文の終わり

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

########### subset oparation
#画像切り取り作業の開始
#veg.imgは植生図データ、pre.bsqはＴＭ画像BSQ形式ファイル、pre.imgはＴＭ画像データ
#
   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 /usr1/DATA/IMAGE/GIS/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       : 6' >> 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  | /usr1/DATA/IMAGE/GIS/6s/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 'Modified Cosine Coefficient Table' >> $area_name.doc         #以上の係数項を$area_name.docに記録
echo '------------------------------------' >> $area_name.doc

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

set numb=0                              #ここからループシェルが始まる。

set cn1=0                               #cn1という変数を０にセットする。
while ( $cn1 < 3 )                      #cn1の値が３未満までループを続ける。以下同様 

set cn2=0
while ( $cn2 < 3 )

set cn3=0
while ( $cn3 < 3 )

set file_name="$area_name$cn1$cn2$cn3"
############# Modified Ccosine Method
#修正コサイン法による大気地形効果の補正計算
#あたかも真上から太陽光が入射したかのような結果のデータを得るのが目的
#これらの係数は、統計解析によって求めたが、正確ではない
#左からバンド１、２、３・・・７の修正コサイン法係数になっている

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

   @ k1 = 59 + $cn1 * 2                                            #59にcn1*2を足した値をk1とする。
   @ k2 = 20 + $cn2 * 2                                            #同様
   @ k3 = 16 + $cn3 * 1　　　　　　　　　　　　　　　　　　　　　　#計算式の前には必ず＠マークをつける。　

echo -n $k1 " " $k2 " " $k3 " "  > pre.crt　　　　　　　　　　　　 #-nとは行替えをしないで次の出力は右側について来る。

#
#echo '63.32  21.94  16.66  51.73  28.87  8.10' > pre.crt          #論文で示す切片aのこと

echo ' 51.73  28.87  8.10' >> pre.crt                              #Change
echo '-0.003 -0.001  -0.001  0  0  0' >> pre.crt                   #論文で示す傾斜b
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.01392 0.00840 0.00897 0.00220 0.00016 0.00009' >> pre.crt  #参考係数（入力時必要だが無視）

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

cat pre.crt >> $file_name.doc
echo '------------------------------------' >> $file_name.doc


############ Correct operation
#実際の修正コサイン法による大気地形効果補正
#これで出力されるpost.bsqが、補正された結果であり、pre.bsqに比べ
#影の部分が平坦になったり、標高差による輝度の違いも少なくなっている
#
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部分;ピラミッドリンキングでフィルタがかかった画像を入力する、という意味
#部分;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 16 10 50 100 18 24' >> $file_name.doc

#echo 'isodata c g 30 10 50 100 18 24' >> $file_name.doc

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

cp g*.cod g.cod


#分類作業の記録
echo 'Unsupervised Clasification 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.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
#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

cp g.cod g.rst

#rm *.cod

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

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

############## カテゴリー対応作業
echo ""  >> $file_name.doc
echo "作成されたay.pgmによりカテゴリーを自動的に割当てる"   >> $file_name.doc
echo ""  >> $file_name.doc

############## 植生図の凡例表示
echo "黒(0)未分類、黄(1)草地、茶(2)針葉、緑(3)広葉、水(4)水田、白(5)都市、青(6)水域" >> $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
cat ay.cls >>  $file_name.doc

reclass ay ay ayr

rm ay.tp1 ay.cls

echo 'カテゴリー対応による再配置画像をayr.imgとし'   >> $file_name.doc

########### カテゴリー対応による色割り当て結果
paint $PIXEL $LINE ayr.img /usr1/DATA/IMAGE/GIS/VEG/iwaveg7.tbl ayr.ppm
cjpeg ayr.ppm  > pyr.jpg
rm ayr.ppm


echo "" >> $file_name.doc
echo "" >> $file_name.doc


########### Confuse 評価の結果を$file_name.docに記録

echo 'CONFUSE マトリクスによる一致度の評価' >> $file_name.doc
confuse9 veg ayr | erchk 7 >> $file_name.doc



 @ cn3 = $cn3 + 1
 end

 @ cn2 = $cn2 + 1
 end

 @ cn1 = $cn1 + 1
 end


  # End of C-Shell Script
  # suguru kikuchi , Ichinoseki National College, 1999.1



date >> $file_name.doc                #全作業終了時間