#!/bin/csh -f
##############################################
#　最適結果への自動調整
#
#1999年06月01
#地形データの統計量を算出
#zstat slp.img 500 500
#
#1999年05月29
# comp7自動対応 & 修正コサイン
#

set area_name="HIROTA2" 

@ x0 = 2500
@ y0 = 5500

set PIXEL=500
set LINE=500

@ x1 = $x0 + $PIXEL - 1
@ y1 = $y0 + $LINE - 1

#Best 231 for Average
set start1=60
set start2=20
set start3=18


########### 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/iwaveg72


set mon=6
set date=16
set time=0.745

############ 
echo Image Center Position Information > $area_name.doc  #タイトルを記録保存
echo "1985_0616A 0:44:43(GST) 40:20(N) 141:24(E) 61.3(SUN_ang) 118(SUN_Zen)" >> $area_name.doc
echo "1985_0616B 0:45:07(GST) 38:54(N) 140:58(E) 61.7(SUN_ang) 115.2(SUN_Zen)" >> $area_name.doc
echo '------------------------------------' >> $area_name.doc  #単なる区切り線
date >> $area_name.doc             #処理開始時間の保存
hostname >> $area_name.doc         #処理しているマシン名の保存
pwd >>  $area_name.doc             #現在位置（ディレクトリ）の記録
whoami >> $area_name.doc           #誰が作業したものか記録

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

   if ( -e veg.img ) then
      rm veg.img
   endif
   if ( -e pre.bsq ) then
      rm pre.bsq
   endif
   if ( -e pre.img ) then
      rm 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$area_name.jpg
rm veg.ppm

subset ~/VEG/original/iwate3+4th.img veg_org.img $y0 $y1 $x0 $x1
paint $PIXEL $LINE veg_org.img ~/VEG/paint5.tbl veg.ppm
cjpeg veg.ppm > veg_new.jpg
rm veg.ppm

echo '------------ Actual Vegetation Items ------------' >> $area_name.doc
freq veg_org | awk ' NR>8 {printf(" %d \n",$1)}' | legendc.bat > zz1
freq veg_org | awk ' NR>8 {print $0}' > zz2
paste zz2 zz1 >> $area_name.doc
echo '-------------------------------------------------' >> $area_name.doc
rm zz1 zz2

########### 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
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
endif

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

############# Sun angle,zenith value
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

zsubset $DEM dem $y0 $y1 $x0 $x1
zsubset $SLP slp $y0 $y1 $x0 $x1
zsubset $ASP asp $y0 $y1 $x0 $x1

############# Incident Data get
zincident slp asp inc $suna $sunz

echo "# DEM.img Statistical data ----------" >> $area_name.doc
zstat dem.img $PIXEL $LINE >> $area_name.doc
echo "# SLP.img Statistical data ----------" >> $area_name.doc
zstat slp.img $PIXEL $LINE >> $area_name.doc
echo "# ASP.img Statistical data ----------" >> $area_name.doc
zstat asp.img $PIXEL $LINE >> $area_name.doc
echo "# INC.img Statistical data ----------" >> $area_name.doc
zstat inc.img $PIXEL $LINE >> $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
echo '------------------------------------' >> $area_name.doc

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

@ k1 = $start1
@ k2 = $start2
@ k3 = $start3

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

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 'Modified Cosine Coefficient Table' >> $area_name.doc
cat pre.crt >> $area_name.doc
echo '------------------------------------' >> $area_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
mv post.bsq a.bsq    #ファイル名の変更。大気地形効果補正画像のBSQファイルを分類用に使う
mv post.doc a.doc    #付属ファイルも名前変更


############ ピラミッドリンキング
pcomp > out

############ ISODATA法によるクラスタ分類
isodata c g 40 10 50 100 18 24

#上の分類コマンドそのものを記録しておく
echo 'isodata c g 40 10 50 100 18 24' >> $area_name.doc
echo "" >> $area_name.doc
echo '------------------------------------' >> $area_name.doc

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

cp g*.cod g.cod

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


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


############# クラスタ分析にもとずいてクラスタによる画像作成・第２段階
prcont cy.img by.img b

############# クラスタ分析にもとずいてクラスタによる画像作成・第３段階
prcont by.img ay.img a

############# 後の作業で使用しない画像データ等の削除
gvc2 dem.img $PIXEL $LINE dem.pgm
cjpeg dem.pgm  > dem.jpg
gvc2 asp.img $PIXEL $LINE asp.pgm
cjpeg asp.pgm  > asp.jpg
gvc2 slp.img $PIXEL $LINE slp.pgm
cjpeg slp.pgm  > slp.jpg
gvc2 inc.img $PIXEL $LINE inc.pgm
cjpeg inc.pgm  > inc.jpg
rm *.pgm


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

comp71 $PIXEL $LINE veg.img ay.img 0

echo "comp7 Process Example" >> $area_name.doc
echo '------------------------------------' >> $area_name.doc
comp71 $PIXEL $LINE veg.img ay.img 1 >> $area_name.doc

set num=`freq ay | tail -n 1 | awk '{print $1}'`

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


reclass ay ay0 ayr0
reclass ay ay1 ayr1
reclass ay ay2 ayr2


########### カテゴリー対応による色割り当て結果
paint $PIXEL $LINE ayr0.img ~/VEG/iwaveg7.tbl ayr.ppm
cjpeg ayr.ppm  > Y0$area_name.jpg
paint $PIXEL $LINE ayr1.img ~/VEG/iwaveg7.tbl ayr.ppm
cjpeg ayr.ppm  > Y1$area_name.jpg
paint $PIXEL $LINE ayr2.img ~/VEG/iwaveg7.tbl ayr.ppm
cjpeg ayr.ppm  > Y2$area_name.jpg
rm ayr.ppm

echo '## Maximum Average Value ##' >> $area_name.doc
cat ay0.cls >> $area_name.doc
confuse9 veg ayr0 | erchk 7 >> $area_name.doc
echo '## Maximum Overall Value ##' >> $area_name.doc
cat ay1.cls >> $area_name.doc
confuse9 veg ayr1 | erchk 7 >> $area_name.doc
echo '## Maximum Kappa Value ##' >> $area_name.doc
cat ay2.cls >> $area_name.doc
confuse9 veg ayr2 | erchk 7 >> $area_name.doc

rm ayr?.img ayr?.doc



# End of C-Shell Script