#!/bin/bash

# Code to interogate the utm grid from Pat at CSUMB.
# Kurt Schwehr, Sept 2003

if [ 0 == 1 ]; then
    xmin_utm2=469701.080
    xmax_utm2=476294.460
    ymin_utm2=3634847.160
    ymax_utm2=3645222.660
    zmin_utm2=1.927
    zmax_utm2=322.301
    numpoints2=7922433
fi

declare -r in=Keep/LaJolla2mxyz.txt
#declare -r in=test.txt
tmpval=${in%%.txt}
declare -r name=${tmpval##*/}
unset tmpval

echo $name


# These take forever to make... only do it if we must
if [ ! -f x ]; then echo x; cat ${in} | awk '{print $1}' > x; fi
if [ ! -f y ]; then echo y; cat ${in} | awk '{print $2}' > y; fi
if [ ! -f z ]; then echo z; cat ${in} | awk '{print $3}' > z; fi

if [ ! -f x.sort ]; then echo x.sort; sort -n x > x.sort; fi
if [ ! -f y.sort ]; then echo y.sort; sort -n y > y.sort; fi
if [ ! -f z.sort ]; then echo z.sort; sort -n z > z.sort; fi

xmin_utm=`head -1 x.sort`
xmax_utm=`tail -1 x.sort`

ymin_utm=`head -1 y.sort`
ymax_utm=`tail -1 y.sort`

zmin_utm=`head -1 z.sort`
zmax_utm=`tail -1 z.sort`

if [ -z $numpoints ]; then
    echo finding numpoints
    numpoints=`wc -l z`
fi

if [ 0 == 1 ]; then
    echo xmin_utm=$xmin_utm
    echo xmax_utm=$xmax_utm
    echo ymin_utm=$ymin_utm
    echo ymax_utm=$ymax_utm
    echo zmin_utm=${zmin_utm}
    echo ymin_utm=${ymin_utm}
    echo numpoints=$numpoints
fi

if [ 0 == 1 ]; then
    echo xmin_utm2=$xmin_utm2
    echo xmax_utm2=$xmax_utm2
    echo ymin_utm2=$ymin_utm2
    echo ymax_utm2=$ymax_utm2
    echo zmin_utm2=${zmin_utm2}
    echo ymin_utm2=${ymin_utm2}
    echo numpoints2=$numpoints2
    echo name=$name
fi

# upper-left = ul
# lower-right= lr
ll=`echo "$xmin_utm $ymin_utm" | invproj +proj=utm +zone=11`
lr=`echo "$xmax_utm $ymin_utm" | invproj +proj=utm +zone=11`

ur=`echo "$xmax_utm $ymax_utm" | invproj +proj=utm +zone=11`
ul=`echo "$xmin_utm $ymax_utm" | invproj +proj=utm +zone=11`

if [ 1 == 0 ]; then
    echo ll="${ll}"
    echo lr="${lr}"
    echo ur="${ur}"
    echo ul="${ul}"
fi

echo "${ll} 0." > bbox.ll-proj
echo "${lr} 0." >> bbox.ll-proj
echo "${ur} 0." >> bbox.ll-proj
echo "${ul} 0." >> bbox.ll-proj

./proj-ll-fix.py bbox.ll-proj bbox.ll

echo "./bbox-xyz.py bbox.ll > ${name}-1-bbox.iv"
./bbox-xyz.py bbox.ll > ${name}-1-bbox.iv


ll=`head -1 bbox.ll | tail -1`
lr=`head -2 bbox.ll | tail -1`
ur=`head -3 bbox.ll | tail -1`
ul=`head -4 bbox.ll | tail -1`

echo "${ll}"
echo "${lr}"
echo "${ur}"
echo "${ul}"


unset xmin xmax ymin ymax
xmin=`echo ${ll} | awk '{print $1}'`
xmin2=`echo ${ul} | awk '{print $1}'`
xmax=`echo ${ur} | awk '{print $1}'`
xmax2=`echo ${lr} | awk '{print $1}'`
ymin=`echo ${ll} | awk '{print $2}'`
ymin2=`echo ${lr} | awk '{print $2}'`
ymax=`echo ${ur} | awk '{print $2}'`
ymax2=`echo ${lr} | awk '{print $2}'`

if [ 1 == 0 ]; then
    echo Lon-Lat
    echo $xmin .. $xmax
    echo $ymin .. $ymax
fi

xmin=`./min.py $xmin $xmin2`
ymin=`./min.py $ymin $ymin2`
xmax=`./max.py $xmax $xmax2`
ymax=`./max.py $ymax $ymax2`

xmin=`printf "%.6f" $xmin`
ymin=`printf "%.3f" $ymin`
xmax=`printf "%.6f" $xmax`
ymax=`printf "%.3f" $ymax`

echo
echo Lon-Lat
echo $xmin .. $xmax
echo $ymin .. $ymax

xrange=`python -c "print $xmax - $xmin"`
yrange=`python -c "print $ymax - $ymin"`
#zrange=`python -c "print $zmax - $zmin"`

if [ 1 == 0 ]; then
    echo xrange=$xrange
    echo yrange=$yrange
fi

xrange=`printf "%.6f" ${xrange}`
yrange=`printf "%.6f" ${yrange}`

echo
echo xrange=$xrange
echo yrange=$yrange
#echo zrange=$zrange

xdelta=`python -c "print ${xrange}/400."`
ydelta=`python -c "print ${yrange}/400."`

#echo xdelta=${xdelta}
xdelta=`printf "%.8f" ${xdelta}`
ydelta=`printf "%.8f" ${ydelta}`
echo xdelta=${xdelta}
echo ydelta=${ydelta}





# Make a bounding box of the data set.

out=${name}-bbox.iv
echo "Creating $out"
cat <<EOF > $out
#Inventor V2.1 ascii
Separator {
        Coordinate3 { point [
                             $ll ,
                             $lr ,
                             $ur ,
                             $ul 
			    ] }
        #IndexedFaceSet { coordIndex [ 0, 1, 2, 3, -1, ]}
        BaseColor { rgb [ 1 0 0] }
        IndexedLineSet { coordIndex [ 0, 1, 2, 3, 0 ] }

	#Translation { translation $XAVE $YAVE 0.05 }
	#Scale { scaleFactor .002 .002 .002 }
	#BaseColor { rgb [ 1 1 1] }
	#Text3 { string "${file%%.grd}" }
}

EOF

if [ ! -f ${name}-ll.xyz ]; then
    echo "./utm2latlon.py ${in} ${name}-ll.xyz"
    ./utm2latlon.py ${in} ${name}-ll.xyz
    exit
fi



if [ ! -f ${name}.grd ]; then
  echo "xyz2grd ${name}-ll.xyz -G${name}.grd -I${xdelta}/${ydelta} -R${xmin}/${xmax}/${ymin}/${ymax} -V"
  xyz2grd ${name}-ll.xyz -G${name}.grd -I${xdelta}/${ydelta} -R${xmin}/${xmax}/${ymin}/${ymax} -V 
  exit
fi







if [ ! -f ${name}_2.grd ]; then
  xdelta2=`python -c "print ${xdelta}*2."`
  ydelta2=`python -c "print ${ydelta}*2."`
  grdsample ${name}.grd -G${name}_2.grd -I${xdelta2}/${ydelta2} -V
fi

if [ ! -f ${name}_4.grd ]; then
  xdelta4=`python -c "print ${xdelta}*4."`
  ydelta4=`python -c "print ${ydelta}*4."`
  grdsample ${name}.grd -G${name}_4.grd -I${xdelta4}/${ydelta4} -V
fi

if [ ! -f ${name}_8.grd ]; then
  xdelta8=`python -c "print ${xdelta}*8."`
  ydelta8=`python -c "print ${ydelta}*8."`
  grdsample ${name}.grd -G${name}_8.grd -I${xdelta8}/${ydelta8} -V
fi

if [ ! -f ${name}_16.grd ]; then
  xdelta16=`python -c "print ${xdelta}*16."`
  ydelta16=`python -c "print ${ydelta}*16."`
  grdsample ${name}.grd -G${name}_16.grd -I${xdelta16}/${ydelta16} -V
fi
